! +-======-+ ! Copyright (c) 2003-2018 United States Government as represented by ! the Admistrator of the National Aeronautics and Space Administration. ! All Rights Reserved. ! ! THIS OPEN SOURCE AGREEMENT ("AGREEMENT") DEFINES THE RIGHTS OF USE, ! REPRODUCTION, DISTRIBUTION, MODIFICATION AND REDISTRIBUTION OF CERTAIN ! COMPUTER SOFTWARE ORIGINALLY RELEASED BY THE UNITED STATES GOVERNMENT AS ! REPRESENTED BY THE GOVERNMENT AGENCY LISTED BELOW ("GOVERNMENT AGENCY"). ! THE UNITED STATES GOVERNMENT, AS REPRESENTED BY GOVERNMENT AGENCY, IS AN ! INTENDED THIRD-PARTY BENEFICIARY OF ALL SUBSEQUENT DISTRIBUTIONS OR ! REDISTRIBUTIONS OF THE SUBJECT SOFTWARE. ANYONE WHO USES, REPRODUCES, ! DISTRIBUTES, MODIFIES OR REDISTRIBUTES THE SUBJECT SOFTWARE, AS DEFINED ! HEREIN, OR ANY PART THEREOF, IS, BY THAT ACTION, ACCEPTING IN FULL THE ! RESPONSIBILITIES AND OBLIGATIONS CONTAINED IN THIS AGREEMENT. ! ! Government Agency: National Aeronautics and Space Administration ! Government Agency Original Software Designation: GSC-15354-1 ! Government Agency Original Software Title: GEOS-5 GCM Modeling Software ! User Registration Requested. Please Visit http://opensource.gsfc.nasa.gov ! Government Agency Point of Contact for Original Software: ! Dale Hithon, SRA Assistant, (301) 286-2691 ! ! +-======-+ module shallow use GEOS_Mod, only: write_parallel use MAPL_Mod, only: MAPL_UNDEF use GEOS_UtilsMod, only: GEOS_QSAT, GEOS_DQSAT use SHLWPARAMS use MAPL_ConstantsMod, only: MAPL_TICE , MAPL_CP , & MAPL_GRAV , MAPL_ALHS , & MAPL_ALHL , MAPL_ALHF , & MAPL_RGAS , MAPL_H2OMW, & MAPL_AIRMW, MAPL_RVAP , & MAPL_PI, r8 => MAPL_R8 implicit none private public compute_uwshcu_inv real*8, parameter :: xlv = MAPL_ALHL ! Latent heat of vaporization real*8, parameter :: xlf = MAPL_ALHF ! Latent heat of fusion real*8, parameter :: xls = MAPL_ALHS ! Latent heat sublimation real*8, parameter :: cp = MAPL_CP ! Specific heat of dry air real*8, parameter :: zvir = 0.609_r8 ! r_H2O/r_air-1 real*8, parameter :: r = MAPL_RGAS ! Gas constant for dry air real*8, parameter :: g = MAPL_GRAV ! Acceleration due to gravity real*8, parameter :: ep2 = MAPL_H2OMW/MAPL_AIRMW ! Ratio of molc wgt H2O to dry air real*8, parameter :: p00 = 1e5_r8 ! Reference pressure real*8, parameter :: rovcp = MAPL_RGAS/MAPL_CP ! Gas constant over specific heat contains real*8 function exnerfn(pressure) real*8, intent(in) :: pressure ! in Pa exnerfn = (pressure/p00)**rovcp return end function exnerfn subroutine compute_uwshcu_inv(idim, k0, ncnst, dt,pifc0_inv, zifc0_inv, & exnifc0_inv, pmid0_inv, zmid0_inv, exnmid0_inv, dp0_inv, & u0_inv, v0_inv, qv0_inv, ql0_inv, & qi0_inv, th0_inv, kpbl_inv, & tke_inv, cush, umf_inv, & qvten_inv, qlten_inv, qiten_inv, thten_inv, & uten_inv, vten_inv, qrten_inv, qsten_inv, & cufrc_inv, qcu_inv, qlu_inv, qiu_inv, cbmf, & qc_inv, cnt_inv, cnb_inv, cin, plcl, plfc, pinv,prel,pbup, & wlcl, qtsrc, thlsrc, thvlsrc, tkeavg, cldtop, & wu_inv, qtu_inv, thlu_inv, thvu_inv, uu_inv, vu_inv, & fer_inv, fdr_inv, xc_inv, qldet_inv, qidet_inv, qlsub_inv, & qisub_inv, ndrop_inv, nice_inv, thlsrc_pert, tr0_inv,shlwparams) implicit none integer, intent(in) :: idim integer, intent(in) :: k0 integer, intent(in) :: ncnst real , intent(in) :: dt real, intent(in) :: pifc0_inv(idim,k0+1) ! Environmental pressure at the interfaces [ Pa ] real, intent(in) :: zifc0_inv(idim,k0+1) ! Environmental height at the interfaces [ m ] real, intent(in) :: exnifc0_inv(idim,k0+1) ! Exner function at the interfaces real, intent(in) :: pmid0_inv(idim,k0) ! Environmental pressure at the layer mid-point [ Pa ] real, intent(in) :: zmid0_inv(idim,k0) ! Environmental height at the layer mid-point [ m ] real, intent(in) :: exnmid0_inv(idim,k0) ! Exner function at the layer mid-point real, intent(in) :: dp0_inv(idim,k0) ! Environmental layer pressure thickness [ Pa ] > 0. ! real*8, intent(in) :: dpdry0_inv(idim,k0) ! Environmental dry layer pressure thickness [ Pa ] real, intent(in) :: u0_inv(idim,k0) ! Environmental zonal wind [ m/s ] real, intent(in) :: v0_inv(idim,k0) ! Environmental meridional wind [ m/s ] real, intent(in) :: qv0_inv(idim,k0) ! Environmental water vapor specific humidity [ kg/kg ] real, intent(in) :: ql0_inv(idim,k0) ! Environmental liquid water specific humidity [ kg/kg ] real, intent(in) :: qi0_inv(idim,k0) ! Environmental ice specific humidity [ kg/kg ] real, intent(in) :: th0_inv(idim,k0) ! Environmental temperature [ K ] real, intent(inout) :: tr0_inv(idim,k0,ncnst) ! Environmental tracers [ #, kg/kg ] real, intent(in) :: tke_inv(idim,k0+1) ! Turbulent kinetic energy at the interfaces [ m2/s2 ] ! at the previous time step [ fraction ] real, intent(in) :: kpbl_inv(idim) ! Height of PBL [ m ] real, intent(in) :: thlsrc_pert real, intent(inout) :: cush(idim) ! Convective scale height [m] real, intent(out) :: umf_inv(idim,k0+1) ! Updraft mass flux at interfaces [kg/m2/s] real, intent(out) :: qvten_inv(idim,k0) ! Tendency of water vapor specific humidity [ kg/kg/s ] real, intent(out) :: qlten_inv(idim,k0) ! Tendency of liquid water specific humidity [ kg/kg/s ] real, intent(out) :: qiten_inv(idim,k0) ! Tendency of ice specific humidity [ kg/kg/s ] real, intent(out) :: thten_inv(idim,k0) ! Tendency of potential temperature [ K/s ] real, intent(out) :: uten_inv(idim,k0) ! Tendency of zonal wind [ m/s2 ] real, intent(out) :: vten_inv(idim,k0) ! Tendency of meridional wind [ m/s2 ] ! real, intent(out) :: trten_inv(idim,k0,ncnst) ! Tendency of tracers [ #/s, kg/kg/s ] real, intent(out) :: qrten_inv(idim,k0) ! Tendency of rain water specific humidity [ kg/kg/s ] real, intent(out) :: qsten_inv(idim,k0) ! Tendency of snow specific humidity [ kg/kg/s ] real, intent(out) :: cufrc_inv(idim,k0) ! Shallow cumulus cloud fraction at the layer mid-point [ fraction ] real, intent(out) :: qcu_inv(idim,k0) ! Liquid+ice specific humidity within cumulus updraft [ kg/kg ] real, intent(out) :: qlu_inv(idim,k0) ! Liquid water specific humidity within cumulus updraft [ kg/kg ] real, intent(out) :: qiu_inv(idim,k0) ! Ice specific humidity within cumulus updraft [ kg/kg ] real, intent(out) :: qc_inv(idim,k0) ! Tendency of cumulus condensate detrained into the environment [ kg/kg/s ] real, intent(out) :: cbmf(idim) ! Cumulus base mass flux [ kg/m2/s ] real, intent(out) :: cnt_inv(idim) ! Cumulus top interface index, cnt = kpen [ no ] real, intent(out) :: cnb_inv(idim) ! Cumulus base interface index, cnb = krel - 1 [ no ] real, intent(out) :: cin(idim) real, intent(out) :: plcl(idim) real, intent(out) :: plfc(idim) real, intent(out) :: pinv(idim) real, intent(out) :: prel(idim) real, intent(out) :: pbup(idim) real, intent(out) :: wlcl(idim) real, intent(out) :: qtsrc(idim) real, intent(out) :: thlsrc(idim) real, intent(out) :: thvlsrc(idim) real, intent(out) :: tkeavg(idim) real, intent(out) :: cldtop(idim) real, intent(out) :: wu_inv(idim,k0+1) real, intent(out) :: qtu_inv(idim,k0+1) real, intent(out) :: thlu_inv(idim,k0+1) real, intent(out) :: thvu_inv(idim,k0+1) real, intent(out) :: uu_inv(idim,k0+1) real, intent(out) :: vu_inv(idim,k0+1) real, intent(out) :: fer_inv(idim,k0) real, intent(out) :: fdr_inv(idim,k0) real, intent(out) :: xc_inv(idim,k0) real, intent(out) :: qldet_inv(idim,k0) real, intent(out) :: qidet_inv(idim,k0) real, intent(out) :: qlsub_inv(idim,k0) real, intent(out) :: qisub_inv(idim,k0) real, intent(out) :: ndrop_inv(idim,k0) real, intent(out) :: nice_inv(idim,k0) type(shlwparam_type), intent(in) :: shlwparams !----- Internal variables ----- real*8 :: pifc0(idim,0:k0) ! Environmental pressure at the interfaces [ Pa ] real*8 :: zifc0(idim,0:k0) ! Environmental height at the interfaces [ m ] real*8 :: exnifc0(idim,0:k0) real*8 :: pmid0(idim,k0) ! Environmental pressure at the layer mid-point [ Pa ] real*8 :: zmid0(idim,k0) ! Environmental height at the layer mid-point [ m ] real*8 :: exnmid0(idim,k0) real*8 :: dp0(idim,k0) ! Environmental layer pressure thickness [ Pa ] > 0. ! real*8 :: dpdry0(idim,k0) ! Environmental dry layer pressure thickness [ Pa ] real*8 :: u0(idim,k0) ! Environmental zonal wind [ m/s ] real*8 :: v0(idim,k0) ! Environmental meridional wind [ m/s ] ! at the previous time step [ fraction ] real*8 :: qv0(idim,k0) ! Environmental water vapor specific humidity [ kg/kg ] real*8 :: ql0(idim,k0) ! Environmental liquid water specific humidity [ kg/kg ] real*8 :: qi0(idim,k0) ! Environmental ice specific humidity [ kg/kg ] real*8 :: th0(idim,k0) ! Environmental temperature [ K ] real*8 :: tke(idim,0:k0) real*8 :: tr0(idim,k0,ncnst) ! Environmental tracers [ #, kg/kg ] real*8 :: umf(idim,0:k0) ! Updraft mass flux at the interfaces [ kg/m2/s ] real*8 :: qvten(idim,k0) ! Tendency of water vapor specific humidity [ kg/kg/s ] real*8 :: qlten(idim,k0) ! Tendency of liquid water specific humidity [ kg/kg/s ] real*8 :: qiten(idim,k0) ! tendency of ice specific humidity [ kg/kg/s ] real*8 :: sten(idim,k0) ! Tendency of static energy [ J/kg/s ] real*8 :: uten(idim,k0) ! Tendency of zonal wind [ m/s2 ] real*8 :: vten(idim,k0) ! Tendency of meridional wind [ m/s2 ] real*8 :: trten(idim,k0,ncnst) ! Tendency of tracers [ #/s, kg/kg/s ] real*8 :: qrten(idim,k0) ! Tendency of rain water specific humidity [ kg/kg/s ] real*8 :: qsten(idim,k0) ! Tendency of snow speficif humidity [ kg/kg/s ] real*8 :: cufrc(idim,k0) ! Shallow cumulus cloud fraction at the layer mid-point [ fraction ] real*8 :: qcu(idim,k0) ! Condensate water specific humidity within cumulus updraft ! at the layer mid-point [ kg/kg ] real*8 :: qlu(idim,k0) ! Liquid water specific humidity within cumulus updraft ! at the layer mid-point [ kg/kg ] real*8 :: qiu(idim,k0) ! Ice specific humidity within cumulus updraft ! at the layer mid-point [ kg/kg ] real*8 :: qc(idim,k0) ! Tendency of cumulus condensate detrained into the environment [ kg/kg/s ] real :: cnt(idim) ! Cumulus top interface index, cnt = kpen [ no ] real :: cnb(idim) ! Cumulus base interface index, cnb = krel - 1 [ no ] real*8 :: wu(idim,0:k0) real*8 :: qtu(idim,0:k0) real*8 :: thlu(idim,0:k0) real*8 :: thvu(idim,0:k0) real*8 :: uu(idim,0:k0) real*8 :: vu(idim,0:k0) real*8 :: fer(idim,k0) real*8 :: fdr(idim,k0) real*8 :: xc(idim,k0) real*8 :: qldet(idim,k0) real*8 :: qidet(idim,k0) real*8 :: qlsub(idim,k0) real*8 :: qisub(idim,k0) real*8 :: ndrop(idim,k0) real*8 :: nice(idim,k0) real :: trten_inv(idim,k0,ncnst) ! Tendency of tracers [ #/s, kg/kg/s ] integer :: k ! Vertical index for local fields [ no ] integer :: k_inv ! Vertical index for incoming fields [ no ] integer :: m ! Tracer index [ no ] integer :: kpbl(idim) do k = 1, k0 k_inv = k0 + 1 - k pmid0(:idim,k) = pmid0_inv(:idim,k_inv) u0(:idim,k) = u0_inv(:idim,k_inv) v0(:idim,k) = v0_inv(:idim,k_inv) zmid0(:idim,k) = zmid0_inv(:idim,k_inv) exnmid0(:idim,k) = exnmid0_inv(:idim,k_inv) dp0(:idim,k) = dp0_inv(:idim,k_inv) qv0(:idim,k) = qv0_inv(:idim,k_inv) ql0(:idim,k) = ql0_inv(:idim,k_inv) qi0(:idim,k) = qi0_inv(:idim,k_inv) th0(:idim,k) = th0_inv(:idim,k_inv) do m = 1, ncnst tr0(:idim,k,m) = tr0_inv(:idim,k_inv,m) enddo enddo tke(:,:) = 0._r8 pifc0(:,:) = 0._r8 zifc0(:,:) = 0._r8 exnifc0(:,:) = 0._r8 do k = 0, k0-1 k_inv = k0 - k tke(:idim,k) = tke_inv(:idim,k_inv) pifc0(:idim,k) = pifc0_inv(:idim,k_inv) zifc0(:idim,k) = zifc0_inv(:idim,k_inv) exnifc0(:idim,k) = exnifc0_inv(:idim,k_inv) end do ! print *,'pifc0_inv=',pifc0_inv ! print *,'pifc0=',pifc0 pifc0(:idim,k0) = pmid0(:idim,k0)+(pmid0(:idim,k0)-pifc0(:idim,k0-1)) zifc0(:idim,k0) = zmid0(:idim,k0)+(zmid0(:idim,k0)-zifc0(:idim,k0-1)) exnifc0(:idim,k0) = exnifc0(:idim,k0-1) tke(:idim,k0) = tke(:idim,k0-1) kpbl = int(kpbl_inv) ! print *,'pifc0=',pifc0 ! print *,"pmid0=",pmid0 ! print *,"exnmid0=",exnmid0 call compute_uwshcu( idim,k0, dt, ncnst,pifc0, zifc0, & exnifc0, pmid0, zmid0, exnmid0, dp0, u0, v0, & qv0, ql0, qi0, th0, tr0, kpbl, tke, cush, umf, & qvten, qlten, qiten, sten, uten, vten, & qrten, qsten, cufrc, qcu, qlu, qiu, cbmf, qc, & cnt, cnb, cin, plcl, plfc, pinv, prel, pbup, & wlcl, qtsrc, thlsrc, thvlsrc, tkeavg, cldtop, & wu, qtu, thlu, thvu, uu, vu, fer, fdr, xc, qldet,& qidet, qlsub, qisub, ndrop, nice, thlsrc_pert, & trten, shlwparams ) ! Reverse again cnt_inv(:idim) = k0 + 1 - cnt(:idim) cnb_inv(:idim) = k0 + 1 - cnb(:idim) do k = 0, k0 k_inv = k0 + 1 - k umf_inv(:idim,k_inv) = umf(:idim,k) wu_inv(:idim,k_inv) = wu(:idim,k) qtu_inv(:idim,k_inv) = qtu(:idim,k) thlu_inv(:idim,k_inv) = thlu(:idim,k) thvu_inv(:idim,k_inv) = thvu(:idim,k) uu_inv(:idim,k_inv) = uu(:idim,k) vu_inv(:idim,k_inv) = vu(:idim,k) end do do k = 1, k0 k_inv = k0 + 1 - k qvten_inv(:idim,k_inv) = qvten(:idim,k) qlten_inv(:idim,k_inv) = qlten(:idim,k) qiten_inv(:idim,k_inv) = qiten(:idim,k) thten_inv(:idim,k_inv) = sten(:idim,k) / (cp*exnmid0(:idim,k)) uten_inv(:idim,k_inv) = uten(:idim,k) vten_inv(:idim,k_inv) = vten(:idim,k) qrten_inv(:idim,k_inv) = qrten(:idim,k) qsten_inv(:idim,k_inv) = qsten(:idim,k) cufrc_inv(:idim,k_inv) = cufrc(:idim,k) qcu_inv(:idim,k_inv) = qcu(:idim,k) qlu_inv(:idim,k_inv) = qlu(:idim,k) qiu_inv(:idim,k_inv) = qiu(:idim,k) qc_inv(:idim,k_inv) = qc(:idim,k) fer_inv(:idim,k_inv) = fer(:idim,k) fdr_inv(:idim,k_inv) = fdr(:idim,k) xc_inv(:idim,k_inv) = xc(:idim,k) qldet_inv(:idim,k_inv) = qldet(:idim,k) qidet_inv(:idim,k_inv) = qidet(:idim,k) qlsub_inv(:idim,k_inv) = qlsub(:idim,k) qisub_inv(:idim,k_inv) = qisub(:idim,k) ndrop_inv(:idim,k_inv) = ndrop(:idim,k) nice_inv(:idim,k_inv) = nice(:idim,k) do m = 1, ncnst trten_inv(:idim,k_inv,m) = trten(:idim,k,m) enddo end do end subroutine compute_uwshcu_inv subroutine compute_uwshcu(idim, k0, dt,ncnst, pifc0_in,zifc0_in,& exnifc0_in, pmid0_in, zmid0_in, exnmid0_in, dp0_in, & u0_in, v0_in, qv0_in, ql0_in, qi0_in, th0_in, tr0_inout, & kpbl_in, tke_in, cush_inout, umf_out, & qvten_out, qlten_out, & qiten_out, sten_out, uten_out, vten_out, qrten_out, & qsten_out, cufrc_out, & qcu_out, qlu_out, qiu_out, cbmf_out, qc_out, & cnt_out, cnb_out, cinh_out, plcl_out, plfc_out, pinv_out, & prel_out, pbup_out, wlcl_out, qtsrc_out, thlsrc_out, thvlsrc_out, & tkeavg_out, cldhgt_out, wu_out, qtu_out, thlu_out, & thvu_out, uu_out, vu_out, fer_out, fdr_out, xc_out, & qldet_out, qidet_out, qlsub_out, qisub_out, ndrop_out, & nice_out, thlsrc_pert, trten_out, shlwparams ) ! ------------------------------------------------------------ ! ! ! ! University of Washington Shallow Convection Scheme ! ! ! ! Described in Park and Bretherton. 2008. J. Climate : ! ! ! ! 'The University of Washington shallow convection and ! ! moist turbulent schemes and their impact on climate ! ! simulations with the Community Atmosphere Model' ! ! ! ! Coded in CESM by Sungsu Park. Oct.2005. ! ! May.2008. ! ! ! ! Coded in GEOS by Nathan Arnold. July 2016. ! ! ! ! ! ! For general questions, email sungsup@ucar.edu or ! ! sungsu@atmos.washington.edu ! ! ! ! For GEOS-specific questions, email nathan.arnold@nasa.gov ! ! ! ! ------------------------------------------------------------ ! integer, intent(in) :: idim ! Number of columns integer, intent(in) :: k0 ! Number of vertical levels integer, intent(in) :: ncnst ! Number of tracers real, intent(in) :: dt ! Timestep [s] real*8, intent(in) :: pifc0_in( idim,0:k0 ) ! Environmental pressure at interfaces [Pa] real*8, intent(in) :: zifc0_in( idim,0:k0 ) ! Environmental height at interfaces [m] real*8, intent(in) :: exnifc0_in( idim,0:k0 ) ! Exner function at interfaces real*8, intent(in) :: pmid0_in( idim,k0 ) ! Environmental pressure at midpoints [Pa] real*8, intent(in) :: zmid0_in( idim,k0 ) ! Environmental height at midpoints [m] real*8, intent(in) :: exnmid0_in( idim,k0 ) ! Exner function at midpoints real*8, intent(in) :: dp0_in( idim,k0 ) ! Environmental layer pressure thickness real*8, intent(in) :: u0_in ( idim,k0 ) ! Environmental zonal wind [m/s] real*8, intent(in) :: v0_in ( idim,k0 ) ! Environmental meridional wind [m/s] real*8, intent(in) :: qv0_in( idim,k0 ) ! Environmental specific humidity real*8, intent(in) :: ql0_in( idim,k0 ) ! Environmental liquid water specific humidity real*8, intent(in) :: qi0_in( idim,k0 ) ! Environmental ice specific humidity real*8, intent(in) :: th0_in ( idim,k0 ) ! Environmental potential temperature [K] real*8, intent(inout) :: tr0_inout(idim,k0,ncnst) ! Environmental tracers [ #, kg/kg ] real*8, intent(in) :: tke_in( idim,0:k0 ) ! Turbulent kinetic energy at interfaces integer, intent(in) :: kpbl_in( idim ) ! Boundary layer top layer index real, intent(in) :: thlsrc_pert real, intent(inout) :: cush_inout( idim ) ! Convective scale height [m] real*8, intent(out) :: umf_out(idim,0:k0) ! Updraft mass flux at the interfaces [ kg/m2/s ] real*8, intent(out) :: qvten_out(idim,k0) ! Tendency of water vapor specific humidity [ kg/kg/s ] real*8, intent(out) :: qlten_out(idim,k0) ! Tendency of liquid water specific humidity [ kg/kg/s ] real*8, intent(out) :: qiten_out(idim,k0) ! Tendency of ice specific humidity [ kg/kg/s ] real*8, intent(out) :: sten_out(idim,k0) ! Tendency of dry static energy [ J/kg/s ] real*8, intent(out) :: uten_out(idim,k0) ! Tendency of zonal wind [ m/s2 ] real*8, intent(out) :: vten_out(idim,k0) ! Tendency of meridional wind [ m/s2 ] real*8, intent(out) :: trten_out(idim,k0,ncnst) ! Tendency of tracers [ #/s, kg/kg/s ] real*8, intent(out) :: qrten_out(idim,k0) ! Tendency of rain water specific humidity [ kg/kg/s ] real*8, intent(out) :: qsten_out(idim,k0) ! Tendency of snow specific humidity [ kg/kg/s ] real*8, intent(out) :: wu_out(idim,0:k0) ! Updraft vertical velocity ! ( defined from the release level to 'kpen-1' interface ) real*8, intent(out) :: qtu_out(idim,0:k0) ! Updraft qt [ kg/kg ] real*8, intent(out) :: thlu_out(idim,0:k0) ! Updraft thl [ K ] real*8, intent(out) :: thvu_out(idim,0:k0) ! Updraft thv [ K ] real*8, intent(out) :: uu_out(idim,0:k0) ! Updraft zonal wind [ m/s ] real*8, intent(out) :: vu_out(idim,0:k0) ! Updraft meridional wind [ m/s ] real*8, intent(out) :: cufrc_out(idim,k0) ! Shallow cumulus cloud fraction at the layer mid-point [ fraction ] real*8, intent(out) :: qcu_out(idim,k0) ! Condensate water specific humidity within cumulus updraft [ kg/kg ] real*8, intent(out) :: qlu_out(idim,k0) ! Liquid water specific humidity within cumulus updraft [ kg/kg ] real*8, intent(out) :: qiu_out(idim,k0) ! Ice specific humidity within cumulus updraft [ kg/kg ] real, intent(out) :: cbmf_out(idim) ! Cloud base mass flux [kg/m2/s] real*8, intent(out) :: qc_out(idim,k0) ! Tendency of detrained cumulus condensate ! into the environment [ kg/kg/s ] real, intent(out) :: cnt_out(idim) ! Cumulus top interface index real, intent(out) :: cnb_out(idim) ! Cumulus base interface index real, intent(out) :: cinh_out(idim) real, intent(out) :: pinv_out(idim) ! PBL top pressure [ Pa ] real, intent(out) :: plfc_out(idim) ! LFC of source air [ Pa ] real, intent(out) :: plcl_out(idim) ! LCL of source air [ Pa ] real, intent(out) :: prel_out(idim) real, intent(out) :: pbup_out(idim) real*8, intent(out) :: qldet_out(idim,k0) real*8, intent(out) :: qidet_out(idim,k0) real*8, intent(out) :: qlsub_out(idim,k0) real*8, intent(out) :: qisub_out(idim,k0) real*8, intent(out) :: ndrop_out(idim,k0) real*8, intent(out) :: nice_out(idim,k0) real, intent(out) :: tkeavg_out(idim) ! Average tke over the PBL [ m2/s2 ] type(shlwparam_type), intent(in) :: shlwparams ! ! Internal Output Variables ! real cldhgt_out(idim) real*8 qtten_out(idim,k0) ! Tendency of qt [ kg/kg/s ] real*8 slten_out(idim,k0) ! Tendency of sl [ J/kg/s ] real*8 ufrc_out(idim,0:k0) ! Updraft fractional area at the interfaces [ fraction ] ! real*8 uflx_out(idim,0:k0) ! Updraft/pen.entrainment zonal momentum flux [ m/s/m2/s ] ! real*8 vflx_out(idim,0:k0) ! Updraft/pen.entrainment meridional momentum flux [ m/s/m2/s ] real*8 fer_out(idim,k0) ! Fractional lateral entrainment rate [ 1/Pa ] real*8 fdr_out(idim,k0) ! Fractional lateral detrainment rate [ 1/Pa ] real*8 xc_out(idim,k0) ! real*8 cinh_out(idim) ! Convective INhibition upto LFC (CIN) [ J/kg ] real*8 trflx_out(idim,0:k0,ncnst) ! Updraft/pen.entrainment tracer flux [ #/m2/s, kg/kg/m2/s ] !----------------------------------------------- ! One-dimensional variables at each grid point !----------------------------------------------- ! Input variables real*8 :: pifc0(0:k0) real*8 :: zifc0(0:k0) real*8 :: pmid0(k0) real*8 :: zmid0(k0) real*8 :: dp0(k0) real*8 :: u0(k0) real*8 :: v0(k0) real*8 :: tke(1:k0) real*8 :: qv0(k0) real*8 :: ql0(k0) real*8 :: qi0(k0) real*8 :: cush real*8 :: tr0(k0,ncnst) ! Environmental variables derived from input variables real*8 :: t0(k0) real*8 :: s0(k0) real*8 :: qt0(k0) real*8 :: thl0(k0) real*8 :: thvl0(k0) real*8 :: ssqt0(k0) real*8 :: ssthl0(k0) real*8 :: ssu0(k0) real*8 :: ssv0(k0) real*8 :: thv0bot(k0) real*8 :: thv0top(k0) real*8 :: thvl0bot(k0) real*8 :: thvl0top(k0) real*8 :: exnmid0(k0) real*8 :: exnifc0(0:k0) real*8 :: sstr0(k0,ncnst) ! 2-1. For preventing negative condensate at the provisional time step real*8 qv0_star(k0) ! Environmental water vapor specific humidity [ kg/kg ] real*8 ql0_star(k0) ! Environmental liquid water specific humidity [ kg/kg ] real*8 qi0_star(k0) ! Environmental ice specific humidity [ kg/kg ] real*8 s0_star(k0) ! Environmental dry static energy [ J/kg ] ! 3. Variables associated with cumulus convection real*8 umf(0:k0) ! Updraft mass real*8 emf(0:k0) ! Penetrative real*8 qvten(k0) ! Tendency of real*8 qlten(k0) ! Tendency of real*8 qiten(k0) ! Tendency of real*8 sten(k0) ! Tendency of real*8 uten(k0) ! Tendency of real*8 vten(k0) ! Tendency of real*8 qrten(k0) ! Tendency of real*8 qsten(k0) ! Tendency of real*8 slflx(0:k0) ! real*8 qtflx(0:k0) ! real*8 uflx(0:k0) ! real*8 vflx(0:k0) ! real*8 cufrc(k0) ! Shallow real*8 qcu(k0) ! Condensate real*8 qlu(k0) ! Liquid water real*8 qiu(k0) ! Ice specific real*8 dwten(k0) ! Detrained real*8 diten(k0) ! Detrained ice real*8 fer(k0) ! Fractional real*8 fdr(k0) ! Fractional real*8 xco(k0) real*8 uf(k0) ! Zonal wind at real*8 vf(k0) ! Meridional real*8 qc(k0) ! Tendency due real*8 qlten_det(k0) real*8 qiten_det(k0) real*8 qldet(k0) real*8 qidet(k0) real*8 qc_l(k0) ! Tendency due real*8 qc_i(k0) ! Tendency due real*8 qc_lm real*8 qc_im real*8 nc_lm real*8 nc_im real*8 totsink real*8 ql_emf_kbup real*8 qi_emf_kbup real*8 nl_emf_kbup real*8 ni_emf_kbup real*8 cnt ! Cumulus top real*8 cnb ! Cumulus base real*8 qtten(k0) ! Tendency of real*8 slten(k0) ! Tendency of real*8 ufrc(0:k0) ! Updraft real*8 trten(k0,ncnst) ! Tendency of real*8 trflx(0:k0,ncnst) ! Flux of real*8 trflx_d(0:k0) ! Adjustive real*8 trflx_u(0:k0) ! Adjustive real*8 trmin real*8 pdelx, dum !----- Variables used for the calculation of condensation sink associated with compensating subsidence ! In the current code, this 'sink' tendency is simply set to be zero. real*8 uemf(0:k0) ! Net updraft mass flux at the interface ( emf + umf ) [ kg/m2/s ] real*8 comsub(k0) ! Compensating subsidence ! at the layer mid-point ( unit of mass flux, umf ) [ kg/m2/s ] real*8 qlten_sink(k0) ! Liquid condensate tendency ! by compensating subsidence/upwelling [ kg/kg/s ] real*8 qiten_sink(k0) ! Ice condensate tendency ! by compensating subsidence/upwelling [ kg/kg/s ] real*8 thlten_sub, qtten_sub ! Tendency of conservative scalars ! by compensating subsidence/upwelling real*8 qlten_sub, qiten_sub ! Tendency of ql0, qi0 ! by compensating subsidence/upwelling real*8 nlten_sub, niten_sub ! Tendency of nl0, ni0 ! by compensating subsidence/upwelling real*8 thl_prog, qt_prog ! Prognosed 'thl, qt' ! by compensating subsidence/upwelling !----- Variables describing cumulus updraft real*8 wu(0:k0) ! Updraft vertical velocity at the interface [ m/s ] real*8 thlu(0:k0) ! Updraft liquid potential temperature at the interface [ K ] real*8 qtu(0:k0) ! Updraft total specific humidity at the interface [ kg/kg ] real*8 uu(0:k0) ! Updraft zonal wind at the interface [ m/s ] real*8 vu(0:k0) ! Updraft meridional wind at the interface [ m/s ] real*8 thvu(0:k0) ! Updraft virtual potential temperature at the interface [ m/s ] real*8 rei(k0) ! Updraft fractional mixing rate with the environment [ 1/Pa ] real*8 tru(0:k0,ncnst) ! Updraft tracers [ #, kg/kg ] !----- Variables describing conservative scalars of entraining downdrafts at the ! entraining interfaces, i.e., 'kbup <= k < kpen-1'. At the other interfaces, ! belows are simply set to equal to those of updraft for simplicity - but it ! does not influence numerical calculation. real*8 thlu_emf(0:k0) ! Penetrative downdraft liquid potential temperature ! at entraining interfaces [ K ] real*8 qtu_emf(0:k0) ! Penetrative downdraft total water ! at entraining interfaces [ kg/kg ] real*8 uu_emf(0:k0) ! Penetrative downdraft zonal wind ! at entraining interfaces [ m/s ] real*8 vu_emf(0:k0) ! Penetrative downdraft meridional wind ! at entraining interfaces [ m/s ] real*8 tru_emf(0:k0,ncnst) ! Penetrative Downdraft tracers ! at entraining interfaces [ #, kg/kg ] ! Other internal variables integer kk, k, i, kp1, km1, mm, m integer iter_scaleh, iter_xc integer id_check, status integer klcl ! Layer containing LCL of source air integer kinv ! Inversion layer with PBL top interface as lower interface integer krel ! Release layer where buoyancy sorting first occurs integer klfc ! LFC layer of cumulus source air integer kbup ! Top layer in which buoyancy is positive at top interface integer kpen ! Highest layer with positive updraft velocity logical id_exit logical forcedCu real*8 cin, cinlcl real*8 thlsrc, qtsrc, usrc, vsrc, thvlsrc real*8 uplus, vplus real*8 trsrc(ncnst), tre(ncnst) real*8 plcl, plfc, prel, wrel real*8 ee2, ud2, wtw, wtwb real*8 xc real*8 cldhgt, scaleh, tscaleh, cridis real*8 sigmaw, tkeavg, dpsum, dpi, thvlmin real*8 thlxsat, qtxsat, thvxsat, x_cu, x_en, thv_x0, thv_x1 real*8 dpe, exne, thvebot, thle, qte, ue, ve, thlue, qtue, wue real*8 mu, mumin0, mumin2, mulcl, mulclstar real*8 cbmf, wcrit, winv, wlcl, ufrcinv, ufrclcl real*8 exql, exqi, ppen real*8 thj, qvj, qlj, qij, thvj, tj, thv0j, rhomid0j, rhoifc0j, qse real*8 thl0top, thl0bot, qt0bot, qt0top, thvubot, thvutop real*8 thl0lcl, qt0lcl, thv0lcl, thv0rel, rho0inv, autodet real*8 thlu_top, qtu_top, qlu_top, qiu_top, qlu_mid, qiu_mid, exntop real*8 aquad, bquad, cquad, xc1, xc2, excessu, excess0, xsat, xs1, xs2 real*8 bogbot, bogtop, delbog, drage, expfac real*8 rcwp, rlwp, riwp, qcubelow, qlubelow, qiubelow real*8 qs real qsat_arg, qsat_pe real*8 pe real*8 xsrc, xmean, xtop, xbot, xflx(0:k0) !----- Some diagnostic internal output variables real*8 ufrcinvbase_out(idim) ! Cumulus updraft fraction at the PBL top [ fraction ] real*8 ufrclcl_out(idim) ! Cumulus updraft fraction at the LCL ! ( or PBL top when LCL is below PBL top ) [ fraction ] real*8 winvbase_out(idim) ! Cumulus updraft velocity at the PBL top [ m/s ] real wlcl_out(idim) ! Cumulus updraft velocity at the LCL ! ( or PBL top when LCL is below PBL top ) [ m/s ] ! real*8 pbup_out(idim) ! Highest interface level of positive buoyancy [ Pa ] real*8 ppen_out(idim) ! Highest interface evel where Cu w = 0 [ Pa ] real qtsrc_out(idim) ! Sourse air qt [ kg/kg ] real thlsrc_out(idim) ! Sourse air thl [ K ] real thvlsrc_out(idim) ! Sourse air thvl [ K ] real*8 emfkbup_out(idim) ! Penetrative downward mass flux at 'kbup' interface [ kg/m2/s ] real*8 cinlclh_out(idim) ! Convective INhibition upto LCL (CIN) [ J/kg = m2/s2 ] real*8 cbmflimit_out(idim) ! Cloud base mass flux limiter [ kg/m2/s ] real*8 zinv_out(idim) ! PBL top height [ m ] real*8 rcwp_out(idim) ! Layer mean Cumulus LWP+IWP [ kg/m2 ] real*8 rlwp_out(idim) ! Layer mean Cumulus LWP [ kg/m2 ] real*8 riwp_out(idim) ! Layer mean Cumulus IWP [ kg/m2 ] real*8 qtu_emf_out(idim,0:k0) ! Penetratively entrained qt [ kg/kg ] real*8 thlu_emf_out(idim,0:k0) ! Penetratively entrained thl [ K ] real*8 uu_emf_out(idim,0:k0) ! Penetratively entrained u [ m/s ] real*8 vu_emf_out(idim,0:k0) ! Penetratively entrained v [ m/s ] real*8 uemf_out(idim,0:k0) ! Net upward mass flux ! including penetrative entrainment (umf+emf) [ kg/m2/s ] real*8 dwten_out(idim,k0) real*8 diten_out(idim,k0) real*8 tru_out(idim,0:k0,ncnst) ! Updraft tracers [ #, kg/kg ] real*8 tru_emf_out(idim,0:k0,ncnst) ! Penetratively entrained tracers [ #, kg/kg ] real*8 wu_s(0:k0) ! Same as above but for implicit CIN real*8 qtu_s(0:k0) real*8 thlu_s(0:k0) real*8 thvu_s(0:k0) real*8 uu_s(0:k0) real*8 vu_s(0:k0) real*8 qtu_emf_s(0:k0) real*8 thlu_emf_s(0:k0) real*8 uu_emf_s(0:k0) real*8 vu_emf_s(0:k0) real*8 uemf_s(0:k0) real*8 tru_s(0:k0,ncnst) real*8 tru_emf_s(0:k0,ncnst) real*8 dwten_s(k0) real*8 diten_s(k0) real*8 excessu_arr_out(idim,k0) real*8 excessu_arr(k0) real*8 excessu_arr_s(k0) real*8 excess0_arr_out(idim,k0) real*8 excess0_arr(k0) real*8 excess0_arr_s(k0) real*8 xc_arr_out(idim,k0) real*8 xc_arr(k0) real*8 xc_arr_s(k0) real*8 aquad_arr_out(idim,k0) real*8 aquad_arr(k0) real*8 aquad_arr_s(k0) real*8 bquad_arr_out(idim,k0) real*8 bquad_arr(k0) real*8 bquad_arr_s(k0) real*8 cquad_arr_out(idim,k0) real*8 cquad_arr(k0) real*8 cquad_arr_s(k0) real*8 bogbot_arr_out(idim,k0) real*8 bogbot_arr(k0) real*8 bogbot_arr_s(k0) real*8 bogtop_arr_out(idim,k0) real*8 bogtop_arr(k0) real*8 bogtop_arr_s(k0) real*8 exit_ufrc(idim) real*8 exit_wtw(idim) real*8 exit_drycore(idim) real*8 exit_wu(idim) real*8 exit_cufilter(idim) real*8 exit_rei(idim) real*8 exit_kinv1(idim) real*8 exit_klfck0(idim) real*8 exit_klclk0(idim) real*8 exit_uwcu(idim) real*8 exit_conden(idim) real*8 limit_cinlcl(idim) real*8 limit_cin(idim) real*8 ind_delcin(idim) real*8 limit_rei(idim) real*8 limit_shcu(idim) real*8 limit_negcon(idim) real*8 limit_ufrc(idim) real*8 limit_ppen(idim) real*8 limit_emf(idim) real*8 limit_cbmf(idim) real*8 :: ufrcinvbase_s, ufrclcl_s, winvbase_s, wlcl_s, plcl_s, pinv_s, prel_s, plfc_s, & qtsrc_s, thlsrc_s, thvlsrc_s, emfkbup_s, cinlcl_s, pbup_s, ppen_s, cbmflimit_s, & tkeavg_s, zinv_s, rcwp_s, rlwp_s, riwp_s real*8 :: ufrcinvbase, winvbase, emfkbup, cbmflimit !----- Variables for implicit CIN computation real*8, dimension(k0) :: qv0_s , ql0_s , qi0_s , s0_s , u0_s , & v0_s , t0_s , & !qt0_s , thl0_s ,thvl0_s , & qvten_s , qlten_s, qiten_s , qrten_s , qsten_s , sten_s ,& uten_s, vten_s, cufrc_s, qcu_s, qlu_s, qiu_s, & fer_s, fdr_s, xc_s, qc_s, qtten_s, slten_s, qldet_s, & qidet_s, qlsub_s, qisub_s real*8, dimension(0:k0) :: umf_s , slflx_s , qtflx_s , ufrc_s , uflx_s , vflx_s real*8 :: cush_s , cin_s , cbmf_s, cnt_s, cnb_s real*8 :: cin_i, cin_f, del_CIN, ke, alpha !, thlj real*8 :: cinlcl_i, cinlcl_f, del_cinlcl integer :: iter real*8, dimension(k0,ncnst) :: tr0_s, trten_s real*8, dimension(0:k0,ncnst) :: trflx_s !----- Variables for temporary storage real*8, dimension(k0) :: qv0_o, ql0_o, qi0_o, t0_o, s0_o, u0_o, v0_o real*8, dimension(k0) :: qt0_o , thl0_o , thvl0_o , & thv0bot_o, thv0top_o, thvl0bot_o, thvl0top_o, & ssthl0_o , ssqt0_o , ssu0_o , ssv0_o real*8 :: tkeavg_o , thvlmin_o, qtsrc_o, thvlsrc_o, thlsrc_o , & usrc_o, vsrc_o, plcl_o, plfc_o, thv0lcl_o integer :: kinv_o , klcl_o , klfc_o real*8, dimension(k0,ncnst) :: tr0_o real*8, dimension(k0,ncnst) :: trten_o, sstr0_o real*8, dimension(0:k0,ncnst) :: trflx_o real*8, dimension(ncnst) :: trsrc_o integer :: ixnumliq, ixnumice, ixcldliq, ixcldice ! ------------------ ! ! ! ! Define Parameters ! ! ! ! ------------------ ! ! ------------------------ ! ! Iterative xc calculation ! ! ------------------------ ! integer :: niter_xc ! ----------------------------------------------------------- ! ! Choice of 'CIN = cin' (.true.) or 'CIN = cinlcl' (.false.). ! ! ----------------------------------------------------------- ! logical :: use_CINcin ! --------------------------------------------------------------- ! ! Choice of 'explicit' ( 1 ) or 'implicit' ( 2 ) CIN. ! ! ! ! When choose 'CIN = cinlcl' above, it is recommended not to use ! ! implicit CIN, i.e., do 'NOT' choose simultaneously : ! ! [ 'use_CINcin=.false. & 'iter_cin=2' ] ! ! since 'cinlcl' will be always set to zero whenever LCL is below ! ! the PBL top interface in the current code. So, averaging cinlcl ! ! of two iter_cin steps is likely not so good. Except that, all ! ! the other combinations of 'use_CINcin' & 'iter_cin' are OK. ! ! --------------------------------------------------------------- ! integer :: iter_cin ! ---------------------------------------------------------------- ! ! Choice of 'self-detrainment' by negative buoyancy in calculating ! ! cumulus updraft mass flux at the top interface in each layer. ! ! ---------------------------------------------------------------- ! logical :: use_self_detrain ! --------------------------------------------------------- ! ! Cumulus momentum flux : turn-on (.true.) or off (.false.) ! ! --------------------------------------------------------- ! logical :: use_momenflx ! ----------------------------------------------------------------------------------------- ! ! Penetrative Entrainment : Cumulative ( .true. , original ) or Non-Cumulative ( .false. ) ! ! This option ( .false. ) is designed to reduce the sensitivity to the vertical resolution. ! ! ----------------------------------------------------------------------------------------- ! logical :: use_cumpenent logical :: scverbose real*8 :: rpen ! penetrative entrainment efficiency ! ----------------------- ! ! For lateral entrainment ! ! ----------------------- ! real*8 :: rle ! For critical stopping distance for lateral entrainment [no unit] real*8 :: rkm ! Determine the amount of air that is involved in buoyancy-sorting [no unit] real*8 :: rkfre ! Vertical velocity variance as fraction of tke. real*8 :: rmaxfrac ! Maximum allowable 'core' updraft fraction real*8 :: mumin1 ! Normalized CIN ('mu') corresponding to 'rmaxfrac' at the PBL top ! obtaind by inverting 'rmaxfrac = 0.5*erfc(mumin1)'. ! [rmaxfrac:mumin1]=[ 0.05:1.163, 0.075:1.018, 0.1:0.906, 0.15:0.733, 0.2:0.595, 0.25:0.477] real*8 :: rbuoy ! For nonhydrostatic pressure effects on updraft [no unit] real*8 :: rdrag ! Drag coefficient [no unit] real*8 :: epsvarw ! Variance of w at PBL top by meso-scale component [m2/s2] real*8 :: PGFc ! This is used for calculating vertical variations cumulus ! 'u' & 'v' by horizontal PGF during upward motion [no unit] !!! TEMPORARY: should be ncnst array of minimum values for all constituents real*8, parameter,dimension(4) :: qmin = [0._r8,0._r8,0._r8,0._r8] ! ---------------------------------------- ! ! Bulk microphysics controlling parameters ! ! --------------------------------------------------------------------------- ! ! criqc : Maximum condensate that can be hold by cumulus updraft [kg/kg] ! ! frc_rasn : Fraction of precipitable condensate in the expelled cloud water ! ! from cumulus updraft. The remaining fraction ('1-frc_rasn') is ! ! 'suspended condensate'. ! ! 0 : all expelled condensate is 'suspended condensate' ! ! 1 : all expelled condensate is 'precipitable condensate' ! ! kevp : Evaporative efficiency ! ! noevap_krelkpen : No evaporation from 'krel' to 'kpen' layers ! ! --------------------------------------------------------------------------- ! real*8 :: criqc real*8 :: frc_rasn real*8 :: rdrop ! ------------------------ ! ! Assign parameter values ! ! ------------------------ ! niter_xc = shlwparams%niter_xc use_CINcin = shlwparams%use_CINcin iter_cin = shlwparams%iter_cin use_self_detrain = shlwparams%use_self_detrain use_momenflx = shlwparams%use_momenflx use_cumpenent = shlwparams%use_cumpenent scverbose = shlwparams%scverbose rpen = shlwparams%rpen rle = shlwparams%rle ! For critical stopping distance for lateral entrainment [no unit] rkm = shlwparams%rkm ! Determine the amount of air that is involved in buoyancy-sorting [no unit] rkfre = shlwparams%rkfre ! Vertical velocity variance as fraction of tke. rmaxfrac = shlwparams%rmaxfrac ! Maximum allowable 'core' updraft fraction mumin1 = shlwparams%mumin1 rbuoy = shlwparams%rbuoy ! For nonhydrostatic pressure effects on updraft [no unit] rdrag = shlwparams%rdrag ! Drag coefficient [no unit] epsvarw = shlwparams%epsvarw ! Variance of w at PBL top by meso-scale component [m2/s2] PGFc = shlwparams%PGFc ! This is used for calculating vertical variations cumulus ! 'u' & 'v' by horizontal PGF during upward motion [no unit] criqc = shlwparams%criqc frc_rasn = shlwparams%frc_rasn rdrop = shlwparams%rdrop !------------------------! ! ! ! Start Main Calculation ! ! ! !------------------------! !!! TEMPORARY test settings ixcldice = 1 ixcldliq = 2 ixnumliq = 3 ixnumice = 4 ! ------------------------------------------------------- ! ! Initialize output variables defined for all grid points ! ! ------------------------------------------------------- ! cldhgt_out(:idim) = 0.0_r8 umf_out(:idim,0:k0) = 0.0_r8 cufrc_out(:idim,:k0) = 0.0_r8 qcu_out(:idim,:k0) = 0.0_r8 qlu_out(:idim,:k0) = 0.0_r8 qiu_out(:idim,:k0) = 0.0_r8 fer_out(:idim,:k0) = 0.0_r8 fdr_out(:idim,:k0) = 0.0_r8 xc_out(:idim,:k0) = 0.0_r8 cinh_out(:idim) = MAPL_UNDEF cinlclh_out(:idim) = MAPL_UNDEF cbmf_out(:idim) = 0.0_r8 qc_out(:idim,:k0) = 0.0_r8 qldet_out(:idim,:k0) = 0.0_r8 qidet_out(:idim,:k0) = 0.0_r8 qlsub_out(:idim,:k0) = 0.0_r8 qisub_out(:idim,:k0) = 0.0_r8 ndrop_out(:idim,:k0) = 0.0_r8 nice_out(:idim,:k0) = 0.0_r8 cnt_out(:idim) = real(k0) cnb_out(:idim) = 0.0_r8 ufrc_out(:idim,0:k0) = 0.0_r8 ! uflx_out(:idim,0:k0) = 0.0_r8 ! vflx_out(:idim,0:k0) = 0.0_r8 ufrcinvbase_out(:idim) = 0.0_r8 ufrclcl_out(:idim) = 0.0_r8 winvbase_out(:idim) = 0.0_r8 wlcl_out(:idim) = 0.0_r8 plcl_out(:idim) = 0.0_r8 pinv_out(:idim) = 0.0_r8 plfc_out(:idim) = 0.0_r8 prel_out(:idim) = 0.0_r8 pbup_out(:idim) = 0.0_r8 ppen_out(:idim) = 0.0_r8 qtsrc_out(:idim) = 0.0_r8 thlsrc_out(:idim) = 0.0_r8 thvlsrc_out(:idim) = 0.0_r8 emfkbup_out(:idim) = 0.0_r8 cbmflimit_out(:idim) = 0.0_r8 tkeavg_out(:idim) = 0.0_r8 zinv_out(:idim) = 0.0_r8 rcwp_out(:idim) = 0.0_r8 rlwp_out(:idim) = 0.0_r8 riwp_out(:idim) = 0.0_r8 wu_out(:idim,0:k0) = MAPL_UNDEF qtu_out(:idim,0:k0) = MAPL_UNDEF thlu_out(:idim,0:k0) = MAPL_UNDEF thvu_out(:idim,0:k0) = MAPL_UNDEF uu_out(:idim,0:k0) = MAPL_UNDEF vu_out(:idim,0:k0) = MAPL_UNDEF qtu_emf_out(:idim,0:k0) = 0.0_r8 thlu_emf_out(:idim,0:k0) = 0.0_r8 uu_emf_out(:idim,0:k0) = 0.0_r8 vu_emf_out(:idim,0:k0) = 0.0_r8 uemf_out(:idim,0:k0) = 0.0_r8 dwten_out(:idim,:k0) = 0.0_r8 diten_out(:idim,:k0) = 0.0_r8 trten_out(:idim,:k0,:ncnst) = 0.0_r8 trflx_out(:idim,0:k0,:ncnst) = 0.0_r8 tru_out(:idim,0:k0,:ncnst) = 0.0_r8 tru_emf_out(:idim,0:k0,:ncnst) = 0.0_r8 excessu_arr_out(:idim,:k0) = 0.0_r8 excess0_arr_out(:idim,:k0) = 0.0_r8 xc_arr_out(:idim,:k0) = 0.0_r8 aquad_arr_out(:idim,:k0) = 0.0_r8 bquad_arr_out(:idim,:k0) = 0.0_r8 cquad_arr_out(:idim,:k0) = 0.0_r8 bogbot_arr_out(:idim,:k0) = 0.0_r8 bogtop_arr_out(:idim,:k0) = 0.0_r8 exit_UWCu(:idim) = 0.0_r8 exit_conden(:idim) = 0.0_r8 exit_klclk0(:idim) = 0.0_r8 exit_klfck0(:idim) = 0.0_r8 exit_ufrc(:idim) = 0.0_r8 exit_wtw(:idim) = 0.0_r8 exit_drycore(:idim) = 0.0_r8 exit_wu(:idim) = 0.0_r8 exit_cufilter(:idim) = 0.0_r8 exit_kinv1(:idim) = 0.0_r8 exit_rei(:idim) = 0.0_r8 limit_shcu(:idim) = 0.0_r8 limit_negcon(:idim) = 0.0_r8 limit_ufrc(:idim) = 0.0_r8 limit_ppen(:idim) = 0.0_r8 limit_emf(:idim) = 0.0_r8 limit_cinlcl(:idim) = 0.0_r8 limit_cin(:idim) = 0.0_r8 limit_cbmf(:idim) = 0.0_r8 limit_rei(:idim) = 0.0_r8 ind_delcin(:idim) = 0.0_r8 !======================== ! Start column loop !======================== do i = 1, idim id_exit = .false. pifc0(0:k0) = pifc0_in(i,0:k0) zifc0(0:k0) = zifc0_in(i,0:k0) pmid0(:k0) = pmid0_in(i,:k0) zmid0(:k0) = zmid0_in(i,:k0) dp0(:k0) = dp0_in(i,:k0) u0(:k0) = u0_in(i,:k0) v0(:k0) = v0_in(i,:k0) qv0(:k0) = qv0_in(i,:k0) ql0(:k0) = ql0_in(i,:k0) qi0(:k0) = qi0_in(i,:k0) tke(1:k0) = tke_in(i,1:k0) ! pblh = pblh_in(i) cush = cush_inout(i) do m = 1,ncnst ! loop over tracers tr0(:k0,m) = tr0_inout(i,:k0,m) end do !------------------------------------------------------! ! Compute basic thermodynamic variables directly from ! ! input variables for each column ! !------------------------------------------------------! ! Compute internal environmental variables exnmid0(:k0) = exnmid0_in(i,:k0) exnifc0(:k0) = exnifc0_in(i,:k0) t0(:k0) = th0_in(i,:k0) * exnmid0(:k0) s0(:k0) = g*zmid0(:k0) + cp*t0(:k0) qt0(:k0) = qv0(:k0) + ql0(:k0) + qi0(:k0) thl0(:k0) = ( t0(:k0) - xlv*ql0(:k0)/cp - xls*qi0(:k0)/cp ) / exnmid0(:k0) thvl0(:k0) = ( 1._r8 + zvir*qt0(:k0) )*thl0(:k0) ! Compute slopes of environmental variables in each layer ssthl0 = slope( k0, thl0, pmid0 ) ssqt0 = slope( k0, qt0 , pmid0 ) ssu0 = slope( k0, u0 , pmid0 ) ssv0 = slope( k0, v0 , pmid0 ) do m = 1,ncnst ! loop over tracers sstr0(:k0,m) = slope( k0, tr0(:k0,m), pmid0 ) end do ! Compute thv0 and thvl0 at top/bottom interfaces in each layer do k = 1,k0 thl0bot = thl0(k) + ssthl0(k)*(pifc0(k-1) - pmid0(k)) qt0bot = qt0(k) + ssqt0(k)*(pifc0(k-1) - pmid0(k)) call conden( pifc0(k-1),thl0bot,qt0bot,thj,qvj,qlj,qij,qse,id_check ) if ( id_check .eq. 1 ) then exit_conden(i) = 1.0 id_exit = .true. if (scverbose) then call write_parallel('------- UW ShCu: Exit, conden') end if go to 333 end if thv0bot(k) = thj*(1._r8 + zvir*qvj - qlj - qij) thvl0bot(k) = thl0bot*(1._r8 + zvir*qt0bot) thl0top = thl0(k) + ssthl0(k)*(pifc0(k) - pmid0(k)) qt0top = qt0(k) + ssqt0(k)*(pifc0(k) - pmid0(k)) if (k.lt.k0) then call conden( pifc0(k),thl0top,qt0top,thj,qvj,qlj,qij,qse,id_check ) if ( id_check .eq. 1 ) then exit_conden(i) = 1.0 id_exit = .true. if (scverbose) then call write_parallel('------- UW ShCu: Exit, conden') end if go to 333 end if thv0top(k) = thj*(1._r8 + zvir*qvj - qlj - qij) thvl0top(k) = thl0top*(1._r8 + zvir*qt0top) else thv0top(k) = thv0bot(k) thvl0top(k) = thvl0bot(k) end if end do ! print *,'thvl0bot=',thvl0bot ! print *,'thvl0top=',thvl0top ! ------------------------------------------------------------ ! ! Save input and related environmental thermodynamic variables ! ! for use at "iter_cin=2" when "del_CIN >= 0" ! ! ------------------------------------------------------------ ! qv0_o(:k0) = qv0(:k0) ql0_o(:k0) = ql0(:k0) qi0_o(:k0) = qi0(:k0) t0_o(:k0) = t0(:k0) s0_o(:k0) = s0(:k0) u0_o(:k0) = u0(:k0) v0_o(:k0) = v0(:k0) qt0_o(:k0) = qt0(:k0) thl0_o(:k0) = thl0(:k0) thvl0_o(:k0) = thvl0(:k0) ssthl0_o(:k0) = ssthl0(:k0) ssqt0_o(:k0) = ssqt0(:k0) thv0bot_o(:k0) = thv0bot(:k0) thv0top_o(:k0) = thv0top(:k0) thvl0bot_o(:k0) = thvl0bot(:k0) thvl0top_o(:k0) = thvl0top(:k0) ssu0_o(:k0) = ssu0(:k0) ssv0_o(:k0) = ssv0(:k0) do m = 1, ncnst tr0_o(:k0,m) = tr0(:k0,m) sstr0_o(:k0,m) = sstr0(:k0,m) enddo ! ---------------------------------------------- ! ! Initialize output variables at each grid point ! ! ---------------------------------------------- ! umf(0:k0) = 0.0_r8 emf(0:k0) = 0.0_r8 slflx(0:k0) = 0.0_r8 qtflx(0:k0) = 0.0_r8 uflx(0:k0) = 0.0_r8 vflx(0:k0) = 0.0_r8 qvten(:k0) = 0.0_r8 qlten(:k0) = 0.0_r8 qiten(:k0) = 0.0_r8 sten(:k0) = 0.0_r8 uten(:k0) = 0.0_r8 vten(:k0) = 0.0_r8 qrten(:k0) = 0.0_r8 qsten(:k0) = 0.0_r8 dwten(:k0) = 0.0_r8 diten(:k0) = 0.0_r8 cufrc(:k0) = 0.0_r8 qcu(:k0) = 0.0_r8 qlu(:k0) = 0.0_r8 qiu(:k0) = 0.0_r8 fer(:k0) = MAPL_UNDEF fdr(:k0) = MAPL_UNDEF xco(:k0) = 0.0_r8 cin = 0.0_r8 cinlcl = 0.0_r8 cbmf = 0.0_r8 qc(:k0) = 0.0_r8 qldet(:k0) = 0.0_r8 qidet(:k0) = 0.0_r8 qc_l(:k0) = 0.0_r8 qc_i(:k0) = 0.0_r8 cnt = real(k0) cnb = 0.0_r8 qtten(:k0) = 0.0_r8 slten(:k0) = 0.0_r8 ufrc(0:k0) = 0.0_r8 thlu(0:k0) = MAPL_UNDEF qtu(0:k0) = MAPL_UNDEF uu(0:k0) = MAPL_UNDEF vu(0:k0) = MAPL_UNDEF wu(0:k0) = MAPL_UNDEF thvu(0:k0) = MAPL_UNDEF thlu_emf(0:k0) = MAPL_UNDEF qtu_emf(0:k0) = MAPL_UNDEF uu_emf(0:k0) = MAPL_UNDEF vu_emf(0:k0) = MAPL_UNDEF ufrcinvbase = 0.0_r8 ufrclcl = 0.0_r8 winvbase = 0.0_r8 wlcl = 0.0_r8 emfkbup = 0.0_r8 cbmflimit = 0.0_r8 excessu_arr(:k0) = 0.0_r8 excess0_arr(:k0) = 0.0_r8 xc_arr(:k0) = 0.0_r8 aquad_arr(:k0) = 0.0_r8 bquad_arr(:k0) = 0.0_r8 cquad_arr(:k0) = 0.0_r8 bogbot_arr(:k0) = 0.0_r8 bogtop_arr(:k0) = 0.0_r8 uemf(0:k0) = 0.0_r8 comsub(:k0) = 0.0_r8 qlten_sink(:k0) = 0.0_r8 qiten_sink(:k0) = 0.0_r8 qlten_det(:k0) = 0.0_r8 qiten_det(:k0) = 0.0_r8 ! nlten_sink(:k0) = 0.0_r8 ! niten_sink(:k0) = 0.0_r8 do m = 1, ncnst trflx(0:k0,m) = 0.0_r8 trten(:k0,m) = 0.0_r8 tru(0:k0,m) = 0.0_r8 tru_emf(0:k0,m) = 0.0_r8 enddo !-----------------------------------------------! ! Below 'iter' loop is for implicit CIN closure ! !-----------------------------------------------! ! ----------------------------------------------------------------------------- ! ! It is important to note that this iterative cin loop is located at the outer ! ! shell of the code. Thus, source air properties can also be changed during the ! ! iterative cin calculation, because cumulus convection induces non-zero fluxes ! ! even at interfaces below PBL top height through 'fluxbelowinv' subroutine. ! ! ----------------------------------------------------------------------------- ! do iter = 1, iter_cin ! ---------------------------------------------------------------------- ! ! Cumulus scale height ! ! In contrast to the premitive code, cumulus scale height is iteratively ! ! calculated at each time step, and at each iterative cin step. ! ! It is not clear whether I should locate below two lines within or out ! ! of the iterative cin loop. ! ! ---------------------------------------------------------------------- ! tscaleh = cush cush = -1._r8 tkeavg = 0._r8 ! ----------------------------------------------------------------------- ! ! Find PBL top height interface index, 'kinv-1' where 'kinv' is the layer ! ! index with PBLH in it. When PBLH is exactly at interface, 'kinv' is the ! ! layer index having PBLH as a lower interface. ! ! In the previous code, I set the lower limit of 'kinv' by 2 in order to ! ! be consistent with the other parts of the code. However in the modified ! ! code, I allowed 'kinv' to be 1 & if 'kinv = 1', I just exit the program ! ! without performing cumulus convection. This new approach seems to be ! ! more reasonable: if PBL height is within 'kinv=1' layer, surface is STL ! ! interface (bflxs <= 0) and interface just above the surface should be ! ! either non-turbulent (Ri>0.19) or stably turbulent (0<=Ri<0.19 but this ! ! interface is identified as a base external interface of upperlying CL. ! ! Thus, when 'kinv=1', PBL scheme guarantees 'bflxs <= 0'. For this case ! ! it is reasonable to assume that cumulus convection does not happen. ! ! When these is SBCL, PBL height from the PBL scheme is likely to be very ! ! close at 'kinv-1' interface, but not exactly, since 'zi' information is ! ! changed between two model time steps. In order to ensure correct identi ! ! fication of 'kinv' for general case including SBCL, I imposed an offset ! ! of 5 [m] in the below 'kinv' finding block. ! ! ----------------------------------------------------------------------- ! ! do k = k0-1, 1, -1 ! if( (pblh + 5. - zifc0(k))*(pblh + 5. - zifc0(k+1)) .lt. 0. ) then ! kinv = k + 1 ! should this be k? ! go to 15 ! endif ! end do ! kinv = 1 if (kpbl_in(i).gt.50) then kinv = k0 - kpbl_in(i) + 1 else kinv = 5 end if ! print *,'kinv=',kinv 15 continue if( kinv .le. 1 ) then exit_kinv1(i) = 1. id_exit = .true. if (scverbose) then call write_parallel('------- UW ShCu: Exit, kinv<=1') end if go to 333 endif if( kinv .ge. k0/4 ) then id_exit = .true. if (scverbose) then call write_parallel('------- UW ShCu: Exit, kinv>k0/4') end if go to 333 endif ! From here, it must be 'kinv >= 2'. ! -------------------------------------------------------------------------- ! ! Find PBL averaged tke ('tkeavg') and minimum 'thvl' ('thvlmin') in the PBL ! ! In the current code, 'tkeavg' is obtained by averaging all interfacial TKE ! ! within the PBL. However, in order to be conceptually consistent with PBL ! ! scheme, 'tkeavg' should be calculated by considering surface buoyancy flux.! ! If surface buoyancy flux is positive ( bflxs >0 ), surface interfacial TKE ! ! should be included in calculating 'tkeavg', while if bflxs <= 0, surface ! ! interfacial TKE should not be included in calculating 'tkeavg'. I should ! ! modify the code when 'bflxs' is available as an input of cumulus scheme. ! ! 'thvlmin' is a minimum 'thvl' within PBL obtained by comparing top & base ! ! interface values of 'thvl' in each layers within the PBL. ! ! -------------------------------------------------------------------------- ! dpsum = 0._r8 thvlmin = 1000._r8 do k = 1,kinv ! max(kinv-1,1) ! Here, 'k' is an interfacial layer index. dpi = pifc0(k-1) - pifc0(k) dpsum = dpsum + dpi tkeavg = tkeavg + dpi*tke(k) if( k .ne. kinv ) thvlmin = min(thvlmin,min(thvl0bot(k),thvl0top(k))) end do tkeavg = tkeavg/dpsum !tkeavg = tkeavg/1e4 ! ------------------------------------------------------------------ ! ! Find characteristics of cumulus source air: qtsrc,thlsrc,usrc,vsrc ! ! Note that 'thlsrc' was concocted using 'thvlsrc' and 'qtsrc'. ! ! 'qtsrc' is defined as the lowest layer mid-point value; 'thlsrc' ! ! is from 'qtsrc' and 'thvlmin=thvlsrc'; 'usrc' & 'vsrc' are defined ! ! as the values just below the PBL top interface. ! ! ------------------------------------------------------------------ ! qtsrc = qt0(1) thvlsrc = thvlmin+thlsrc_pert thlsrc = thvlsrc / ( 1._r8 + zvir * qtsrc ) usrc = u0(kinv-1) + ssu0(kinv-1) * ( pifc0(kinv-1) - pmid0(kinv-1) ) vsrc = v0(kinv-1) + ssv0(kinv-1) * ( pifc0(kinv-1) - pmid0(kinv-1) ) do m = 1, ncnst trsrc(m) = tr0(1,m) enddo ! ------------------------------------------------------------------ ! ! Find LCL of the source air and a layer index containing LCL (klcl) ! ! When the LCL is exactly at the interface, 'klcl' is a layer index ! ! having 'plcl' as the lower interface similar to the 'kinv' case. ! ! In the previous code, I assumed that if LCL is located within the ! ! lowest model layer ( 1 ) or the top model layer ( k0 ), then no ! ! convective adjustment is performed and just exited. However, in ! ! the revised code, I relaxed the first constraint and even though ! ! LCL is at the lowest model layer, I allowed cumulus convection to ! ! be initiated. For this case, cumulus convection should be started ! ! from the PBL top height, as shown in the following code. ! ! When source air is already saturated even at the surface, klcl is ! ! set to 1. ! ! ------------------------------------------------------------------ ! if ( pifc0(0).lt.70000.or.pifc0(0).gt.115000. ) then if (scverbose) then print *,"pifc(0) outside valid range! pifc=",pifc0(0) end if id_exit = .true. go to 333 end if if (qtsrc.gt.0.1.or.qtsrc.lt.1e-8 ) then print *,"qtsrc outside valid range! qtsrc=",qtsrc id_exit = .true. go to 333 end if if (thlsrc.gt.400..or.thlsrc.lt.100. ) then print *,"thlsrc outside valid range! thlsrc=",thlsrc print *,"pifc0(0)=",pifc0(0) id_exit = .true. go to 333 end if plcl = qsinvert(qtsrc,thlsrc,pifc0(0)) do k = 0, k0 if( pifc0(k) .lt. plcl ) then klcl = k go to 25 endif end do klcl = k0 25 continue klcl = max(1,klcl) if( plcl .lt. 60000. ) then exit_klclk0(i) = 1. id_exit = .true. if (scverbose) then call write_parallel('------- UW ShCu: exit, plcl<600mb') end if go to 333 endif ! ------------------------------------------------------------- ! ! Calculate environmental virtual potential temperature at LCL, ! !'thv0lcl' which is solely used in the 'cin' calculation. Note ! ! that 'thv0lcl' is calculated first by calculating 'thl0lcl' ! ! and 'qt0lcl' at the LCL, and performing 'conden' afterward, ! ! in fully consistent with the other parts of the code. ! ! ------------------------------------------------------------- ! thl0lcl = thl0(klcl) + ssthl0(klcl) * ( plcl - pmid0(klcl) ) qt0lcl = qt0(klcl) + ssqt0(klcl) * ( plcl - pmid0(klcl) ) call conden(plcl,thl0lcl,qt0lcl,thj,qvj,qlj,qij,qse,id_check) if( id_check .eq. 1 ) then exit_conden(i) = 1. id_exit = .true. go to 333 end if thv0lcl = thj * ( 1._r8 + zvir * qvj - qlj - qij ) ! ------------------------------------------------------------------------ ! ! Compute Convective Inhibition, 'cin' & 'cinlcl' [J/kg]=[m2/s2] TKE unit. ! ! ! ! 'cin' (cinlcl) is computed from the PBL top interface to LFC (LCL) using ! ! piecewisely reconstructed environmental profiles, assuming environmental ! ! buoyancy profile within each layer ( or from LCL to upper interface in ! ! each layer ) is simply a linear profile. For the purpose of cin (cinlcl) ! ! calculation, we simply assume that lateral entrainment does not occur in ! ! updrafting cumulus plume, i.e., cumulus source air property is conserved.! ! Below explains some rules used in the calculations of cin (cinlcl). In ! ! general, both 'cin' and 'cinlcl' are calculated from a PBL top interface ! ! to LCL and LFC, respectively : ! ! 1. If LCL is lower than the PBL height, cinlcl = 0 and cin is calculated ! ! from PBL height to LFC. ! ! 2. If LCL is higher than PBL height, 'cinlcl' is calculated by summing ! ! both positive and negative cloud buoyancy up to LCL using 'single_cin'! ! From the LCL to LFC, however, only negative cloud buoyancy is counted ! ! to calculate final 'cin' upto LFC. ! ! 3. If either 'cin' or 'cinlcl' is negative, they are set to be zero. ! ! In the below code, 'klfc' is the layer index containing 'LFC' similarto ! ! 'kinv' and 'klcl'. ! ! ------------------------------------------------------------------------! cin = 0._r8 cinlcl = 0._r8 plfc = 0._r8 klfc = k0 !------------------------------------------------------------------------! ! Case 1. LCL height is higher than PBL interface ( 'pLCL <=ps0(kinv-1)' ) ! !------------------------------------------------------------------------! if( klcl .ge. kinv ) then do k = kinv, k0 - 1 if( k .lt. klcl ) then thvubot = thvlsrc thvutop = thvlsrc cin = cin + single_cin(pifc0(k-1),thv0bot(k),pifc0(k),thv0top(k),thvubot,thvutop) elseif( k .eq. klcl ) then !----- Bottom to LCL thvubot = thvlsrc thvutop = thvlsrc cin = cin + single_cin(pifc0(k-1),thv0bot(k),plcl,thv0lcl,thvubot,thvutop) if( cin .lt. 0._r8 ) limit_cinlcl(i) = 1. cinlcl = max(cin,0._r8) cin = cinlcl !----- LCL to Top thvubot = thvlsrc call conden(pifc0(k),thlsrc,qtsrc,thj,qvj,qlj,qij,qse,id_check) if( id_check .eq. 1 ) then exit_conden(i) = 1. id_exit = .true. go to 333 end if thvutop = thj * ( 1._r8 + zvir*qvj - qlj - qij ) call getbuoy(plcl,thv0lcl,pifc0(k),thv0top(k),thvubot,thvutop,plfc,cin) if( plfc .gt. 0._r8 ) then klfc = k go to 35 end if else thvubot = thvutop call conden(pifc0(k),thlsrc,qtsrc,thj,qvj,qlj,qij,qse,id_check) if( id_check .eq. 1 ) then exit_conden(i) = 1. id_exit = .true. go to 333 end if thvutop = thj * ( 1._r8 + zvir*qvj - qlj - qij ) call getbuoy(pifc0(k-1),thv0bot(k),pifc0(k),thv0top(k),thvubot,thvutop,plfc,cin) if( plfc .gt. 0._r8 ) then klfc = k go to 35 end if endif end do ! -----------------------------------------------------------------------! ! Case 2. LCL height is lower than PBL interface ( 'pLCL > ps0(kinv-1)') ! ! -----------------------------------------------------------------------! else cinlcl = 0._r8 do k = kinv, k0 - 1 call conden(pifc0(k-1),thlsrc,qtsrc,thj,qvj,qlj,qij,qse,id_check) if( id_check .eq. 1 ) then exit_conden(i) = 1. id_exit = .true. go to 333 end if thvubot = thj * ( 1._r8 + zvir*qvj - qlj - qij ) call conden(pifc0(k),thlsrc,qtsrc,thj,qvj,qlj,qij,qse,id_check) if( id_check .eq. 1 ) then exit_conden(i) = 1. id_exit = .true. go to 333 end if thvutop = thj * ( 1._r8 + zvir*qvj - qlj - qij ) call getbuoy(pifc0(k-1),thv0bot(k),pifc0(k),thv0top(k),thvubot,thvutop,plfc,cin) if( plfc .gt. 0._r8 ) then klfc = k go to 35 end if end do endif ! End of CIN case selection 35 continue if( cin .lt. 0._r8 ) limit_cin(i) = 1. cin = max(0._r8,cin) if( klfc .ge. k0 ) then klfc = k0 if (scverbose) then call write_parallel('------ UWShCu: klfc >= k0') end if exit_klfck0(i) = 1. id_exit = .true. go to 333 endif !----------------------------------------------------------------------! ! In order to calculate implicit 'cin' (or 'cinlcl'), save the initially ! ! calculated 'cin' and 'cinlcl', and other related variables. These will ! ! be restored after calculating implicit CIN. ! !----------------------------------------------------------------------! if( iter .eq. 1 ) then cin_i = cin cinlcl_i = cinlcl ke = rbuoy / ( rkfre * tkeavg + epsvarw ) kinv_o = kinv klcl_o = klcl klfc_o = klfc plcl_o = plcl plfc_o = plfc tkeavg_o = tkeavg thvlmin_o = thvlmin qtsrc_o = qtsrc thvlsrc_o = thvlsrc thlsrc_o = thlsrc usrc_o = usrc vsrc_o = vsrc thv0lcl_o = thv0lcl do m = 1, ncnst trsrc_o(m) = trsrc(m) enddo endif ! Modification : If I impose w = max(0.1, w) up to the top interface of ! klfc, I should only use cinlfc. That is, if I want to ! use cinlcl, I should not impose w = max(0.1, w). ! Using cinlcl is equivalent to treating only 'saturated' ! moist convection. Note that in this sense, I should keep ! the functionality of both cinlfc and cinlcl. ! However, the treatment of penetrative entrainment level becomes ! ambiguous if I choose 'cinlcl'. Thus, the best option is to use ! 'cinlfc'. ! -------------------------------------------------------------------------- ! ! Calculate implicit 'cin' by averaging initial and final cins. Note that ! ! implicit CIN is adopted only when cumulus convection stabilized the system,! ! i.e., only when 'del_CIN >0'. If 'del_CIN<=0', just use explicit CIN. Note ! ! also that since 'cinlcl' is set to zero whenever LCL is below the PBL top, ! ! (see above CIN calculation part), the use of 'implicit CIN=cinlcl' is not ! ! good. Thus, when using implicit CIN, always try to only use 'implicit CIN= ! ! cin', not 'implicit CIN=cinlcl'. However, both 'CIN=cin' and 'CIN=cinlcl' ! ! are good when using explicit CIN. ! ! -------------------------------------------------------------------------- ! if( iter .ne. 1 ) then cin_f = cin cinlcl_f = cinlcl if( use_CINcin ) then del_CIN = cin_f - cin_i else del_CIN = cinlcl_f - cinlcl_i endif if( del_CIN .gt. 0. ) then ! -------------------------------------------------------------- ! ! Calculate implicit 'cin' and 'cinlcl'. Note that when we chose ! ! to use 'implicit CIN = cin', choose 'cinlcl = cinlcl_i' below: ! ! because iterative CIN only aims to obtain implicit CIN, once ! ! we obtained 'implicit CIN=cin', it is good to use the original ! ! profiles information for all the other variables after that. ! ! Note 'cinlcl' will be explicitly used in calculating 'wlcl' & ! ! 'ufrclcl' after calculating 'winv' & 'ufrcinv' at the PBL top ! ! interface later, after calculating 'cbmf'. ! ! -------------------------------------------------------------- ! alpha = compute_alpha( del_CIN, ke ) cin = cin_i + alpha * del_CIN if( use_CINcin ) then cinlcl = cinlcl_i else cinlcl = cinlcl_i + alpha * del_cinlcl endif ! ----------------------------------------------------------------- ! ! Restore the original values from the previous 'iter_cin' step (1) ! ! to compute correct tendencies for (n+1) time step by implicit CIN ! ! ----------------------------------------------------------------- ! kinv = kinv_o klcl = klcl_o klfc = klfc_o plcl = plcl_o plfc = plfc_o tkeavg = tkeavg_o thvlmin = thvlmin_o qtsrc = qtsrc_o thvlsrc = thvlsrc_o thlsrc = thlsrc_o usrc = usrc_o vsrc = vsrc_o thv0lcl = thv0lcl_o do m = 1, ncnst trsrc(m) = trsrc_o(m) enddo qv0(:k0) = qv0_o(:k0) ql0(:k0) = ql0_o(:k0) qi0(:k0) = qi0_o(:k0) t0(:k0) = t0_o(:k0) s0(:k0) = s0_o(:k0) u0(:k0) = u0_o(:k0) v0(:k0) = v0_o(:k0) qt0(:k0) = qt0_o(:k0) thl0(:k0) = thl0_o(:k0) thvl0(:k0) = thvl0_o(:k0) ssthl0(:k0) = ssthl0_o(:k0) ssqt0(:k0) = ssqt0_o(:k0) thv0bot(:k0) = thv0bot_o(:k0) thv0top(:k0) = thv0top_o(:k0) thvl0bot(:k0) = thvl0bot_o(:k0) thvl0top(:k0) = thvl0top_o(:k0) ssu0(:k0) = ssu0_o(:k0) ssv0(:k0) = ssv0_o(:k0) do m = 1, ncnst tr0(:k0,m) = tr0_o(:k0,m) sstr0(:k0,m) = sstr0_o(:k0,m) enddo ! ------------------------------------------------------ ! ! Initialize all fluxes, tendencies, and other variables ! ! in association with cumulus convection. ! ! ------------------------------------------------------ ! umf(0:k0) = 0.0_r8 emf(0:k0) = 0.0_r8 slflx(0:k0) = 0.0_r8 qtflx(0:k0) = 0.0_r8 uflx(0:k0) = 0.0_r8 vflx(0:k0) = 0.0_r8 qvten(:k0) = 0.0_r8 qlten(:k0) = 0.0_r8 qiten(:k0) = 0.0_r8 sten(:k0) = 0.0_r8 uten(:k0) = 0.0_r8 vten(:k0) = 0.0_r8 qrten(:k0) = 0.0_r8 qsten(:k0) = 0.0_r8 dwten(:k0) = 0.0_r8 diten(:k0) = 0.0_r8 cufrc(:k0) = 0.0_r8 qcu(:k0) = 0.0_r8 qlu(:k0) = 0.0_r8 qiu(:k0) = 0.0_r8 fer(:k0) = 0.0_r8 fdr(:k0) = 0.0_r8 xco(:k0) = 0.0_r8 qc(:k0) = 0.0_r8 qldet(:k0) = 0.0_r8 qidet(:k0) = 0.0_r8 qlten_sub = 0.0_r8 qiten_sub = 0.0_r8 qc_l(:k0) = 0.0_r8 qc_i(:k0) = 0.0_r8 cbmf = 0.0_r8 cnt = k0 cnb = 0.0_r8 qtten(:k0) = 0.0_r8 slten(:k0) = 0.0_r8 ufrc(0:k0) = 0.0_r8 thlu(0:k0) = MAPL_UNDEF qtu(0:k0) = MAPL_UNDEF uu(0:k0) = MAPL_UNDEF vu(0:k0) = MAPL_UNDEF wu(0:k0) = MAPL_UNDEF thvu(0:k0) = MAPL_UNDEF thlu_emf(0:k0) = MAPL_UNDEF qtu_emf(0:k0) = MAPL_UNDEF uu_emf(0:k0) = MAPL_UNDEF vu_emf(0:k0) = MAPL_UNDEF do m = 1, ncnst trflx(0:k0,m) = 0.0_r8 trten(:k0,m) = 0.0_r8 tru(0:k0,m) = 0.0_r8 tru_emf(0:k0,m) = 0.0_r8 enddo ! -------------------------------------------------- ! ! Below are diagnostic output variables for detailed ! ! analysis of cumulus scheme. ! ! -------------------------------------------------- ! ufrcinvbase = 0.0_r8 ufrclcl = 0.0_r8 winvbase = 0.0_r8 wlcl = 0.0_r8 emfkbup = 0.0_r8 cbmflimit = 0.0_r8 excessu_arr(:k0) = 0.0_r8 excess0_arr(:k0) = 0.0_r8 xc_arr(:k0) = 0.0_r8 aquad_arr(:k0) = 0.0_r8 bquad_arr(:k0) = 0.0_r8 cquad_arr(:k0) = 0.0_r8 bogbot_arr(:k0) = 0.0_r8 bogtop_arr(:k0) = 0.0_r8 else ! When 'del_CIN < 0', use explicit CIN instead of implicit CIN. ! ----------------------------------------------------------- ! ! Identifier showing whether explicit or implicit CIN is used ! ! ----------------------------------------------------------- ! ind_delcin(i) = 1. if (scverbose) then call write_parallel('------ UWShCu: del_CIN<0') end if ! --------------------------------------------------------- ! ! Restore original output values of "iter_cin = 1" and exit ! ! --------------------------------------------------------- ! umf_out(i,0:k0) = umf_s(0:k0) qvten_out(i,:k0) = qvten_s(:k0) qlten_out(i,:k0) = qlten_s(:k0) qiten_out(i,:k0) = qiten_s(:k0) sten_out(i,:k0) = sten_s(:k0) uten_out(i,:k0) = uten_s(:k0) vten_out(i,:k0) = vten_s(:k0) qrten_out(i,:k0) = qrten_s(:k0) qsten_out(i,:k0) = qsten_s(:k0) qldet_out(i,:k0) = qldet_s(:k0) qidet_out(i,:k0) = qidet_s(:k0) qlsub_out(i,:k0) = qlsub_s(:k0) qisub_out(i,:k0) = qisub_s(:k0) cush_inout(i) = cush_s cufrc_out(i,:k0) = cufrc_s(:k0) qcu_out(i,:k0) = qcu_s(:k0) qlu_out(i,:k0) = qlu_s(:k0) qiu_out(i,:k0) = qiu_s(:k0) cbmf_out(i) = cbmf_s qc_out(i,:k0) = qc_s(:k0) cnt_out(i) = cnt_s cnb_out(i) = cnb_s do m = 1, ncnst trten_out(i,:k0,m) = trten_s(:k0,m) enddo ! ------------------------------------------------------------------------------ ! ! Below are diagnostic output variables for detailed analysis of cumulus scheme. ! ! The order of vertical index is reversed for this internal diagnostic output. ! ! ------------------------------------------------------------------------------ ! fer_out(i,1:k0) = fer_s(:k0) fdr_out(i,1:k0) = fdr_s(:k0) xc_out(i,1:k0) = xc_s(:k0) cinh_out(i) = cin_s cinlclh_out(i) = cinlcl_s !!! ufrc_out(i,k0:0:-1) = ufrc_s(0:k0) ! uflx_out(i,k0:0:-1) = uflx_s(0:k0) ! vflx_out(i,k0:0:-1) = vflx_s(0:k0) ufrcinvbase_out(i) = ufrcinvbase_s ufrclcl_out(i) = ufrclcl_s winvbase_out(i) = winvbase_s wlcl_out(i) = wlcl_s plcl_out(i) = plcl_s pinv_out(i) = pinv_s prel_out(i) = prel_s plfc_out(i) = plfc_s pbup_out(i) = pbup_s ppen_out(i) = ppen_s qtsrc_out(i) = qtsrc_s thlsrc_out(i) = thlsrc_s thvlsrc_out(i) = thvlsrc_s emfkbup_out(i) = emfkbup_s cbmflimit_out(i) = cbmflimit_s tkeavg_out(i) = tkeavg_s zinv_out(i) = zinv_s rcwp_out(i) = rcwp_s rlwp_out(i) = rlwp_s riwp_out(i) = riwp_s wu_out(i,k0:0:-1) = wu_s(0:k0) qtu_out(i,k0:0:-1) = qtu_s(0:k0) thlu_out(i,k0:0:-1) = thlu_s(0:k0) thvu_out(i,k0:0:-1) = thvu_s(0:k0) uu_out(i,k0:0:-1) = uu_s(0:k0) vu_out(i,k0:0:-1) = vu_s(0:k0) qtu_emf_out(i,k0:0:-1) = qtu_emf_s(0:k0) thlu_emf_out(i,k0:0:-1) = thlu_emf_s(0:k0) uu_emf_out(i,k0:0:-1) = uu_emf_s(0:k0) vu_emf_out(i,k0:0:-1) = vu_emf_s(0:k0) uemf_out(i,k0:0:-1) = uemf_s(0:k0) excessu_arr_out(i,k0:1:-1) = excessu_arr_s(:k0) excess0_arr_out(i,k0:1:-1) = excess0_arr_s(:k0) xc_arr_out(i,k0:1:-1) = xc_arr_s(:k0) aquad_arr_out(i,k0:1:-1) = aquad_arr_s(:k0) bquad_arr_out(i,k0:1:-1) = bquad_arr_s(:k0) cquad_arr_out(i,k0:1:-1) = cquad_arr_s(:k0) bogbot_arr_out(i,k0:1:-1) = bogbot_arr_s(:k0) bogtop_arr_out(i,k0:1:-1) = bogtop_arr_s(:k0) do m = 1, ncnst trflx_out(i,k0:0:-1,m) = trflx_s(0:k0,m) tru_out(i,k0:0:-1,m) = tru_s(0:k0,m) tru_emf_out(i,k0:0:-1,m) = tru_emf_s(0:k0,m) enddo id_exit = .false. go to 333 endif endif ! ------------------------------------------------------------------ ! ! Define a release level, 'prel' and release layer, 'krel'. ! ! 'prel' is the lowest level from which buoyancy sorting occurs, and ! ! 'krel' is the layer index containing 'prel' in it, similar to the ! ! previous definitions of 'kinv', 'klcl', and 'klfc'. In order to ! ! ensure that only PBL scheme works within the PBL, if LCL is below ! ! PBL top height, then 'krel = kinv', while if LCL is above PBL top ! ! height, then 'krel = klcl'. Note however that regardless of the ! ! definition of 'krel', cumulus convection induces fluxes within PBL ! ! through 'fluxbelowinv'. We can make cumulus convection start from ! ! any level, even within the PBL by appropriately defining 'krel' & ! ! 'prel' here. Then it must be accompanied by appropriate definition ! ! of source air properties, CIN, and re-setting of 'fluxbelowinv', & ! ! many other stuffs. ! ! Note that even when 'prel' is located above the PBL top height, we ! ! still have cumulus convection between PBL top height and 'prel': ! ! we simply assume that no lateral mixing occurs in this range. ! ! ------------------------------------------------------------------ ! ! print *,'klcl=',klcl if( klcl .lt. kinv ) then krel = kinv prel = pifc0(krel-1) thv0rel = thv0bot(krel) else krel = klcl prel = plcl thv0rel = thv0lcl endif ! print *,'krel=',krel ! --------------------------------------------------------------------------- ! ! Calculate cumulus base mass flux ('cbmf'), fractional area ('ufrcinv'), and ! ! and mean vertical velocity (winv) of cumulus updraft at PBL top interface. ! ! Also, calculate updraft fractional area (ufrclcl) and vertical velocity at ! ! the LCL (wlcl). When LCL is below PBLH, cinlcl = 0 and 'ufrclcl = ufrcinv', ! ! and 'wlcl = winv. ! ! Only updrafts strong enough to overcome CIN can rise over PBL top interface.! ! Thus, in order to calculate cumulus mass flux at PBL top interface, 'cbmf',! ! we need to know 'CIN' ( the strength of potential energy barrier ) and ! ! 'sigmaw' ( a standard deviation of updraft vertical velocity at the PBL top ! ! interface, a measure of turbulentce strength in the PBL ). Naturally, the ! ! ratio of these two variables, 'mu' - normalized CIN by TKE- is key variable ! ! controlling 'cbmf'. If 'mu' becomes large, only small fraction of updrafts ! ! with very strong TKE can rise over the PBL - both 'cbmf' and 'ufrc' becomes ! ! small, but 'winv' becomes large ( this can be easily understood by PDF of w ! ! at PBL top ). If 'mu' becomes small, lots of updraft can rise over the PBL ! ! top - both 'cbmf' and 'ufrc' becomes large, but 'winv' becomes small. Thus, ! ! all of the key variables associated with cumulus convection at the PBL top ! ! - 'cbmf', 'ufrc', 'winv' where 'cbmf = rho*ufrc*winv' - are a unique functi ! ! ons of 'mu', normalized CIN. Although these are uniquely determined by 'mu',! ! we usually impose two comstraints on 'cbmf' and 'ufrc': (1) because we will ! ! simply assume that subsidence warming and drying of 'kinv-1' layer in assoc ! ! iation with 'cbmf' at PBL top interface is confined only in 'kinv-1' layer, ! ! cbmf must not be larger than the mass within the 'kinv-1' layer. Otherwise, ! ! instability will occur due to the breaking of stability con. If we consider ! ! semi-Lagrangian vertical advection scheme and explicitly consider the exten ! ! t of vertical movement of each layer in association with cumulus mass flux, ! ! we don't need to impose this constraint. However, using a semi-Lagrangian ! ! scheme is a future research subject. Note that this constraint should be ap ! ! plied for all interfaces above PBL top as well as PBL top interface. As a ! ! result, this 'cbmf' constraint impose a 'lower' limit on mu - 'mumin0'. (2) ! ! in order for mass flux parameterization - rho*(w'a')= M*(a_c-a_e) - to be ! ! valid, cumulus updraft fractional area should be much smaller than 1. In ! ! current code, we impose 'rmaxfrac = 0.1 ~ 0.2' through the whole vertical ! ! layers where cumulus convection occurs. At the PBL top interface, the same ! ! constraint is made by imposing another lower 'lower' limit on mu, 'mumin1'. ! ! After that, also limit 'ufrclcl' to be smaller than 'rmaxfrac' by 'mumin2'. ! ! --------------------------------------------------------------------------- ! ! --------------------------------------------------------------------------- ! ! Calculate normalized CIN, 'mu' satisfying all the three constraints imposed ! ! on 'cbmf'('mumin0'), 'ufrc' at the PBL top - 'ufrcinv' - ( by 'mumin1' from ! ! a parameter sentence), and 'ufrc' at the LCL - 'ufrclcl' ( by 'mumin2'). ! ! Note that 'cbmf' does not change between PBL top and LCL because we assume ! ! that buoyancy sorting does not occur when cumulus updraft is unsaturated. ! ! --------------------------------------------------------------------------- if( use_CINcin ) then wcrit = sqrt( 2._r8 * cin * rbuoy ) else wcrit = sqrt( 2._r8 * cinlcl * rbuoy ) endif sigmaw = sqrt( rkfre * tkeavg + epsvarw ) mu = wcrit/sigmaw/1.4142_r8 if( mu .ge. 3._r8 ) then if (scverbose) then call write_parallel('mu >= 3') end if id_exit = .true. go to 333 endif rho0inv = pifc0(kinv-1)/(r*thv0top(kinv-1)*exnifc0(kinv-1)) cbmf = (rho0inv*sigmaw/2.5066_r8)*exp(-mu**2) ! 1. 'cbmf' constraint cbmflimit = 0.9_r8*dp0(kinv-1)/g/dt mumin0 = 0._r8 if( cbmf .gt. cbmflimit ) mumin0 = sqrt(-log(2.5066_r8*cbmflimit/rho0inv/sigmaw)) ! 2. 'ufrcinv' constraint mu = max(max(mu,mumin0),mumin1) ! 3. 'ufrclcl' constraint mulcl = sqrt(2._r8*cinlcl*rbuoy)/1.4142_r8/sigmaw mulclstar = sqrt(max(0._r8,2._r8*(exp(-mu**2)/2.5066_r8)**2*(1._r8/erfc(mu)**2-0.25_r8/rmaxfrac**2))) if( mulcl .gt. 1.e-8_r8 .and. mulcl .gt. mulclstar ) then mumin2 = compute_mumin2(mulcl,rmaxfrac,mu) if( mu .gt. mumin2 ) then call write_parallel('Critical error in mu calculation in UW_ShCu') ! call endrun endif mu = max(mu,mumin2) if( mu .eq. mumin2 ) limit_ufrc(i) = 1. endif if( mu .eq. mumin0 ) limit_cbmf(i) = 1. if( mu .eq. mumin1 ) limit_ufrc(i) = 1. ! ------------------------------------------------------------------- ! ! Calculate final ['cbmf','ufrcinv','winv'] at the PBL top interface. ! ! Note that final 'cbmf' here is obtained in such that 'ufrcinv' and ! ! 'ufrclcl' are smaller than ufrcmax with no instability. ! ! ------------------------------------------------------------------- ! cbmf = (rho0inv*sigmaw/2.5066_r8)*exp(-mu**2) winv = sigmaw*(2._r8/2.5066_r8)*exp(-mu**2)/erfc(mu) ufrcinv = cbmf/winv/rho0inv ! ------------------------------------------------------------------- ! ! Calculate ['ufrclcl','wlcl'] at the LCL. When LCL is below PBL top, ! ! it automatically becomes 'ufrclcl = ufrcinv' & 'wlcl = winv', since ! ! it was already set to 'cinlcl=0' if LCL is below PBL top interface. ! ! Note 'cbmf' at the PBL top is the same as 'cbmf' at the LCL. Note ! ! also that final 'cbmf' here is obtained in such that 'ufrcinv' and ! ! 'ufrclcl' are smaller than ufrcmax and there is no instability. ! ! By construction, it must be 'wlcl > 0' but for assurance, I checked ! ! this again in the below block. If 'ufrclcl < 0.1%', just exit. ! ! ------------------------------------------------------------------- ! wtw = winv * winv - 2._r8 * cinlcl * rbuoy if( wtw .le. 0._r8 ) then if (scverbose) then call write_parallel('wlcl < 0 at the LCL') end if exit_wtw(i) = 1. id_exit = .true. go to 333 endif wlcl = sqrt(wtw) ufrclcl = cbmf/wlcl/rho0inv wrel = wlcl if( ufrclcl .le. 0.0001_r8 ) then if (scverbose) then call write_parallel( 'ufrclcl <= 0.0001' ) end if exit_ufrc(i) = 1. id_exit = .true. go to 333 endif ufrc(krel-1) = ufrclcl ! ----------------------------------------------------------------------- ! ! Below is just diagnostic output for detailed analysis of cumulus scheme ! ! ----------------------------------------------------------------------- ! ufrcinvbase = ufrcinv winvbase = winv umf(kinv-1:krel-1) = cbmf wu(kinv-1:krel-1) = winv ! -------------------------------------------------------------------------- ! ! Define updraft properties at the level where buoyancy sorting starts to be ! ! happening, i.e., by definition, at 'prel' level within the release layer. ! ! Because no lateral entrainment occurs upto 'prel', conservative scalars of ! ! cumulus updraft at release level is same as those of source air. However, ! ! horizontal momentums of source air are modified by horizontal PGF forcings ! ! from PBL top interface to 'prel'. For this case, we should add additional ! ! horizontal momentum from PBL top interface to 'prel' as will be done below ! ! to 'usrc' and 'vsrc'. Note that below cumulus updraft properties - umf, wu,! ! thlu, qtu, thvu, uu, vu - are defined all interfaces not at the layer mid- ! ! point. From the index notation of cumulus scheme, wu(k) is the cumulus up- ! ! draft vertical velocity at the top interface of k layer. ! ! Diabatic horizontal momentum forcing should be treated as a kind of 'body' ! ! forcing without actual mass exchange between convective updraft and ! ! environment, but still taking horizontal momentum from the environment to ! ! the convective updrafts. Thus, diabatic convective momentum transport ! ! vertically redistributes environmental horizontal momentum. ! ! -------------------------------------------------------------------------- ! emf(krel-1) = 0._r8 umf(krel-1) = cbmf wu(krel-1) = wrel thlu(krel-1) = thlsrc qtu(krel-1) = qtsrc call conden(prel,thlsrc,qtsrc,thj,qvj,qlj,qij,qse,id_check) if( id_check .eq. 1 ) then exit_conden(i) = 1. id_exit = .true. if (scverbose) then call write_parallel('------- UW ShCu: exit, conden') end if go to 333 endif thvu(krel-1) = thj * ( 1._r8 + zvir*qvj - qlj - qij ) uplus = 0._r8 vplus = 0._r8 if( krel .eq. kinv ) then uplus = PGFc * ssu0(kinv) * ( prel - pifc0(kinv-1) ) vplus = PGFc * ssv0(kinv) * ( prel - pifc0(kinv-1) ) else do k = kinv, max(krel-1,kinv) uplus = uplus + PGFc * ssu0(k) * ( pifc0(k) - pifc0(k-1) ) vplus = vplus + PGFc * ssv0(k) * ( pifc0(k) - pifc0(k-1) ) end do uplus = uplus + PGFc * ssu0(krel) * ( prel - pifc0(krel-1) ) vplus = vplus + PGFc * ssv0(krel) * ( prel - pifc0(krel-1) ) end if uu(krel-1) = usrc + uplus vu(krel-1) = vsrc + vplus do m = 1, ncnst tru(krel-1,m) = trsrc(m) enddo ! -------------------------------------------------------------------------- ! ! Define environmental properties at the level where buoyancy sorting occurs ! ! ('pe', normally, layer midpoint except in the 'krel' layer). In the 'krel' ! ! layer where buoyancy sorting starts to occur, however, 'pe' is defined ! ! differently because LCL is regarded as lower interface for mixing purpose. ! ! -------------------------------------------------------------------------- ! pe = 0.5_r8 * ( prel + pifc0(krel) ) qsat_pe = 0.5_r8 * ( prel + pifc0(krel) ) dpe = prel - pifc0(krel) exne = exnerfn(pe) thvebot = thv0rel thle = thl0(krel) + ssthl0(krel) * ( pe - pmid0(krel) ) qte = qt0(krel) + ssqt0(krel) * ( pe - pmid0(krel) ) ue = u0(krel) + ssu0(krel) * ( pe - pmid0(krel) ) ve = v0(krel) + ssv0(krel) * ( pe - pmid0(krel) ) do m = 1, ncnst tre(m) = tr0(krel,m) + sstr0(krel,m) * ( pe - pmid0(krel) ) enddo !-------------------------! ! Buoyancy-Sorting Mixing ! !-------------------------!------------------------------------------------ ! ! ! ! In order to complete buoyancy-sorting mixing at layer mid-point, and so ! ! calculate 'updraft mass flux, updraft w velocity, conservative scalars' ! ! at the upper interface of each layer, we need following 3 information. ! ! ! ! 1. Pressure where mixing occurs ('pe'), and temperature at 'pe' which is ! ! necessary to calculate various thermodynamic coefficients at pe. This ! ! temperature is obtained by undiluted cumulus properties lifted to pe. ! ! 2. Undiluted updraft properties at pe - conservative scalar and vertical ! ! velocity -which are assumed to be the same as the properties at lower ! ! interface only for calculation of fractional lateral entrainment and ! ! detrainment rate ( fer(k) and fdr(k) [Pa-1] ), respectively. Final ! ! values of cumulus conservative scalars and w at the top interface are ! ! calculated afterward after obtaining fer(k) & fdr(k). ! ! 3. Environmental properties at pe. ! ! ------------------------------------------------------------------------- ! ! ------------------------------------------------------------------------ ! ! Define cumulus scale height. ! ! Cumulus scale height is defined as the maximum height cumulus can reach. ! ! In case of premitive code, cumulus scale height ('cush') at the current ! ! time step was assumed to be the same as 'cush' of previous time step. ! ! However, I directly calculated cush at each time step using an iterative ! ! method. Note that within the cumulus scheme, 'cush' information is used ! ! only at two places during buoyancy-sorting process: ! ! (1) Even negatively buoyancy mixtures with strong vertical velocity ! ! enough to rise up to 'rle*scaleh' (rle = 0.1) from pe are entrained ! ! into cumulus updraft, ! ! (2) The amount of mass that is involved in buoyancy-sorting mixing ! ! process at pe is rei(k) = rkm/scaleh/rho*g [Pa-1] ! ! In terms of (1), I think critical stopping distance might be replaced by ! ! layer thickness. In future, we will use rei(k) = (0.5*rkm/z0(k)/rho/g). ! ! In the premitive code, 'scaleh' was largely responsible for the jumping ! ! variation of precipitation amount. ! ! ------------------------------------------------------------------------ ! scaleh = tscaleh if( tscaleh .lt. 0.0_r8 ) scaleh = 1000._r8 ! Save time : Set iter_scaleh = 1. This will automatically use 'cush' from the previous time step ! at the first implicit iteration. At the second implicit iteration, it will use ! the updated 'cush' by the first implicit cin. So, this updating has an effect of ! doing one iteration for cush calculation, which is good. ! So, only this setting of 'iter_scaleh = 1' is sufficient-enough to save computation time. ! OK do iter_scaleh = 1, 3 ! iter_scaleh = 1 ! ---------------------------------------------------------------- ! ! Initialization of 'kbup' and 'kpen' ! ! ---------------------------------------------------------------- ! ! 'kbup' is the top-most layer in which cloud buoyancy is positive ! ! both at the top and bottom interface of the layer. 'kpen' is the ! ! layer upto which cumulus panetrates ,i.e., cumulus w at the base ! ! interface is positive, but becomes negative at the top interface.! ! Here, we initialize 'kbup' and 'kpen'. These initializations are ! ! not trivial but important, expecially in calculating turbulent ! ! fluxes without confliction among several physics as explained in ! ! detail in the part of turbulent fluxes calculation later. Note ! ! that regardless of whether 'kbup' and 'kpen' are updated or not ! ! during updraft motion, penetrative entrainments are dumped down ! ! across the top interface of 'kbup' later. More specifically,! ! penetrative entrainment heat and moisture fluxes are calculated ! ! from the top interface of 'kbup' layer to the base interface of ! ! 'kpen' layer. Because of this, initialization of 'kbup' & 'kpen' ! ! influence the convection system when there are not updated. The ! ! below initialization of 'kbup = krel' assures that penetrative ! ! entrainment fluxes always occur at interfaces above the PBL top ! ! interfaces (i.e., only at interfaces k >=kinv ), which seems to ! ! be attractable considering that the most correct fluxes at the ! ! PBL top interface can be ontained from the 'fluxbelowinv' using ! ! reconstructed PBL height. ! ! The 'kbup = krel'(after going through the whole buoyancy sorting ! ! proces during updraft motion) implies that cumulus updraft from ! ! the PBL top interface can not reach to the LFC,so that 'kbup' is ! ! not updated during upward. This means that cumulus updraft did ! ! not fully overcome the buoyancy barrier above just the PBL top. ! ! If 'kpen' is not updated either ( i.e., cumulus cannot rise over ! ! the top interface of release layer),penetrative entrainment will ! ! not happen at any interfaces. If cumulus updraft can rise above ! ! the release layer but cannot fully overcome the buoyancy barrier ! ! just above PBL top interface, penetratve entrainment occurs at ! ! several above interfaces, including the top interface of release ! ! layer. In the latter case, warming and drying tendencies will be ! ! be initiated in 'krel' layer. Note current choice of 'kbup=krel' ! ! is completely compatible with other flux physics without double ! ! or miss counting turbulent fluxes at any interface. However, the ! ! alternative choice of 'kbup=krel-1' also has itw own advantage - ! ! when cumulus updraft cannot overcome buoyancy barrier just above ! ! PBL top, entrainment warming and drying are concentrated in the ! ! 'kinv-1' layer instead of 'kinv' layer for this case. This might ! ! seems to be more dynamically reasonable, but I will choose the ! ! 'kbup = krel' choice since it is more compatible with the other ! ! parts of the code, expecially, when we chose ' use_emf=.false. ' ! ! as explained in detail in turbulent flux calculation part. ! ! ---------------------------------------------------------------- ! kbup = krel kpen = krel ! ------------------------------------------------------------ ! ! Since 'wtw' is continuously updated during vertical motion, ! ! I need below initialization command within this 'iter_scaleh'! ! do loop. Similarily, I need initializations of environmental ! ! properties at 'krel' layer as below. ! ! ------------------------------------------------------------ ! wtw = wlcl * wlcl pe = 0.5_r8 * ( prel + pifc0(krel) ) dpe = prel - pifc0(krel) exne = exnerfn(pe) thvebot = thv0rel thle = thl0(krel) + ssthl0(krel) * ( pe - pmid0(krel) ) qte = qt0(krel) + ssqt0(krel) * ( pe - pmid0(krel) ) ue = u0(krel) + ssu0(krel) * ( pe - pmid0(krel) ) ve = v0(krel) + ssv0(krel) * ( pe - pmid0(krel) ) do m = 1, ncnst tre(m) = tr0(krel,m) + sstr0(krel,m) * ( pe - pmid0(krel) ) enddo ! ----------------------------------------------------------------------- ! ! Cumulus rises upward from 'prel' ( or base interface of 'krel' layer ) ! ! until updraft vertical velocity becomes zero. ! ! Buoyancy sorting is performed via two stages. (1) Using cumulus updraft ! ! properties at the base interface of each layer,perform buoyancy sorting ! ! at the layer mid-point, 'pe', and update cumulus properties at the top ! ! interface, and then (2) by averaging updated cumulus properties at the ! ! top interface and cumulus properties at the base interface, calculate ! ! cumulus updraft properties at pe that will be used in buoyancy sorting ! ! mixing - thlue, qtue and, wue. Using this averaged properties, perform ! ! buoyancy sorting again at pe, and re-calculate fer(k) and fdr(k). Using ! ! this recalculated fer(k) and fdr(k), finally calculate cumulus updraft ! ! properties at the top interface - thlu, qtu, thvu, uu, vu. In the below,! ! 'iter_xc = 1' performs the first stage, while 'iter_xc= 2' performs the ! ! second stage. We can increase the number of iterations, 'nter_xc'.as we ! ! want, but a sample test indicated that about 3 - 5 iterations produced ! ! satisfactory converent solution. Finally, identify 'kbup' and 'kpen'. ! ! ----------------------------------------------------------------------- ! do k = krel, k0 - 1 ! Here, 'k' is a layer index. km1 = k - 1 thlue = thlu(km1) qtue = qtu(km1) wue = wu(km1) wtwb = wtw do iter_xc = 1, niter_xc wtw = wu(km1) * wu(km1) ! ---------------------------------------------------------------- ! ! Calculate environmental and cumulus saturation 'excess' at 'pe'. ! ! Note that in order to calculate saturation excess, we should use ! ! liquid water temperature instead of temperature as the argument ! ! of "qsat". But note normal argument of "qsat" is temperature. ! ! ---------------------------------------------------------------- ! call conden(pe,thle,qte,thj,qvj,qlj,qij,qse,id_check) if( id_check .eq. 1 ) then exit_conden(i) = 1. id_exit = .true. if (scverbose) then call write_parallel('------- UW ShCu: exit, conden') end if go to 333 end if thv0j = thj * ( 1._r8 + zvir*qvj - qlj - qij ) rhomid0j = pe / ( r * thv0j * exne ) qsat_arg = thle*exne qs = GEOS_QSAT(qsat_arg,qsat_pe/100.) excess0 = qte - qs call conden(pe,thlue,qtue,thj,qvj,qlj,qij,qse,id_check) if( id_check .eq. 1 ) then exit_conden(i) = 1. id_exit = .true. if (scverbose) then call write_parallel('------- UW ShCu: exit, conden') end if go to 333 end if ! ----------------------------------------------------------------- ! ! Detrain excessive condensate larger than 'criqc' from the cumulus ! ! updraft before performing buoyancy sorting. All I should to do is ! ! to update 'thlue' & 'que' here. Below modification is completely ! ! compatible with the other part of the code since 'thule' & 'qtue' ! ! are used only for buoyancy sorting. I found that as long as I use ! ! 'niter_xc >= 2', detraining excessive condensate before buoyancy ! ! sorting has negligible influence on the buoyancy sorting results. ! ! ----------------------------------------------------------------- ! if( (qlj + qij) .gt. criqc ) then exql = ( ( qlj + qij ) - criqc ) * qlj / ( qlj + qij ) exqi = ( ( qlj + qij ) - criqc ) * qij / ( qlj + qij ) qtue = qtue - exql - exqi thlue = thlue + (xlv/cp/exne)*exql + (xls/cp/exne)*exqi endif call conden(pe,thlue,qtue,thj,qvj,qlj,qij,qse,id_check) if( id_check .eq. 1 ) then exit_conden(i) = 1. id_exit = .true. if (scverbose) then call write_parallel('------- UW ShCu: exit, conden') end if go to 333 end if thvj = thj * ( 1. + zvir * qvj - qlj - qij ) tj = thj * exne ! This 'tj' is used for computing thermo. coeffs. below qsat_arg = thlue*exne qs = GEOS_QSAT(qsat_arg,qsat_pe/100.) excessu = qtue - qs ! ------------------------------------------------------------------- ! ! Calculate critical mixing fraction, 'xc'. Mixture with mixing ratio ! ! smaller than 'xc' will be entrained into cumulus updraft. Both the ! ! saturated updrafts with 'positive buoyancy' or 'negative buoyancy + ! ! strong vertical velocity enough to rise certain threshold distance' ! ! are kept into the updraft in the below program. If the core updraft ! ! is unsaturated, we can set 'xc = 0' and let the cumulus convection ! ! still works or we may exit. ! ! Current below code does not entrain unsaturated mixture. However it ! ! should be modified such that it also entrain unsaturated mixture. ! ! ------------------------------------------------------------------- ! ! ----------------------------------------------------------------- ! ! cridis : Critical stopping distance for buoyancy sorting purpose. ! ! scaleh is only used here. ! ! ----------------------------------------------------------------- ! cridis = rle*scaleh ! Original code ! cridis = 1._r8*(zifc0(k) - zifc0(k-1)) ! New code ! ---------------- ! ! Buoyancy Sorting ! ! ---------------- ! ! ----------------------------------------------------------------- ! ! Case 1 : When both cumulus and env. are unsaturated or saturated. ! ! ----------------------------------------------------------------- ! if( ( excessu .le. 0._r8 .and. excess0 .le. 0._r8 ) .or. ( excessu .ge. 0._r8 .and. excess0 .ge. 0._r8 ) ) then xc = min(1._r8,max(0._r8,1._r8-2._r8*rbuoy*g*cridis/wue**2.*(1._r8-thvj/thv0j))) ! Below 3 lines are diagnostic output not influencing ! numerical calculations. aquad = 0._r8 bquad = 0._r8 cquad = 0._r8 else ! -------------------------------------------------- ! ! Case 2 : When either cumulus or env. is saturated. ! ! -------------------------------------------------- ! xsat = excessu / ( excessu - excess0 ); thlxsat = thlue + xsat * ( thle - thlue ); qtxsat = qtue + xsat * ( qte - qtue ); call conden(pe,thlxsat,qtxsat,thj,qvj,qlj,qij,qse,id_check) if( id_check .eq. 1 ) then exit_conden(i) = 1. id_exit = .true. if (scverbose) then call write_parallel('------- UW ShCu: exit, conden') end if go to 333 end if thvxsat = thj * ( 1._r8 + zvir * qvj - qlj - qij ) ! -------------------------------------------------- ! ! kk=1 : Cumulus Segment, kk=2 : Environment Segment ! ! -------------------------------------------------- ! do kk = 1, 2 if( kk .eq. 1 ) then thv_x0 = thvj; thv_x1 = ( 1._r8 - 1._r8/xsat ) * thvj + ( 1._r8/xsat ) * thvxsat; else thv_x1 = thv0j; thv_x0 = ( xsat / ( xsat - 1._r8 ) ) * thv0j + ( 1._r8/( 1._r8 - xsat ) ) * thvxsat; endif aquad = wue**2; bquad = 2._r8*rbuoy*g*cridis*(thv_x1 - thv_x0)/thv0j - 2._r8*wue**2; cquad = 2._r8*rbuoy*g*cridis*(thv_x0 - thv0j)/thv0j + wue**2; if( kk .eq. 1 ) then if( ( bquad**2-4._r8*aquad*cquad ) .ge. 0._r8 ) then call roots(aquad,bquad,cquad,xs1,xs2,status) x_cu = min(1._r8,max(0._r8,min(xsat,min(xs1,xs2)))) else x_cu = xsat; endif else if( ( bquad**2-4._r8*aquad*cquad) .ge. 0._r8 ) then call roots(aquad,bquad,cquad,xs1,xs2,status) x_en = min(1._r8,max(0._r8,max(xsat,min(xs1,xs2)))) else x_en = 1._r8 endif endif enddo if( x_cu .eq. xsat ) then xc = max(x_cu, x_en) else xc = x_cu endif endif ! ------------------------------------------------------------------------ ! ! Compute fractional lateral entrainment & detrainment rate in each layers.! ! The unit of rei(k), fer(k), and fdr(k) is [Pa-1]. Alternative choice of ! ! 'rei(k)' is also shown below, where coefficient 0.5 was from approximate ! ! tuning against the BOMEX case. ! ! In order to prevent the onset of instability in association with cumulus ! ! induced subsidence advection, cumulus mass flux at the top interface in ! ! any layer should be smaller than ( 90% of ) total mass within that layer.! ! I imposed limits on 'rei(k)' as below, in such that stability condition ! ! is always satisfied. ! ! Below limiter of 'rei(k)' becomes negative for some cases, causing error.! ! So, for the time being, I came back to the original limiter. ! ! ------------------------------------------------------------------------ ! ee2 = xc**2 ud2 = 1._r8 - 2._r8*xc + xc**2 ! rei(k) = ( rkm / scaleh / g / rho0j ) ! Default. rei(k) = ( 0.5_r8 * rkm / zmid0(k) / g /rhomid0j ) ! Alternative. if( xc .gt. 0.5_r8 ) rei(k) = min(rei(k),0.9_r8*log(dp0(k)/g/dt/umf(km1) + 1._r8)/dpe/(2._r8*xc-1.)) fer(k) = rei(k) * ee2 fdr(k) = rei(k) * ud2 xco(k) = xc ! ------------------------------------------------------------------------------ ! ! Iteration Start due to 'maxufrc' constraint [ ****************************** ] ! ! ------------------------------------------------------------------------------ ! ! -------------------------------------------------------------------------- ! ! Calculate cumulus updraft mass flux and penetrative entrainment mass flux. ! ! Note that non-zero penetrative entrainment mass flux will be asigned only ! ! to interfaces from the top interface of 'kbup' layer to the base interface ! ! of 'kpen' layer as will be shown later. ! ! -------------------------------------------------------------------------- ! umf(k) = umf(km1) * exp( dpe * ( fer(k) - fdr(k) ) ) emf(k) = 0._r8 ! --------------------------------------------------------- ! ! Compute cumulus updraft properties at the top interface. ! ! Also use Tayler expansion in order to treat limiting case ! ! --------------------------------------------------------- ! if( fer(k)*dpe .lt. 1.e-4_r8 ) then thlu(k) = thlu(km1) + ( thle + ssthl0(k) * dpe / 2._r8 - thlu(km1) ) * fer(k) * dpe qtu(k) = qtu(km1) + ( qte + ssqt0(k) * dpe / 2._r8 - qtu(km1) ) * fer(k) * dpe uu(k) = uu(km1) + ( ue + ssu0(k) * dpe / 2._r8 - uu(km1) ) * fer(k) * dpe - PGFc * ssu0(k) * dpe vu(k) = vu(km1) + ( ve + ssv0(k) * dpe / 2._r8 - vu(km1) ) * fer(k) * dpe - PGFc * ssv0(k) * dpe do m = 1, ncnst tru(k,m) = tru(km1,m) + ( tre(m) + sstr0(k,m) * dpe / 2._r8 - tru(km1,m) ) * fer(k) * dpe enddo else thlu(k) = ( thle + ssthl0(k) / fer(k) - ssthl0(k) * dpe / 2._r8 ) - & ( thle + ssthl0(k) * dpe / 2._r8 - thlu(km1) + ssthl0(k) / fer(k) ) * exp(-fer(k) * dpe) qtu(k) = ( qte + ssqt0(k) / fer(k) - ssqt0(k) * dpe / 2._r8 ) - & ( qte + ssqt0(k) * dpe / 2._r8 - qtu(km1) + ssqt0(k) / fer(k) ) * exp(-fer(k) * dpe) uu(k) = ( ue + ( 1._r8 - PGFc ) * ssu0(k) / fer(k) - ssu0(k) * dpe / 2._r8 ) - & ( ue + ssu0(k) * dpe / 2._r8 - uu(km1) + ( 1._r8 - PGFc ) * ssu0(k) / fer(k) ) * exp(-fer(k) * dpe) vu(k) = ( ve + ( 1._r8 - PGFc ) * ssv0(k) / fer(k) - ssv0(k) * dpe / 2._r8 ) - & ( ve + ssv0(k) * dpe / 2._r8 - vu(km1) + ( 1._r8 - PGFc ) * ssv0(k) / fer(k) ) * exp(-fer(k) * dpe) do m = 1, ncnst tru(k,m) = ( tre(m) + sstr0(k,m) / fer(k) - sstr0(k,m) * dpe / 2._r8 ) - & ( tre(m) + sstr0(k,m) * dpe / 2._r8 - tru(km1,m) + sstr0(k,m) / fer(k) ) * exp(-fer(k) * dpe) enddo end if !------------------------------------------------------------------- ! ! Expel some of cloud water and ice from cumulus updraft at the top ! ! interface. Note that this is not 'detrainment' term but a 'sink' ! ! term of cumulus updraft qt ( or one part of 'source' term of mean ! ! environmental qt ). At this stage, as the most simplest choice, if ! ! condensate amount within cumulus updraft is larger than a critical ! ! value, 'criqc', expels the surplus condensate from cumulus updraft ! ! to the environment. A certain fraction ( e.g., 'frc_sus' ) of this ! ! expelled condesnate will be in a form that can be suspended in the ! ! layer k where it was formed, while the other fraction, '1-frc_sus' ! ! will be in a form of precipitatble (e.g.,can potentially fall down ! ! across the base interface of layer k ). In turn we should describe ! ! subsequent falling of precipitable condensate ('1-frc_sus') across ! ! the base interface of the layer k, & evaporation of precipitating ! ! water in the below layer k-1 and associated evaporative cooling of ! ! the later, k-1, and falling of 'non-evaporated precipitating water ! ! ( which was initially formed in layer k ) and a newly-formed preci ! ! pitable water in the layer, k-1', across the base interface of the ! ! lower layer k-1. Cloud microphysics should correctly describe all ! ! of these process. In a near future, I should significantly modify ! ! this cloud microphysics, including precipitation-induced downdraft ! ! also. ! ! ------------------------------------------------------------------ ! call conden(pifc0(k),thlu(k),qtu(k),thj,qvj,qlj,qij,qse,id_check) if( id_check .eq. 1 ) then exit_conden(i) = 1. id_exit = .true. if (scverbose) then call write_parallel('------- UW ShCu: exit, conden') end if go to 333 end if if( (qlj + qij) .gt. criqc ) then exql = ( ( qlj + qij ) - criqc ) * qlj / ( qlj + qij ) exqi = ( ( qlj + qij ) - criqc ) * qij / ( qlj + qij ) ! ---------------------------------------------------------------- ! ! It is very important to re-update 'qtu' and 'thlu' at the upper ! ! interface after expelling condensate from cumulus updraft at the ! ! top interface of the layer. As mentioned above, this is a 'sink' ! ! of cumulus qt (or equivalently, a 'source' of environmentasl qt),! ! not a regular convective'detrainment'. ! ! ---------------------------------------------------------------- ! qtu(k) = qtu(k) - exql - exqi thlu(k) = thlu(k) + (xlv/exnifc0(k)/cp)*exql + (xls/exnifc0(k)/cp)*exqi ! ---------------------------------------------------------------- ! ! Expelled cloud condensate into the environment from the updraft. ! ! After all the calculation later, 'dwten' and 'diten' will have a ! ! unit of [ kg/kg/s ], because it is a tendency of qt. Restoration ! ! of 'dwten' and 'diten' to this correct unit through multiplying ! ! 'umf(k)*g/dp0(k)' will be performed later after finally updating ! ! 'umf' using a 'rmaxfrac' constraint near the end of this updraft ! ! buoyancy sorting loop. ! ! ---------------------------------------------------------------- ! dwten(k) = exql diten(k) = exqi else dwten(k) = 0._r8 diten(k) = 0._r8 endif ! ----------------------------------------------------------------- ! ! Update 'thvu(k)' after detraining condensate from cumulus updraft.! ! ----------------------------------------------------------------- ! call conden(pifc0(k),thlu(k),qtu(k),thj,qvj,qlj,qij,qse,id_check) if( id_check .eq. 1 ) then exit_conden(i) = 1. id_exit = .true. if (scverbose) then call write_parallel('------- UW ShCu: exit, conden') end if go to 333 end if thvu(k) = thj * ( 1._r8 + zvir * qvj - qlj - qij ) ! ----------------------------------------------------------- ! ! Calculate updraft vertical velocity at the upper interface. ! ! In order to calculate 'wtw' at the upper interface, we use ! ! 'wtw' at the lower interface. Note 'wtw' is continuously ! ! updated as cumulus updraft rises. ! ! ----------------------------------------------------------- ! bogbot = rbuoy * ( thvu(km1) / thvebot - 1._r8 ) ! Cloud buoyancy at base interface bogtop = rbuoy * ( thvu(k) / thv0top(k) - 1._r8 ) ! Cloud buoyancy at top interface delbog = bogtop - bogbot drage = fer(k) * ( 1._r8 + rdrag ) expfac = exp(-2._r8*drage*dpe) wtwb = wtw if( drage*dpe .gt. 1.e-3_r8 ) then wtw = wtw*expfac + (delbog + (1._r8-expfac)*(bogbot + delbog/(-2._r8*drage*dpe)))/(rhomid0j*drage) else wtw = wtw + dpe * ( bogbot + bogtop ) / rhomid0j endif ! Force the plume rise at least to klfc of the undiluted plume. ! Because even the below is not complete, I decided not to include this. ! if( k .le. klfc ) then ! wtw = max( 1.e-2, wtw ) ! endif ! -------------------------------------------------------------- ! ! Repeat 'iter_xc' iteration loop until 'iter_xc = niter_xc'. ! ! Also treat the case even when wtw < 0 at the 'kpen' interface. ! ! -------------------------------------------------------------- ! if( wtw .gt. 0. ) then thlue = 0.5_r8 * ( thlu(km1) + thlu(k) ) qtue = 0.5_r8 * ( qtu(km1) + qtu(k) ) wue = 0.5_r8 * sqrt( max( wtwb + wtw, 0._r8 ) ) else go to 111 endif end do ! iter_xc loop 111 continue ! --------------------------------------------------------------------------- ! ! Add the contribution of self-detrainment to vertical variations of cumulus ! ! updraft mass flux. The reason why we are trying to include self-detrainment ! ! is as follows. In current scheme, vertical variation of updraft mass flux ! ! is not fully consistent with the vertical variation of updraft vertical w. ! ! For example, within a given layer, let's assume that cumulus w is positive ! ! at the base interface, while negative at the top interface. This means that ! ! cumulus updraft cannot reach to the top interface of the layer. However, ! ! cumulus updraft mass flux at the top interface is not zero according to the ! ! vertical tendency equation of cumulus mass flux. Ideally, cumulus updraft ! ! mass flux at the top interface should be zero for this case. In order to ! ! assures that cumulus updraft mass flux goes to zero when cumulus updraft ! ! vertical velocity goes to zero, we are imposing self-detrainment term as ! ! below by considering layer-mean cloud buoyancy and cumulus updraft vertical ! ! velocity square at the top interface. Use of auto-detrainment term will be ! ! determined by setting 'use_self_detrain=.true.' in the parameter sentence. ! ! --------------------------------------------------------------------------- ! if( use_self_detrain ) then autodet = min( 0.5_r8*g*(bogbot+bogtop)/(max(wtw,0._r8)+1.e-4_r8), 0._r8 ) umf(k) = umf(k) * exp( 0.637_r8*(dpe/rhomid0j/g) * autodet ) end if if( umf(k) .eq. 0._r8 ) wtw = -1._r8 ! -------------------------------------- ! ! Below block is just a dignostic output ! ! -------------------------------------- ! excessu_arr(k) = excessu excess0_arr(k) = excess0 xc_arr(k) = xc aquad_arr(k) = aquad bquad_arr(k) = bquad cquad_arr(K) = cquad bogbot_arr(k) = bogbot bogtop_arr(k) = bogtop ! ------------------------------------------------------------------- ! ! 'kbup' is the upper most layer in which cloud buoyancy is positive ! ! both at the base and top interface. 'kpen' is the upper most layer ! ! up to cumulus can reach. Usually, 'kpen' is located higher than the ! ! 'kbup'. Note we initialized these by 'kbup = krel' & 'kpen = krel'. ! ! As explained before, it is possible that only 'kpen' is updated, ! ! while 'kbup' keeps its initialization value. For this case, current ! ! scheme will simply turns-off penetrative entrainment fluxes and use ! ! normal buoyancy-sorting fluxes for 'kbup <= k <= kpen-1' interfaces,! ! in order to describe shallow continental cumulus convection. ! ! ------------------------------------------------------------------- ! ! if( bogbot .gt. 0. .and. bogtop .gt. 0. ) then ! if( bogtop .gt. 0. ) then if( bogtop .gt. 0._r8 .and. wtw .gt. 0._r8 ) then kbup = k end if if( wtw .le. 0._r8 ) then kpen = k go to 45 end if wu(k) = sqrt(wtw) if( wu(k) .gt. 100._r8 ) then exit_wu(i) = 1. id_exit = .true. if (scverbose) then call write_parallel('------- UW ShCu: exited, wu>100') end if go to 333 endif ! ---------------------------------------------------------------------------- ! ! Iteration end due to 'rmaxfrac' constraint [ ***************************** ] ! ! ---------------------------------------------------------------------------- ! ! ---------------------------------------------------------------------- ! ! Calculate updraft fractional area at the upper interface and set upper ! ! limit to 'ufrc' by 'rmaxfrac'. In order to keep the consistency among ! ! ['ufrc','umf','wu (or wtw)'], if ufrc is limited by 'rmaxfrac', either ! ! 'umf' or 'wu' should be changed. Although both 'umf' and 'wu (wtw)' at ! ! the current upper interface are used for updating 'umf' & 'wu' at the ! ! next upper interface, 'umf' is a passive variable not influencing the ! ! buoyancy sorting process in contrast to 'wtw'. This is a reason why we ! ! adjusted 'umf' instead of 'wtw'. In turn we updated 'fdr' here instead ! ! of 'fer', which guarantees that all previously updated thermodynamic ! ! variables at the upper interface before applying 'rmaxfrac' constraint ! ! are already internally consistent, even though 'ufrc' is limited by ! ! 'rmaxfrac'. Thus, we don't need to go through interation loop again.If ! ! If we update 'fer' however, we should go through above iteration loop. ! ! ---------------------------------------------------------------------- ! rhoifc0j = pifc0(k) / ( r * 0.5_r8 * ( thv0bot(k+1) + thv0top(k) )*exnifc0(k) ) ufrc(k) = umf(k) / ( rhoifc0j * wu(k) ) if( ufrc(k) .gt. rmaxfrac ) then limit_ufrc(i) = 1. ufrc(k) = rmaxfrac umf(k) = rmaxfrac * rhoifc0j * wu(k) fdr(k) = fer(k) - log( umf(k) / umf(km1) ) / dpe endif ! ------------------------------------------------------------ ! ! Update environmental properties for at the mid-point of next ! ! upper layer for use in buoyancy sorting. ! ! ------------------------------------------------------------ ! pe = pmid0(k+1) dpe = dp0(k+1) exne = exnmid0(k+1) thvebot = thv0bot(k+1) thle = thl0(k+1) qte = qt0(k+1) ue = u0(k+1) ve = v0(k+1) do m = 1, ncnst tre(m) = tr0(k+1,m) enddo end do ! k loop ! ------------------------------------------------------------------------------- ! ! Up to this point, we finished all of buoyancy sorting processes from the 'krel' ! ! layer to 'kpen' layer: at the top interface of individual layers, we calculated ! ! updraft and penetrative mass fluxes [ umf(k) & emf(k) = 0 ], updraft fractional ! ! area [ ufrc(k) ], updraft vertical velocity [ wu(k) ], updraft thermodynamic ! ! variables [thlu(k),qtu(k),uu(k),vu(k),thvu(k)]. In the layer,we also calculated ! ! fractional entrainment-detrainment rate [ fer(k), fdr(k) ], and detrainment ten ! ! dency of water and ice from cumulus updraft [ dwten(k), diten(k) ]. In addition,! ! we updated and identified 'krel' and 'kpen' layer index, if any. In the 'kpen' ! ! layer, we calculated everything mentioned above except the 'wu(k)' and 'ufrc(k)'! ! since a real value of updraft vertical velocity is not defined at the kpen top ! ! interface (note 'ufrc' at the top interface of layer is calculated from 'umf(k)'! ! and 'wu(k)'). As mentioned before, special treatment is required when 'kbup' is ! ! not updated and so 'kbup = krel'. ! ! ------------------------------------------------------------------------------- ! ! ------------------------------------------------------------------------------ ! ! During the 'iter_scaleh' iteration loop, non-physical ( with non-zero values ) ! ! values can remain in the variable arrays above (also 'including' in case of wu ! ! and ufrc at the top interface) the 'kpen' layer. This can happen when the kpen ! ! layer index identified from the 'iter_scaleh = 1' iteration loop is located at ! ! above the kpen layer index identified from 'iter_scaleh = 3' iteration loop. ! ! Thus, in the following calculations, we should only use the values in each ! ! variables only up to finally identified 'kpen' layer & 'kpen' interface except ! ! 'wu' and 'ufrc' at the top interface of 'kpen' layer. Note that in order to ! ! prevent any problems due to these non-physical values, I re-initialized the ! ! values of [ umf(kpen:k0), emf(kpen:k0), dwten(kpen+1:k0), diten(kpen+1:k0),! ! fer(kpen:k0), fdr(kpen+1:k0), ufrc(kpen:k0) ] to be zero after 'iter_scaleh'! ! do loop. ! ! ------------------------------------------------------------------------------ ! 45 continue ! ------------------------------------------------------------------------------ ! ! Calculate 'ppen( < 0 )', updarft penetrative distance from the lower interface ! ! of 'kpen' layer. Note that bogbot & bogtop at the 'kpen' layer either when fer ! ! is zero or non-zero was already calculated above. ! ! It seems that below qudarature solving formula is valid only when bogbot < 0. ! ! Below solving equation is clearly wrong ! I should revise this ! ! ! ------------------------------------------------------------------------------ ! if( drage .eq. 0._r8 ) then aquad = ( bogtop - bogbot ) / ( pifc0(kpen) - pifc0(kpen-1) ) bquad = 2._r8 * bogbot cquad = -wu(kpen-1)**2 * rhomid0j call roots(aquad,bquad,cquad,xc1,xc2,status) if( status .eq. 0 ) then if( xc1 .le. 0._r8 .and. xc2 .le. 0._r8 ) then ppen = max( xc1, xc2 ) ppen = min( 0._r8,max( -dp0(kpen), ppen ) ) elseif( xc1 .gt. 0._r8 .and. xc2 .gt. 0._r8 ) then ppen = -dp0(kpen) if (scverbose) then call write_parallel('Warning : UW-Cumulus penetrates up to kpen interface') end if else ppen = min( xc1, xc2 ) ppen = min( 0._r8,max( -dp0(kpen), ppen ) ) endif else ppen = -dp0(kpen) if (scverbose) then call write_parallel('Warning : UW-Cumulus penetrates up to kpen interface') end if endif else ppen = compute_ppen(wtwb,drage,bogbot,bogtop,rhomid0j,dp0(kpen)) endif if( ppen .eq. -dp0(kpen) .or. ppen .eq. 0._r8 ) limit_ppen(i) = 1. ! -------------------------------------------------------------------- ! ! Re-calculate the amount of expelled condensate from cloud updraft ! ! at the cumulus top. This is necessary for refined calculations of ! ! bulk cloud microphysics at the cumulus top. Note that ppen < 0. ! ! In the below, I explicitly calculate 'thlu_top' & 'qtu_top' by ! ! using non-zero 'fer(kpen)'. ! ! -------------------------------------------------------------------- ! if( fer(kpen)*(-ppen) .lt. 1.e-4_r8 ) then thlu_top = thlu(kpen-1) + ( thl0(kpen) + ssthl0(kpen) * (-ppen) / 2._r8 - thlu(kpen-1) ) * fer(kpen) * (-ppen) qtu_top = qtu(kpen-1) + ( qt0(kpen) + ssqt0(kpen) * (-ppen) / 2._r8 - qtu(kpen-1) ) * fer(kpen) * (-ppen) else thlu_top = ( thl0(kpen) + ssthl0(kpen) / fer(kpen) - ssthl0(kpen) * (-ppen) / 2._r8 ) - & ( thl0(kpen) + ssthl0(kpen) * (-ppen) / 2._r8 - thlu(kpen-1) + ssthl0(kpen) / fer(kpen) ) & * exp(-fer(kpen) * (-ppen)) qtu_top = ( qt0(kpen) + ssqt0(kpen) / fer(kpen) - ssqt0(kpen) * (-ppen) / 2._r8 ) - & ( qt0(kpen) + ssqt0(kpen) * (-ppen) / 2._r8 - qtu(kpen-1) + ssqt0(kpen) / fer(kpen) ) & * exp(-fer(kpen) * (-ppen)) end if call conden(pifc0(kpen-1)+ppen,thlu_top,qtu_top,thj,qvj,qlj,qij,qse,id_check) if( id_check .eq. 1 ) then exit_conden(i) = 1. id_exit = .true. if (scverbose) then call write_parallel('------- UW ShCu: exit, conden') end if go to 333 end if exntop = ((pifc0(kpen-1)+ppen)/p00)**rovcp if( (qlj + qij) .gt. criqc ) then dwten(kpen) = ( ( qlj + qij ) - criqc ) * qlj / ( qlj + qij ) diten(kpen) = ( ( qlj + qij ) - criqc ) * qij / ( qlj + qij ) qtu_top = qtu_top - dwten(kpen) - diten(kpen) thlu_top = thlu_top + (xlv/cp/exntop)*dwten(kpen) + (xls/cp/exntop)*diten(kpen) else dwten(kpen) = 0._r8 diten(kpen) = 0._r8 endif ! ----------------------------------------------------------------------- ! ! Calculate cumulus scale height as the top height that cumulus can reach.! ! ----------------------------------------------------------------------- ! rhoifc0j = pifc0(kpen-1)/(r*0.5_r8*(thv0bot(kpen)+thv0top(kpen-1))*exnifc0(kpen-1)) cush = zifc0(kpen-1) - ppen/rhoifc0j/g scaleh = cush end do ! iter_scaleh loop ! -------------------------------------------------------------------- ! ! The 'forcedCu' is logical identifier saying whether cumulus updraft ! ! overcome the buoyancy barrier just above the PBL top. If it is true, ! ! cumulus did not overcome the barrier - this is a shallow convection ! ! with negative cloud buoyancy, mimicking shallow continental cumulus ! ! convection. Depending on 'forcedCu' parameter, treatment of heat & ! ! moisture fluxes at the entraining interfaces, 'kbup <= k < kpen - 1' ! ! will be set up in a different ways, as will be shown later. ! ! -------------------------------------------------------------------- ! if( kbup .eq. krel ) then forcedCu = .true. limit_shcu(i) = 1. else forcedCu = .false. limit_shcu(i) = 0. endif ! ------------------------------------------------------------------ ! ! Filtering of unerasonable cumulus adjustment here. This is a very ! ! important process which should be done cautiously. Various ways of ! ! filtering are possible depending on cases mainly using the indices ! ! of key layers - 'klcl','kinv','krel','klfc','kbup','kpen'. At this ! ! stage, the followings are all possible : 'kinv >= 2', 'klcl >= 1', ! ! 'krel >= kinv', 'kbup >= krel', 'kpen >= krel'. I must design this ! ! filtering very cautiously, in such that none of realistic cumulus ! ! convection is arbitrarily turned-off. Potentially, I might turn-off! ! cumulus convection if layer-mean 'ql > 0' in the 'kinv-1' layer,in ! ! order to suppress cumulus convection growing, based at the Sc top. ! ! This is one of potential future modifications. Note that ppen < 0. ! ! ------------------------------------------------------------------ ! cldhgt = pifc0(kpen-1) + ppen if( forcedCu ) then if (scverbose) then call write_parallel( 'forcedCu - did not overcome initial buoyancy barrier') end if exit_cufilter(i) = 1. id_exit = .true. go to 333 end if ! Limit 'additional shallow cumulus' for DYCOMS simulation. ! if( cldhgt.ge.88000._r8 ) then ! id_exit = .true. ! go to 333 ! end if ! ------------------------------------------------------------------------------ ! ! Re-initializing some key variables above the 'kpen' layer in order to suppress ! ! the influence of non-physical values above 'kpen', in association with the use ! ! of 'iter_scaleh' loop. Note that umf, emf, ufrc are defined at the interfaces ! ! (0:k0), while 'dwten','diten', 'fer', 'fdr' are defined at layer mid-points. ! ! Initialization of 'fer' and 'fdr' is for correct writing purpose of diagnostic ! ! output. Note that we set umf(kpen)=emf(kpen)=ufrc(kpen)=0, in consistent with ! ! wtw < 0 at the top interface of 'kpen' layer. However, we still have non-zero ! ! expelled cloud condensate in the 'kpen' layer. ! ! ------------------------------------------------------------------------------ ! umf(kpen:k0) = 0._r8 emf(kpen:k0) = 0._r8 ufrc(kpen:k0) = 0._r8 dwten(kpen+1:k0) = 0._r8 diten(kpen+1:k0) = 0._r8 fer(kpen+1:k0) = 0._r8 fdr(kpen+1:k0) = 0._r8 xco(kpen+1:k0) = 0._r8 ! ------------------------------------------------------------------------ ! ! Calculate downward penetrative entrainment mass flux, 'emf(k) < 0', and ! ! thermodynamic properties of penetratively entrained airs at entraining ! ! interfaces. emf(k) is defined from the top interface of the layer kbup ! ! to the bottom interface of the layer 'kpen'. Note even when kbup = krel,! ! i.e.,even when 'kbup' was not updated in the above buoyancy sorting do ! ! loop (i.e., 'kbup' remains as the initialization value), below do loop ! ! of penetrative entrainment flux can be performed without any conceptual ! ! or logical problems, because we have already computed all the variables ! ! necessary for performing below penetrative entrainment block. ! ! In the below 'do' loop, 'k' is an interface index at which non-zero 'emf'! ! (penetrative entrainment mass flux) is calculated. Since cumulus updraft ! ! is negatively buoyant in the layers between the top interface of 'kbup' ! ! layer (interface index, kbup) and the top interface of 'kpen' layer, the ! ! fractional lateral entrainment, fer(k) within these layers will be close ! ! to zero - so it is likely that only strong lateral detrainment occurs in ! ! thses layers. Under this situation,we can easily calculate the amount of ! ! detrainment cumulus air into these negatively buoyanct layers by simply ! ! comparing cumulus updraft mass fluxes between the base and top interface ! ! of each layer: emf(k) = emf(k-1)*exp(-fdr(k)*dp0(k)) ! ! ~ emf(k-1)*(1-rei(k)*dp0(k)) ! ! emf(k-1)-emf(k) ~ emf(k-1)*rei(k)*dp0(k) ! ! Current code assumes that about 'rpen~10' times of these detrained mass ! ! are penetratively re-entrained down into the 'k-1' interface. And all of ! ! these detrained masses are finally dumped down into the top interface of ! ! 'kbup' layer. Thus, the amount of penetratively entrained air across the ! ! top interface of 'kbup' layer with 'rpen~10' becomes too large. ! ! Note that this penetrative entrainment part can be completely turned-off ! ! and we can simply use normal buoyancy-sorting involved turbulent fluxes ! ! by modifying 'penetrative entrainment fluxes' part below. ! ! ------------------------------------------------------------------------ ! ! -----------------------------------------------------------------------! ! Calculate entrainment mass flux and conservative scalars of entraining ! ! free air at interfaces of 'kbup <= k < kpen - 1' ! ! ---------------------------------------------------------------------- ! do k = 0, k0 thlu_emf(k) = thlu(k) qtu_emf(k) = qtu(k) uu_emf(k) = uu(k) vu_emf(k) = vu(k) do m = 1, ncnst tru_emf(k,m) = tru(k,m) enddo end do do k = kpen - 1, kbup, -1 ! Here, 'k' is an interface index at which ! penetrative entrainment fluxes are calculated. rhoifc0j = pifc0(k) / ( r * 0.5_r8 * ( thv0bot(k+1) + thv0top(k) ) * exnifc0(k) ) if( k .eq. kpen - 1 ) then ! ------------------------------------------------------------------------ ! ! Note that 'ppen' has already been calculated in the above 'iter_scaleh' ! ! loop assuming zero lateral entrainmentin the layer 'kpen'. ! ! ------------------------------------------------------------------------ ! ! -------------------------------------------------------------------- ! ! Calculate returning mass flux, emf ( < 0 ) ! ! Current penetrative entrainment rate with 'rpen~10' is too large and ! ! future refinement is necessary including the definition of 'thl','qt'! ! of penetratively entrained air. Penetratively entrained airs across ! ! the 'kpen-1' interface is assumed to have the properties of the base ! ! interface of 'kpen' layer. Note that 'emf ~ - umf/ufrc = - w * rho'. ! ! Thus, below limit sets an upper limit of |emf| to be ~ 10cm/s, which ! ! is very loose constraint. Here, I used more restricted constraint on ! ! the limit of emf, assuming 'emf' cannot exceed a net mass within the ! ! layer above the interface. Similar to the case of warming and drying ! ! due to cumulus updraft induced compensating subsidence, penetrative ! ! entrainment induces compensating upwelling - in order to prevent ! ! numerical instability in association with compensating upwelling, we ! ! should similarily limit the amount of penetrative entrainment at the ! ! interface by the amount of masses within the layer just above the ! ! penetratively entraining interface. ! ! -------------------------------------------------------------------- ! if( ( umf(k)*ppen*rei(kpen)*rpen ) .lt. -0.1_r8*rhoifc0j ) limit_emf(i) = 1. if( ( umf(k)*ppen*rei(kpen)*rpen ) .lt. -0.9_r8*dp0(kpen)/g/dt ) limit_emf(i) = 1. emf(k) = max( max( umf(k)*ppen*rei(kpen)*rpen, -0.1_r8*rhoifc0j), -0.9_r8*dp0(kpen)/g/dt) thlu_emf(k) = thl0(kpen) + ssthl0(kpen) * ( pifc0(k) - pmid0(kpen) ) qtu_emf(k) = qt0(kpen) + ssqt0(kpen) * ( pifc0(k) - pmid0(kpen) ) uu_emf(k) = u0(kpen) + ssu0(kpen) * ( pifc0(k) - pmid0(kpen) ) vu_emf(k) = v0(kpen) + ssv0(kpen) * ( pifc0(k) - pmid0(kpen) ) do m = 1, ncnst tru_emf(k,m) = tr0(kpen,m) + sstr0(kpen,m) * ( pifc0(k) - pmid0(kpen) ) enddo else ! if(k.lt.kpen-1). ! --------------------------------------------------------------------------- ! ! Note we are coming down from the higher interfaces to the lower interfaces. ! ! Also note that 'emf < 0'. So, below operation is a summing not subtracting. ! ! In order to ensure numerical stability, I imposed a modified correct limit ! ! of '-0.9*dp0(k+1)/g/dt' on emf(k). ! ! --------------------------------------------------------------------------- ! if( use_cumpenent ) then ! Original Cumulative Penetrative Entrainment if( ( emf(k+1)-umf(k)*dp0(k+1)*rei(k+1)*rpen ) .lt. -0.1_r8*rhoifc0j ) limit_emf(i) = 1 if( ( emf(k+1)-umf(k)*dp0(k+1)*rei(k+1)*rpen ) .lt. -0.9_r8*dp0(k+1)/g/dt ) limit_emf(i) = 1 emf(k) = max(max(emf(k+1)-umf(k)*dp0(k+1)*rei(k+1)*rpen, -0.1_r8*rhoifc0j), -0.9_r8*dp0(k+1)/g/dt ) if( abs(emf(k)) .gt. abs(emf(k+1)) ) then thlu_emf(k) = ( thlu_emf(k+1) * emf(k+1) + thl0(k+1) * ( emf(k) - emf(k+1) ) ) / emf(k) qtu_emf(k) = ( qtu_emf(k+1) * emf(k+1) + qt0(k+1) * ( emf(k) - emf(k+1) ) ) / emf(k) uu_emf(k) = ( uu_emf(k+1) * emf(k+1) + u0(k+1) * ( emf(k) - emf(k+1) ) ) / emf(k) vu_emf(k) = ( vu_emf(k+1) * emf(k+1) + v0(k+1) * ( emf(k) - emf(k+1) ) ) / emf(k) do m = 1, ncnst tru_emf(k,m) = ( tru_emf(k+1,m) * emf(k+1) + tr0(k+1,m) * ( emf(k) - emf(k+1) ) ) / emf(k) enddo else thlu_emf(k) = thl0(k+1) qtu_emf(k) = qt0(k+1) uu_emf(k) = u0(k+1) vu_emf(k) = v0(k+1) do m = 1, ncnst tru_emf(k,m) = tr0(k+1,m) enddo endif else ! Alternative Non-Cumulative Penetrative Entrainment if( ( -umf(k)*dp0(k+1)*rei(k+1)*rpen ) .lt. -0.1_r8*rhoifc0j ) limit_emf(i) = 1 if( ( -umf(k)*dp0(k+1)*rei(k+1)*rpen ) .lt. -0.9_r8*dp0(k+1)/g/dt ) limit_emf(i) = 1 emf(k) = max(max(-umf(k)*dp0(k+1)*rei(k+1)*rpen, -0.1_r8*rhoifc0j), -0.9_r8*dp0(k+1)/g/dt ) thlu_emf(k) = thl0(k+1) qtu_emf(k) = qt0(k+1) uu_emf(k) = u0(k+1) vu_emf(k) = v0(k+1) do m = 1, ncnst tru_emf(k,m) = tr0(k+1,m) enddo endif endif ! ---------------------------------------------------------------------------- ! ! In this GCM modeling framework, all what we should do is to calculate heat ! ! and moisture fluxes at the given geometrically-fixed height interfaces - we ! ! don't need to worry about movement of material height surface in association ! ! with compensating subsidence or unwelling, in contrast to the bulk modeling. ! ! In this geometrically fixed height coordinate system, heat and moisture flux ! ! at the geometrically fixed height handle everything - a movement of material ! ! surface is implicitly treated automatically. Note that in terms of turbulent ! ! heat and moisture fluxes at model interfaces, both the cumulus updraft mass ! ! flux and penetratively entraining mass flux play the same role -both of them ! ! warms and dries the 'kbup' layer, cools and moistens the 'kpen' layer, and ! ! cools and moistens any intervening layers between 'kbup' and 'kpen' layers. ! ! It is important to note these identical roles on turbulent heat and moisture ! ! fluxes of 'umf' and 'emf'. ! ! When 'kbup' is a stratocumulus-topped PBL top interface, increase of 'rpen' ! ! is likely to strongly diffuse stratocumulus top interface, resulting in the ! ! reduction of cloud fraction. In this sense, the 'kbup' interface has a very ! ! important meaning and role : across the 'kbup' interface, strong penetrative ! ! entrainment occurs, thus any sharp gradient properties across that interface ! ! are easily diffused through strong mass exchange. Thus, an initialization of ! ! 'kbup' (and also 'kpen') should be done very cautiously as mentioned before. ! ! In order to prevent this stron diffusion for the shallow cumulus convection ! ! based at the Sc top, it seems to be good to initialize 'kbup = krel', rather ! ! that 'kbup = krel-1'. ! ! ---------------------------------------------------------------------------- ! end do ! k loop !------------------------------------------------------------------ ! ! ! ! Compute turbulent heat, moisture, momentum flux at all interfaces ! ! ! !------------------------------------------------------------------ ! ! It is very important to note that in calculating turbulent fluxes ! ! below, we must not double count turbulent flux at any interefaces.! ! In the below, turbulent fluxes at the interfaces (interface index ! ! k) are calculated by the following 4 blocks in consecutive order: ! ! ! ! (1) " 0 <= k <= kinv - 1 " : PBL fluxes. ! ! From 'fluxbelowinv' using reconstructed PBL height. Currently,! ! the reconstructed PBLs are independently calculated for each ! ! individual conservative scalar variables ( qt, thl, u, v ) in ! ! each 'fluxbelowinv', instead of being uniquely calculated by ! ! using thvl. Turbulent flux at the surface is assumed to be 0. ! ! (2) " kinv <= k <= krel - 1 " : Non-buoyancy sorting fluxes ! ! Assuming cumulus mass flux and cumulus updraft thermodynamic ! ! properties (except u, v which are modified by the PGFc during ! ! upward motion) are conserved during a updraft motion from the ! ! PBL top interface to the release level. If these layers don't ! ! exist (e,g, when 'krel = kinv'), then current routine do not ! ! perform this routine automatically. So I don't need to modify ! ! anything. ! ! (3) " krel <= k <= kbup - 1 " : Buoyancy sorting fluxes ! ! From laterally entraining-detraining buoyancy sorting plumes. ! ! (4) " kbup <= k < kpen-1 " : Penetrative entrainment fluxes ! ! From penetratively entraining plumes, ! ! ! ! In case of normal situation, turbulent interfaces in each groups ! ! are mutually independent of each other. Thus double flux counting ! ! or ambiguous flux counting requiring the choice among the above 4 ! ! groups do not occur normally. However, in case that cumulus plume ! ! could not completely overcome the buoyancy barrier just above the ! ! PBL top interface and so 'kbup = krel' (.forcedCu=.true.) ( here, ! ! it can be either 'kpen = krel' as the initialization, or ' kpen > ! ! krel' if cumulus updraft just penetrated over the top of release ! ! layer ). If this happens, we should be very careful in organizing ! ! the sequence of the 4 calculation routines above - note that the ! ! routine located at the later has the higher priority. Additional ! ! feature I must consider is that when 'kbup = kinv - 1' (this is a ! ! combined situation of 'kbup=krel-1' & 'krel = kinv' when I chose ! ! 'kbup=krel-1' instead of current choice of 'kbup=krel'), a strong ! ! penetrative entrainment fluxes exists at the PBL top interface, & ! ! all of these fluxes are concentrated (deposited) within the layer ! ! just below PBL top interface (i.e., 'kinv-1' layer). On the other ! ! hand, in case of 'fluxbelowinv', only the compensating subsidence ! ! effect is concentrated in the 'kinv-1' layer and 'pure' turbulent ! ! heat and moisture fluxes ( 'pure' means the fluxes not associated ! ! with compensating subsidence) are linearly distributed throughout ! ! the whole PBL. Thus different choice of the above flux groups can ! ! produce very different results. Output variable should be written ! ! consistently to the choice of computation sequences. ! ! When the case of 'kbup = krel(-1)' happens,another way to dealing ! ! with this case is to simply ' exit ' the whole cumulus convection ! ! calculation without performing any cumulus convection. We can ! ! choose this approach by specifying a condition in the 'Filtering ! ! of unreasonable cumulus adjustment' just after 'iter_scaleh'. But ! ! this seems not to be a good choice (although this choice was used ! ! previous code ), since it might arbitrary damped-out the shallow ! ! cumulus convection over the continent land, where shallow cumulus ! ! convection tends to be negatively buoyant. ! ! ----------------------------------------------------------------- ! ! --------------------------------------------------- ! ! 1. PBL fluxes : 0 <= k <= kinv - 1 ! ! All the information necessary to reconstruct PBL ! ! height are passed to 'fluxbelowinv'. ! ! --------------------------------------------------- ! xsrc = qtsrc xmean = qt0(kinv) xtop = qt0(kinv+1) + ssqt0(kinv+1) * ( pifc0(kinv) - pmid0(kinv+1) ) xbot = qt0(kinv-1) + ssqt0(kinv-1) * ( pifc0(kinv-1) - pmid0(kinv-1) ) call fluxbelowinv( cbmf, pifc0(0:k0), k0, kinv, dt, xsrc, xmean, xtop, xbot, xflx ) qtflx(0:kinv-1) = xflx(0:kinv-1) xsrc = thlsrc xmean = thl0(kinv) xtop = thl0(kinv+1) + ssthl0(kinv+1) * ( pifc0(kinv) - pmid0(kinv+1) ) xbot = thl0(kinv-1) + ssthl0(kinv-1) * ( pifc0(kinv-1) - pmid0(kinv-1) ) call fluxbelowinv( cbmf, pifc0(0:k0), k0, kinv, dt, xsrc, xmean, xtop, xbot, xflx ) slflx(0:kinv-1) = cp * exnifc0(0:kinv-1) * xflx(0:kinv-1) xsrc = usrc xmean = u0(kinv) xtop = u0(kinv+1) + ssu0(kinv+1) * ( pifc0(kinv) - pmid0(kinv+1) ) xbot = u0(kinv-1) + ssu0(kinv-1) * ( pifc0(kinv-1) - pmid0(kinv-1) ) call fluxbelowinv( cbmf, pifc0(0:k0), k0, kinv, dt, xsrc, xmean, xtop, xbot, xflx ) uflx(0:kinv-1) = xflx(0:kinv-1) xsrc = vsrc xmean = v0(kinv) xtop = v0(kinv+1) + ssv0(kinv+1) * ( pifc0(kinv) - pmid0(kinv+1) ) xbot = v0(kinv-1) + ssv0(kinv-1) * ( pifc0(kinv-1) - pmid0(kinv-1) ) call fluxbelowinv( cbmf, pifc0(0:k0), k0, kinv, dt, xsrc, xmean, xtop, xbot, xflx ) vflx(0:kinv-1) = xflx(0:kinv-1) do m = 1, ncnst xsrc = trsrc(m) xmean = tr0(kinv,m) xtop = tr0(kinv+1,m) + sstr0(kinv+1,m) * ( pifc0(kinv) - pmid0(kinv+1) ) xbot = tr0(kinv-1,m) + sstr0(kinv-1,m) * ( pifc0(kinv-1) - pmid0(kinv-1) ) call fluxbelowinv( cbmf, pifc0(0:k0), k0, kinv, dt, xsrc, xmean, xtop, xbot, xflx ) trflx(0:kinv-1,m) = xflx(0:kinv-1) enddo ! -------------------------------------------------------------- ! ! 2. Non-buoyancy sorting fluxes : kinv <= k <= krel - 1 ! ! Note that when 'krel = kinv', below block is never executed ! ! as in a desirable, expected way ( but I must check if this ! ! is the case ). The non-buoyancy sorting fluxes are computed ! ! only when 'krel > kinv'. ! ! -------------------------------------------------------------- ! uplus = 0._r8 vplus = 0._r8 do k = kinv, krel - 1 kp1 = k + 1 qtflx(k) = cbmf * ( qtsrc - ( qt0(kp1) + ssqt0(kp1) * ( pifc0(k) - pmid0(kp1) ) ) ) slflx(k) = cbmf * ( thlsrc - ( thl0(kp1) + ssthl0(kp1) * ( pifc0(k) - pmid0(kp1) ) ) ) * cp * exnifc0(k) uplus = uplus + PGFc * ssu0(k) * ( pifc0(k) - pifc0(k-1) ) vplus = vplus + PGFc * ssv0(k) * ( pifc0(k) - pifc0(k-1) ) uflx(k) = cbmf * ( usrc + uplus - ( u0(kp1) + ssu0(kp1) * ( pifc0(k) - pmid0(kp1) ) ) ) vflx(k) = cbmf * ( vsrc + vplus - ( v0(kp1) + ssv0(kp1) * ( pifc0(k) - pmid0(kp1) ) ) ) do m = 1, ncnst trflx(k,m) = cbmf * ( trsrc(m) - ( tr0(kp1,m) + sstr0(kp1,m) * ( pifc0(k) - pmid0(kp1) ) ) ) enddo end do ! ------------------------------------------------------------------------ ! ! 3. Buoyancy sorting fluxes : krel <= k <= kbup - 1 ! ! In case that 'kbup = krel - 1 ' ( or even in case 'kbup = krel' ), ! ! buoyancy sorting fluxes are not calculated, which is consistent, ! ! desirable feature. ! ! ------------------------------------------------------------------------ ! do k = krel, kbup - 1 kp1 = k + 1 slflx(k) = cp * exnifc0(k) * umf(k) * ( thlu(k) - ( thl0(kp1) + ssthl0(kp1) * ( pifc0(k) - pmid0(kp1) ) ) ) qtflx(k) = umf(k) * ( qtu(k) - ( qt0(kp1) + ssqt0(kp1) * ( pifc0(k) - pmid0(kp1) ) ) ) uflx(k) = umf(k) * ( uu(k) - ( u0(kp1) + ssu0(kp1) * ( pifc0(k) - pmid0(kp1) ) ) ) vflx(k) = umf(k) * ( vu(k) - ( v0(kp1) + ssv0(kp1) * ( pifc0(k) - pmid0(kp1) ) ) ) do m = 1, ncnst trflx(k,m) = umf(k) * ( tru(k,m) - ( tr0(kp1,m) + sstr0(kp1,m) * ( pifc0(k) - pmid0(kp1) ) ) ) enddo end do ! ------------------------------------------------------------------------- ! ! 4. Penetrative entrainment fluxes : kbup <= k <= kpen - 1 ! ! The only confliction that can happen is when 'kbup = kinv-1'. For this ! ! case, turbulent flux at kinv-1 is calculated both from 'fluxbelowinv' ! ! and here as penetrative entrainment fluxes. Since penetrative flux is ! ! calculated later, flux at 'kinv - 1 ' will be that of penetrative flux.! ! However, turbulent flux calculated at 'kinv - 1' from penetrative entr.! ! is less attractable, since more reasonable turbulent flux at 'kinv-1' ! ! should be obtained from 'fluxbelowinv', by considering re-constructed ! ! inversion base height. This conflicting problem can be solved if we can! ! initialize 'kbup = krel', instead of kbup = krel - 1. This choice seems! ! to be more reasonable since it is not conflicted with 'fluxbelowinv' in! ! calculating fluxes at 'kinv - 1' ( for this case, flux at 'kinv-1' is ! ! always from 'fluxbelowinv' ), and flux at 'krel-1' is calculated from ! ! the non-buoyancy sorting flux without being competed with penetrative ! ! entrainment fluxes. Even when we use normal cumulus flux instead of ! ! penetrative entrainment fluxes at 'kbup <= k <= kpen-1' interfaces, ! ! the initialization of kbup=krel perfectly works without any conceptual ! ! confliction. Thus it seems to be much better to choose 'kbup = krel' ! ! initialization of 'kbup', which is current choice. ! ! Note that below formula uses conventional updraft cumulus fluxes for ! ! shallow cumulus which did not overcome the first buoyancy barrier above! ! PBL top while uses penetrative entrainment fluxes for the other cases ! ! 'kbup <= k <= kpen-1' interfaces. Depending on cases, however, I can ! ! selelct different choice. ! ! ------------------------------------------------------------------------------------------------------------------ ! ! if( forcedCu ) then ! ! slflx(k) = cp * exns0(k) * umf(k) * ( thlu(k) - ( thl0(kp1) + ssthl0(kp1) * ( ps0(k) - p0(kp1) ) ) ) ! ! qtflx(k) = umf(k) * ( qtu(k) - ( qt0(kp1) + ssqt0(kp1) * ( ps0(k) - p0(kp1) ) ) ) ! ! uflx(k) = umf(k) * ( uu(k) - ( u0(kp1) + ssu0(kp1) * ( ps0(k) - p0(kp1) ) ) ) ! ! vflx(k) = umf(k) * ( vu(k) - ( v0(kp1) + ssv0(kp1) * ( ps0(k) - p0(kp1) ) ) ) ! ! do m = 1, ncnst ! ! trflx(k,m) = umf(k) * ( tru(k,m) - ( tr0(kp1,m) + sstr0(kp1,m) * ( ps0(k) - p0(kp1) ) ) ) ! ! enddo ! ! else ! ! slflx(k) = cp * exns0(k) * emf(k) * ( thlu_emf(k) - ( thl0(k) + ssthl0(k) * ( ps0(k) - p0(k) ) ) ) ! ! qtflx(k) = emf(k) * ( qtu_emf(k) - ( qt0(k) + ssqt0(k) * ( ps0(k) - p0(k) ) ) ) ! ! uflx(k) = emf(k) * ( uu_emf(k) - ( u0(k) + ssu0(k) * ( ps0(k) - p0(k) ) ) ) ! ! vflx(k) = emf(k) * ( vu_emf(k) - ( v0(k) + ssv0(k) * ( ps0(k) - p0(k) ) ) ) ! ! do m = 1, ncnst ! ! trflx(k,m) = emf(k) * ( tru_emf(k,m) - ( tr0(k,m) + sstr0(k,m) * ( ps0(k) - p0(k) ) ) ) ! ! enddo ! ! endif ! ! ! ! if( use_uppenent ) then ! Combined Updraft + Penetrative Entrainment Flux ! ! slflx(k) = cp * exns0(k) * umf(k) * ( thlu(k) - ( thl0(kp1) + ssthl0(kp1) * ( ps0(k) - p0(kp1) ) ) ) + & ! ! cp * exns0(k) * emf(k) * ( thlu_emf(k) - ( thl0(k) + ssthl0(k) * ( ps0(k) - p0(k) ) ) ) ! ! qtflx(k) = umf(k) * ( qtu(k) - ( qt0(kp1) + ssqt0(kp1) * ( ps0(k) - p0(kp1) ) ) ) + & ! ! emf(k) * ( qtu_emf(k) - ( qt0(k) + ssqt0(k) * ( ps0(k) - p0(k) ) ) ) ! ! uflx(k) = umf(k) * ( uu(k) - ( u0(kp1) + ssu0(kp1) * ( ps0(k) - p0(kp1) ) ) ) + & ! ! emf(k) * ( uu_emf(k) - ( u0(k) + ssu0(k) * ( ps0(k) - p0(k) ) ) ) ! ! vflx(k) = umf(k) * ( vu(k) - ( v0(kp1) + ssv0(kp1) * ( ps0(k) - p0(kp1) ) ) ) + & ! ! emf(k) * ( vu_emf(k) - ( v0(k) + ssv0(k) * ( ps0(k) - p0(k) ) ) ) ! ! do m = 1, ncnst ! ! trflx(k,m) = umf(k) * ( tru(k,m) - ( tr0(kp1,m) + sstr0(kp1,m) * ( ps0(k) - p0(kp1) ) ) ) + & ! ! emf(k) * ( tru_emf(k,m) - ( tr0(k,m) + sstr0(k,m) * ( ps0(k) - p0(k) ) ) ) ! ! enddo ! ! ------------------------------------------------------------------------------------------------------------------ ! do k = kbup, kpen - 1 kp1 = k + 1 slflx(k) = cp * exnifc0(k) * emf(k) * ( thlu_emf(k) - ( thl0(k) + ssthl0(k) * ( pifc0(k) - pmid0(k) ) ) ) qtflx(k) = emf(k) * ( qtu_emf(k) - ( qt0(k) + ssqt0(k) * ( pifc0(k) - pmid0(k) ) ) ) uflx(k) = emf(k) * ( uu_emf(k) - ( u0(k) + ssu0(k) * ( pifc0(k) - pmid0(k) ) ) ) vflx(k) = emf(k) * ( vu_emf(k) - ( v0(k) + ssv0(k) * ( pifc0(k) - pmid0(k) ) ) ) do m = 1, ncnst trflx(k,m) = emf(k) * ( tru_emf(k,m) - ( tr0(k,m) + sstr0(k,m) * ( pifc0(k) - pmid0(k) ) ) ) enddo end do ! ------------------------------------------- ! ! Turn-off cumulus momentum flux as an option ! ! ------------------------------------------- ! if( .not. use_momenflx ) then uflx(0:k0) = 0._r8 vflx(0:k0) = 0._r8 endif ! -------------------------------------------------------- ! ! Condensate tendency by compensating subsidence/upwelling ! ! -------------------------------------------------------- ! uemf(0:k0) = 0._r8 do k = 0, kinv - 2 ! Assume linear updraft mass flux within the PBL. uemf(k) = cbmf * ( pifc0(0) - pifc0(k) ) / ( pifc0(0) - pifc0(kinv-1) ) end do uemf(kinv-1:krel-1) = cbmf uemf(krel:kbup-1) = umf(krel:kbup-1) uemf(kbup:kpen-1) = emf(kbup:kpen-1) ! Only use penetrative entrainment flux consistently. comsub(1:k0) = 0._r8 do k = 1, kpen comsub(k) = 0.5 * ( uemf(k) + uemf(k-1) ) ! comsub defined on interfaces end do do k = 1, kpen if( comsub(k) .ge. 0._r8 ) then if( k .eq. k0 ) then thlten_sub = 0._r8 qtten_sub = 0._r8 qlten_sub = 0._r8 qiten_sub = 0._r8 nlten_sub = 0._r8 niten_sub = 0._r8 else thlten_sub = g * comsub(k) * ( thl0(k+1) - thl0(k) ) / ( pmid0(k) - pmid0(k+1) ) qtten_sub = g * comsub(k) * ( qt0(k+1) - qt0(k) ) / ( pmid0(k) - pmid0(k+1) ) qlten_sub = g * comsub(k) * ( ql0(k+1) - ql0(k) ) / ( pmid0(k) - pmid0(k+1) ) qiten_sub = g * comsub(k) * ( qi0(k+1) - qi0(k) ) / ( pmid0(k) - pmid0(k+1) ) ! nlten_sub = g * comsub(k) * ( tr0(k+1,ixnumliq) - tr0(k,ixnumliq) ) / ( pmid0(k) - pmid0(k+1) ) ! niten_sub = g * comsub(k) * ( tr0(k+1,ixnumice) - tr0(k,ixnumice) ) / ( pmid0(k) - pmid0(k+1) ) endif else if( k .eq. 1 ) then thlten_sub = 0._r8 qtten_sub = 0._r8 qlten_sub = 0._r8 qiten_sub = 0._r8 nlten_sub = 0._r8 niten_sub = 0._r8 else thlten_sub = g * comsub(k) * ( thl0(k) - thl0(k-1) ) / ( pmid0(k-1) - pmid0(k) ) qtten_sub = g * comsub(k) * ( qt0(k) - qt0(k-1) ) / ( pmid0(k-1) - pmid0(k) ) qlten_sub = g * comsub(k) * ( ql0(k) - ql0(k-1) ) / ( pmid0(k-1) - pmid0(k) ) qiten_sub = g * comsub(k) * ( qi0(k) - qi0(k-1) ) / ( pmid0(k-1) - pmid0(k) ) ! nlten_sub = g * comsub(k) * ( tr0(k,ixnumliq) - tr0(k-1,ixnumliq) ) / ( pmid0(k-1) - pmid0(k) ) ! niten_sub = g * comsub(k) * ( tr0(k,ixnumice) - tr0(k-1,ixnumice) ) / ( pmid0(k-1) - pmid0(k) ) endif endif thl_prog = thl0(k) + thlten_sub * dt qt_prog = max( qt0(k) + qtten_sub * dt, 1.e-12_r8 ) call conden(pmid0(k),thl_prog,qt_prog,thj,qvj,qlj,qij,qse,id_check) if( id_check .eq. 1 ) then id_exit = .true. if (scverbose) then call write_parallel('------- UW ShCu: exit, conden L3302') end if go to 333 endif ! qlten_sink(k) = ( qlj - ql0(k) ) / dt ! qiten_sink(k) = ( qij - qi0(k) ) / dt qlten_sink(k) = max( qlten_sub, - ql0(k) / dt ) ! For consistency with prognostic macrophysics scheme qiten_sink(k) = max( qiten_sub, - qi0(k) / dt ) ! For consistency with prognostic macrophysics scheme ! nlten_sink(k) = max( nlten_sub, - tr0(k,ixnumliq) / dt ) ! niten_sink(k) = max( niten_sub, - tr0(k,ixnumice) / dt ) end do ! --------------------------------------------- ! ! ! ! Calculate convective tendencies at each layer ! ! ! ! --------------------------------------------- ! ! ----------------- ! ! Momentum tendency ! ! ----------------- ! do k = 1, kpen km1 = k - 1 uten(k) = ( uflx(km1) - uflx(k) ) * g / dp0(k) vten(k) = ( vflx(km1) - vflx(k) ) * g / dp0(k) uf(k) = u0(k) + uten(k) * dt vf(k) = v0(k) + vten(k) * dt ! do m = 1, ncnst ! trten(k,m) = ( trflx(km1,m) - trflx(k,m) ) * g / dp0(k) ! ! Limit trten(k,m) such that negative value is not developed. ! ! This limitation does not conserve grid-mean tracers and future ! ! refinement is required for tracer-conserving treatment. ! trten(k,m) = max(trten(k,m),-tr0(k,m)/dt) ! enddo end do ! ----------------------------------------------------------------- ! ! Tendencies of thermodynamic variables. ! ! This part requires a careful treatment of bulk cloud microphysics.! ! Relocations of 'precipitable condensates' either into the surface ! ! or into the tendency of 'krel' layer will be performed just after ! ! finishing the below 'do-loop'. ! ! ----------------------------------------------------------------- ! do k = 1, kpen km1 = k - 1 ! ------------------------------------------------------------------------------ ! ! Compute 'slten', 'qtten', 'qvten', 'qlten', 'qiten', and 'sten' ! ! ! ! Key assumptions made in this 'cumulus scheme' are : ! ! 1. Cumulus updraft expels condensate into the environment at the top interface ! ! of each layer. Note that in addition to this expel process ('source' term), ! ! cumulus updraft can modify layer mean condensate through normal detrainment ! ! forcing or compensating subsidence. ! ! 2. Expelled water can be either 'sustaining' or 'precipitating' condensate. By ! ! definition, 'suataining condensate' will remain in the layer where it was ! ! formed, while 'precipitating condensate' will fall across the base of the ! ! layer where it was formed. ! ! 3. All precipitating condensates are assumed to fall into the release layer or ! ! ground as soon as it was formed without being evaporated during the falling ! ! process down to the desinated layer ( either release layer of surface ). ! ! ------------------------------------------------------------------------------ ! ! ------------------------------------------------------------------------- ! ! 'dwten(k)','diten(k)' : Production rate of condensate within the layer k ! ! [ kg/kg/s ] by the expels of condensate from cumulus updraft. ! ! It is important to note that in terms of moisture tendency equation, this ! ! is a 'source' term of enviromental 'qt'. More importantly, these source ! ! are already counted in the turbulent heat and moisture fluxes we computed ! ! until now, assuming all the expelled condensate remain in the layer where ! ! it was formed. Thus, in calculation of 'qtten' and 'slten' below, we MUST ! ! NOT add or subtract these terms explicitly in order not to double or miss ! ! count, unless some expelled condensates fall down out of the layer. Note ! ! this falling-down process ( i.e., precipitation process ) and associated ! ! 'qtten' and 'slten' and production of surface precipitation flux will be ! ! treated later in 'zm_conv_evap' in 'convect_shallow_tend' subroutine. ! ! In below, we are converting expelled cloud condensate into correct unit. ! ! I found that below use of '0.5 * (umf(k-1) + umf(k))' causes conservation ! ! errors at some columns in global simulation. So, I returned to originals. ! ! This will cause no precipitation flux at 'kpen' layer since umf(kpen)=0. ! ! ------------------------------------------------------------------------- ! dwten(k) = dwten(k) * 0.5_r8 * ( umf(k-1) + umf(k) ) * g / dp0(k) ! [ kg/kg/s ] diten(k) = diten(k) * 0.5_r8 * ( umf(k-1) + umf(k) ) * g / dp0(k) ! [ kg/kg/s ] ! dwten(k) = dwten(k) * umf(k) * g / dp0(k) ! [ kg/kg/s ] ! diten(k) = diten(k) * umf(k) * g / dp0(k) ! [ kg/kg/s ] ! --------------------------------------------------------------------------- ! ! 'qrten(k)','qsten(k)' : Production rate of rain and snow within the layer k ! ! [ kg/kg/s ] by cumulus expels of condensates to the environment.! ! --------------------------------------------------------------------------- ! qrten(k) = frc_rasn * dwten(k) qsten(k) = frc_rasn * diten(k) ! ------------------------------------------------------------------------ ! ! 'slten(k)','qtten(k)' ! ! Note that 'slflx(k)' and 'qtflx(k)' we have calculated already included ! ! all the contributions of (1) expels of condensate (dwten(k), diten(k)), ! ! (2) mass detrainment ( delta * umf * ( qtu - qt ) ), & (3) compensating ! ! subsidence ( M * dqt / dz ). Thus 'slflx(k)' and 'qtflx(k)' we computed ! ! is a hybrid turbulent flux containing one part of 'source' term - expel ! ! of condensate. In order to calculate 'slten' and 'qtten', we should add ! ! additional 'source' term, if any. If the expelled condensate falls down ! ! across the base of the layer, it will be another sink (negative source) ! ! term. Note also that we included frictional heating terms in the below ! ! calculation of 'slten'. ! ! ------------------------------------------------------------------------ ! slten(k) = ( slflx(km1) - slflx(k) ) * g / dp0(k) if( k .eq. 1 ) then slten(k) = slten(k) - g / 4._r8 / dp0(k) * ( & uflx(k)*(uf(k+1) - uf(k) + u0(k+1) - u0(k)) + & vflx(k)*(vf(k+1) - vf(k) + v0(k+1) - v0(k))) elseif( k .ge. 2 .and. k .le. kpen-1 ) then slten(k) = slten(k) - g / 4._r8 / dp0(k) * ( & uflx(k)*(uf(k+1) - uf(k) + u0(k+1) - u0(k)) + & uflx(k-1)*(uf(k) - uf(k-1) + u0(k) - u0(k-1)) + & vflx(k)*(vf(k+1) - vf(k) + v0(k+1) - v0(k)) + & vflx(k-1)*(vf(k) - vf(k-1) + v0(k) - v0(k-1))) elseif( k .eq. kpen ) then slten(k) = slten(k) - g / 4._r8 / dp0(k) * ( & uflx(k-1)*(uf(k) - uf(k-1) + u0(k) - u0(k-1)) + & vflx(k-1)*(vf(k) - vf(k-1) + v0(k) - v0(k-1))) endif qtten(k) = ( qtflx(km1) - qtflx(k) ) * g / dp0(k) ! ---------------------------------------------------------------------------- ! ! Compute condensate tendency, including reserved condensate ! ! We assume that eventual detachment and detrainment occurs in kbup layer due ! ! to downdraft buoyancy sorting. In the layer above the kbup, only penetrative ! ! entrainment exists. Penetrative entrained air is assumed not to contain any ! ! condensate. ! ! ---------------------------------------------------------------------------- ! ! Compute in-cumulus condensate at the layer mid-point. if( k .lt. krel .or. k .gt. kpen ) then qlu_mid = 0._r8 qiu_mid = 0._r8 qlj = 0._r8 qij = 0._r8 elseif( k .eq. krel ) then call conden(prel,thlu(krel-1),qtu(krel-1),thj,qvj,qlj,qij,qse,id_check) if( id_check .eq. 1 ) then exit_conden(i) = 1. id_exit = .true. if (scverbose) then call write_parallel('------- UW ShCu: exit, conden') end if go to 333 endif qlubelow = qlj qiubelow = qij call conden(pifc0(k),thlu(k),qtu(k),thj,qvj,qlj,qij,qse,id_check) if( id_check .eq. 1 ) then exit_conden(i) = 1. id_exit = .true. if (scverbose) then call write_parallel('------- UW ShCu: exit, conden') end if go to 333 end if qlu_mid = 0.5_r8 * ( qlubelow + qlj ) * ( prel - pifc0(k) )/( pifc0(k-1) - pifc0(k) ) qiu_mid = 0.5_r8 * ( qiubelow + qij ) * ( prel - pifc0(k) )/( pifc0(k-1) - pifc0(k) ) elseif( k .eq. kpen ) then call conden(pifc0(k-1)+ppen,thlu_top,qtu_top,thj,qvj,qlj,qij,qse,id_check) if( id_check .eq. 1 ) then exit_conden(i) = 1. id_exit = .true. if (scverbose) then call write_parallel('------- UW ShCu: exit, conden') end if go to 333 end if qlu_mid = 0.5_r8 * ( qlubelow + qlj ) * ( -ppen ) /( pifc0(k-1) - pifc0(k) ) qiu_mid = 0.5_r8 * ( qiubelow + qij ) * ( -ppen ) /( pifc0(k-1) - pifc0(k) ) qlu_top = qlj qiu_top = qij else call conden(pifc0(k),thlu(k),qtu(k),thj,qvj,qlj,qij,qse,id_check) if( id_check .eq. 1 ) then exit_conden(i) = 1. id_exit = .true. if (scverbose) then call write_parallel('------- UW ShCu: exit, conden') end if go to 333 end if qlu_mid = 0.5_r8 * ( qlubelow + qlj ) qiu_mid = 0.5_r8 * ( qiubelow + qij ) endif qlubelow = qlj qiubelow = qij ! 1. Non-precipitating portion of expelled condensate qc_l(k) = ( 1._r8 - frc_rasn ) * dwten(k) ! [ kg/kg/s ] qc_i(k) = ( 1._r8 - frc_rasn ) * diten(k) ! [ kg/kg/s ] ! 2. Detrained Condensate if( k .le. kbup ) then qc_l(k) = qc_l(k) + g * 0.5_r8 * ( umf(k-1) + umf(k) ) * fdr(k) * qlu_mid ! [ kg/kg/s ] qc_i(k) = qc_i(k) + g * 0.5_r8 * ( umf(k-1) + umf(k) ) * fdr(k) * qiu_mid ! [ kg/kg/s ] qc_lm = - g * 0.5_r8 * ( umf(k-1) + umf(k) ) * fdr(k) * ql0(k) qc_im = - g * 0.5_r8 * ( umf(k-1) + umf(k) ) * fdr(k) * qi0(k) ! Below 'nc_lm', 'nc_im' should be used only when frc_rasn = 1. ! nc_lm = - g * 0.5_r8 * ( umf(k-1) + umf(k) ) * fdr(k) * tr0(k,ixnumliq) ! nc_im = - g * 0.5_r8 * ( umf(k-1) + umf(k) ) * fdr(k) * tr0(k,ixnumice) else qc_lm = 0._r8 qc_im = 0._r8 nc_lm = 0._r8 nc_im = 0._r8 endif ! 3. Detached Updraft if( k .eq. kbup ) then qc_l(k) = qc_l(k) + g * umf(k) * qlj / ( pifc0(k-1) - pifc0(k) ) ! [ kg/kg/s ] qc_i(k) = qc_i(k) + g * umf(k) * qij / ( pifc0(k-1) - pifc0(k) ) ! [ kg/kg/s ] qc_lm = qc_lm - g * umf(k) * ql0(k) / ( pifc0(k-1) - pifc0(k) ) ! [ kg/kg/s ] qc_im = qc_im - g * umf(k) * qi0(k) / ( pifc0(k-1) - pifc0(k) ) ! [ kg/kg/s ] ! nc_lm = nc_lm - g * umf(k) * tr0(k,ixnumliq) / ( pifc0(k-1) - pifc0(k) ) ! [ kg/kg/s ] ! nc_im = nc_im - g * umf(k) * tr0(k,ixnumice) / ( pifc0(k-1) - pifc0(k) ) ! [ kg/kg/s ] endif ! 4. Cumulative Penetrative entrainment detrained in the 'kbup' layer ! Explicitly compute the properties detrained penetrative entrained airs in k = kbup layer. if( k .eq. kbup ) then call conden(pmid0(k),thlu_emf(k),qtu_emf(k),thj,qvj,ql_emf_kbup,qi_emf_kbup,qse,id_check) if( id_check .eq. 1 ) then id_exit = .true. if (scverbose) then call write_parallel('------- UW ShCu: exit, conden') end if go to 333 endif if( ql_emf_kbup .gt. 0._r8 ) then ! nl_emf_kbup = tru_emf(k,ixnumliq) else nl_emf_kbup = 0._r8 endif if( qi_emf_kbup .gt. 0._r8 ) then ! ni_emf_kbup = tru_emf(k,ixnumice) else ni_emf_kbup = 0._r8 endif qc_lm = qc_lm - g * emf(k) * ( ql_emf_kbup - ql0(k) ) / ( pifc0(k-1) - pifc0(k) ) ! [ kg/kg/s ] qc_im = qc_im - g * emf(k) * ( qi_emf_kbup - qi0(k) ) / ( pifc0(k-1) - pifc0(k) ) ! [ kg/kg/s ] ! nc_lm = nc_lm - g * emf(k) * ( nl_emf_kbup - tr0(k,ixnumliq) ) / ( pifc0(k-1) - pifc0(k) ) ! [ kg/kg/s ] ! nc_im = nc_im - g * emf(k) * ( ni_emf_kbup - tr0(k,ixnumice) ) / ( pifc0(k-1) - pifc0(k) ) ! [ kg/kg/s ] endif qlten_det(k) = qc_l(k) + qc_lm qiten_det(k) = qc_i(k) + qc_im ! --------------------------------------------------------------------------------- ! ! 'qlten(k)','qiten(k)','qvten(k)','sten(k)' ! ! ! ! --------------------------------------------------------------------------------- ! ! if( use_expconten ) then ! if( use_unicondet ) then ! qc_l(k) = 0._r8 ! qc_i(k) = 0._r8 ! qlten(k) = frc_rasn * dwten(k) + qlten_sink(k) + qlten_det ! qiten(k) = frc_rasn * diten(k) + qiten_sink(k) + qiten_det ! else ! qlten(k) = qc_l(k) + frc_rasn * dwten(k) + ( max( 0._r8, ql0(k) + ( qc_lm + qlten_sink(k) ) * dt ) - ql0(k) ) / dt ! qiten(k) = qc_i(k) + frc_rasn * diten(k) + ( max( 0._r8, qi0(k) + ( qc_im + qiten_sink(k) ) * dt ) - qi0(k) ) / dt ! trten(k,ixnumliq) = max( nc_lm + nlten_sink(k), - tr0(k,ixnumliq) / dt ) ! trten(k,ixnumice) = max( nc_im + niten_sink(k), - tr0(k,ixnumice) / dt ) ! endif ! else ! if( use_unicondet ) then ! qc_l(k) = 0._r8 ! qc_i(k) = 0._r8 ! endif ! qlten(k) = dwten(k) + ( qtten(k) - dwten(k) - diten(k) ) * ( ql0(k) / qt0(k) ) ! qiten(k) = diten(k) + ( qtten(k) - dwten(k) - diten(k) ) * ( qi0(k) / qt0(k) ) ! endif ! ----------------------------------------------------------------- ! ! If combined sinks would result in negative liquid/ice, rescale ! ! such that liquid/ice is zero. ! ! ! ----------------------------------------------------------------- ! if ( ((qc_lm+qlten_sink(k))*dt+ql0(k)).lt.0._r8 ) then totsink = qc_lm+qlten_sink(k) if (totsink.ne.0._r8) then qc_lm = -(ql0(k)/dt) * qc_lm/totsink qlten_sink(k) = -(ql0(k)/dt) * qlten_sink(k)/totsink qlten_det(k) = qc_l(k) + qc_lm end if end if if ( ((qc_im+qiten_sink(k))*dt+qi0(k)).lt.0._r8 ) then totsink = qc_im+qiten_sink(k) if (totsink.ne.0._r8) then qc_im = -(qi0(k)/dt) * qc_im/totsink qiten_sink(k) = -(qi0(k)/dt) * qiten_sink(k)/totsink qiten_det(k) = qc_i(k) + qc_im end if end if qlten(k) = qrten(k) + qlten_sink(k) + qlten_det(k) qiten(k) = qsten(k) + qiten_sink(k) + qiten_det(k) qvten(k) = qtten(k) - qlten(k) - qiten(k) sten(k) = slten(k) + xlv * qlten(k) + xls * qiten(k) qldet(k) = qc_l(k) qidet(k) = qc_i(k) qc(k) = qc_l(k) + qc_i(k) qlten(k) = qlten(k) - qrten(k) qiten(k) = qiten(k) - qsten(k) qtten(k) = qlten(k) + qiten(k) + qvten(k) ! if( ( qv0(k) + qvten(k)*dt ) .lt. qmin(1) .or. & ! ( ql0(k) + qlten(k)*dt ) .lt. qmin(ixcldliq) .or. & ! ( qi0(k) + qiten(k)*dt ) .lt. qmin(ixcldice) ) then ! limit_negcon(i) = 1. ! end if slten(k) = sten(k) - xlv*qlten(k) - xls*qiten(k) slten(k) = slten(k) + xlv * qrten(k) + xls * qsten(k) sten(k) = slten(k) + xlv * qlten(k) + xls * qiten(k) end do ! ---------------------------------------------------------------- ! ! Now treats the 'evaporation' and 'melting' of rain ( qrten ) and ! ! snow ( qsten ) during falling process. ! ! ---------------------------------------------------------------- ! do k = k0, 1, -1 ! 'k' is a layer index : 'k0'('1') is the top ('bottom') layer ! ---------------------------------- ! ! Calculate thermodynamic tendencies ! ! --------------------------------------------------------------------------- ! ! Note that equivalently, we can write tendency formula of 'sten' and 'slten' ! ! by 'sten(k) = sten(k) - xlv*evprain - xls*evpsnow - (xls-xlv)*snowmlt' & ! ! 'slten(k) = sten(k) - xlv*qlten(k) - xls*qiten(k)'. ! ! The above formula is equivalent to the below formula. However below formula ! ! is preferred since we have already imposed explicit constraint on 'ntraprd' ! ! and 'ntsnprd' in case that flxrain(k-1) < 0 & flxsnow(k-1) < 0. ! ! Note : In future, I can elborate the limiting of 'qlten','qvten','qiten' ! ! such that that energy and moisture conservation error is completely ! ! suppressed. ! ! Re-storation to the positive condensate will be performed later below ! ! --------------------------------------------------------------------------- ! end do ! --------------------------------------------------------------- ! ! Prevent the onset-of negative condensate at the next time step ! ! Potentially, this block can be moved just in front of the above ! ! block. ! ! --------------------------------------------------------------- ! ! Modification : I should check whether this 'positive_moisture_single' routine is ! consistent with the one used in UW PBL and cloud macrophysics schemes. ! Modification : Below may overestimate resulting 'ql, qi' if we use the new 'qc_l', 'qc_i' ! in combination with the original computation of qlten, qiten. However, ! if we use new 'qlten,qiten', there is no problem. qv0_star(:k0) = qv0(:k0) + qvten(:k0) * dt ql0_star(:k0) = ql0(:k0) + qlten(:k0) * dt qi0_star(:k0) = qi0(:k0) + qiten(:k0) * dt s0_star(:k0) = s0(:k0) + sten(:k0) * dt call positive_moisture_single( xlv, xls, k0, dt, qmin(1), qmin(ixcldliq), qmin(ixcldice), & dp0, qv0_star, ql0_star, qi0_star, s0_star, qvten, qlten, qiten, sten ) qtten(:k0) = qvten(:k0) + qlten(:k0) + qiten(:k0) slten(:k0) = sten(:k0) - xlv * qlten(:k0) - xls * qiten(:k0) ! --------------------- ! ! Tendencies of tracers ! ! --------------------- ! do m = 1, ncnst ! if( m .ne. ixnumliq .and. m .ne. ixnumice ) then trmin = 0._r8 !qmin(m) !#ifdef MODAL_AERO ! do mm = 1, ntot_amode ! if( m .eq. numptr_amode(mm) ) then ! trmin = 1.e-5_r8 ! goto 55 ! endif ! enddo ! 55 continue !#endif trflx_d(0:k0) = 0._r8 trflx_u(0:k0) = 0._r8 do k = 1, k0-1 ! if( cnst_get_type_byind(m) .eq. 'wet' ) then pdelx = dp0(k) ! else ! pdelx = dpdry0(k) ! endif km1 = k - 1 dum = ( tr0(k,m) - trmin ) * pdelx / g / dt + trflx(km1,m) - trflx(k,m) + trflx_d(km1) trflx_d(k) = min( 0._r8, dum ) enddo do k = k0, 2, -1 ! if( cnst_get_type_byind(m) .eq. 'wet' ) then pdelx = dp0(k) ! else ! pdelx = dpdry0(k) ! endif km1 = k - 1 dum = ( tr0(k,m) - trmin ) * pdelx / g / dt + trflx(km1,m) - trflx(k,m) + & trflx_d(km1) - trflx_d(k) - trflx_u(k) trflx_u(km1) = max( 0._r8, -dum ) enddo do k = 1, k0 ! if( cnst_get_type_byind(m) .eq. 'wet' ) then pdelx = dp0(k) ! else ! pdelx = dpdry0(k) ! endif km1 = k - 1 ! Check : I should re-check whether '_u', '_d' are correctly ordered in ! the below tendency computation. trten(k,m) = ( trflx(km1,m) - trflx(k,m) + & trflx_d(km1) - trflx_d(k) + & trflx_u(km1) - trflx_u(k) ) * g / pdelx enddo ! endif enddo ! ---------------------------------------------------------------- ! ! Cumpute default diagnostic outputs ! ! Note that since 'qtu(krel-1:kpen-1)' & 'thlu(krel-1:kpen-1)' has ! ! been adjusted after detraining cloud condensate into environment ! ! during cumulus updraft motion, below calculations will exactly ! ! reproduce in-cloud properties as shown in the output analysis. ! ! ---------------------------------------------------------------- ! call conden(prel,thlu(krel-1),qtu(krel-1),thj,qvj,qlj,qij,qse,id_check) if( id_check .eq. 1 ) then exit_conden(i) = 1. id_exit = .true. if (scverbose) then call write_parallel('------- UW ShCu: exit, conden') end if go to 333 end if qcubelow = qlj + qij qlubelow = qlj qiubelow = qij rcwp = 0._r8 rlwp = 0._r8 riwp = 0._r8 ! --------------------------------------------------------------------- ! ! In the below calculations, I explicitly considered cloud base ( LCL ) ! ! and cloud top height ( pifc0(kpen-1) + ppen ) ! ! ----------------------------------------------------------------------! do k = krel, kpen ! This is a layer index ! ------------------------------------------------------------------ ! ! Calculate cumulus condensate at the upper interface of each layer. ! ! Note 'ppen < 0' and at 'k=kpen' layer, I used 'thlu_top'&'qtu_top' ! ! which explicitly considered zero or non-zero 'fer(kpen)'. ! ! ------------------------------------------------------------------ ! if( k .eq. kpen ) then call conden(pifc0(k-1)+ppen,thlu_top,qtu_top,thj,qvj,qlj,qij,qse,id_check) else call conden(pifc0(k),thlu(k),qtu(k),thj,qvj,qlj,qij,qse,id_check) endif if( id_check .eq. 1 ) then exit_conden(i) = 1. id_exit = .true. if (scverbose) then call write_parallel('------- UW ShCu: exit, conden') end if go to 333 end if ! ---------------------------------------------------------------- ! ! Calculate in-cloud mean LWC ( qlu(k) ), IWC ( qiu(k) ), & layer ! ! mean cumulus fraction ( cufrc(k) ), vertically-integrated layer ! ! mean LWP and IWP. Expel some of in-cloud condensate at the upper ! ! interface if it is largr than criqc. Note cumulus cloud fraction ! ! is assumed to be twice of core updraft fractional area. Thus LWP ! ! and IWP will be twice of actual value coming from our scheme. ! ! ---------------------------------------------------------------- ! qcu(k) = 0.5_r8 * ( qcubelow + qlj + qij ) qlu(k) = 0.5_r8 * ( qlubelow + qlj ) qiu(k) = 0.5_r8 * ( qiubelow + qij ) cufrc(k) = ( ufrc(k-1) + ufrc(k) ) if( k .eq. krel ) then cufrc(k) = ( ufrclcl + ufrc(k) )*( prel - pifc0(k) )/( pifc0(k-1) - pifc0(k) ) else if( k .eq. kpen ) then cufrc(k) = ( ufrc(k-1) + 0. )*( -ppen ) /( pifc0(k-1) - pifc0(k) ) if( (qlj + qij) .gt. criqc ) then qcu(k) = 0.5_r8 * ( qcubelow + criqc ) qlu(k) = 0.5_r8 * ( qlubelow + criqc * qlj / ( qlj + qij ) ) qiu(k) = 0.5_r8 * ( qiubelow + criqc * qij / ( qlj + qij ) ) endif endif rcwp = rcwp + ( qlu(k) + qiu(k) ) * ( pifc0(k-1) - pifc0(k) ) / g * cufrc(k) rlwp = rlwp + qlu(k) * ( pifc0(k-1) - pifc0(k) ) / g * cufrc(k) riwp = riwp + qiu(k) * ( pifc0(k-1) - pifc0(k) ) / g * cufrc(k) qcubelow = qlj + qij qlubelow = qlj qiubelow = qij end do ! ------------------------------------ ! ! Cloud top and base interface indices ! ! ------------------------------------ ! cnt = real( kpen ) cnb = real( krel - 1 ) ! ------------------------------------------------------------------------- ! ! End of formal calculation. Below blocks are for implicit CIN calculations ! ! with re-initialization and save variables at iter_cin = 1._r8 ! ! ------------------------------------------------------------------------- ! ! --------------------------------------------------------------- ! ! Adjust the original input profiles for implicit CIN calculation ! ! --------------------------------------------------------------- ! if( iter .ne. iter_cin ) then ! ------------------------------------------------------------------- ! ! Save the output from "iter_cin = 1" ! ! These output will be writed-out if "iter_cin = 1" was not performed ! ! for some reasons. ! ! ------------------------------------------------------------------- ! qv0_s(:k0) = qv0(:k0) + qvten(:k0) * dt ql0_s(:k0) = ql0(:k0) + qlten(:k0) * dt qi0_s(:k0) = qi0(:k0) + qiten(:k0) * dt s0_s(:k0) = s0(:k0) + sten(:k0) * dt u0_s(:k0) = u0(:k0) + uten(:k0) * dt v0_s(:k0) = v0(:k0) + vten(:k0) * dt ! qt0_s(:k0) = qv0_s(:k0) + ql0_s(:k0) + qi0_s(:k0) t0_s(:k0) = t0(:k0) + sten(:k0) * dt / cp do m = 1, ncnst tr0_s(:k0,m) = tr0(:k0,m) + trten(:k0,m) * dt enddo umf_s(0:k0) = umf(0:k0) qvten_s(:k0) = qvten(:k0) qlten_s(:k0) = qlten(:k0) qiten_s(:k0) = qiten(:k0) sten_s(:k0) = sten(:k0) uten_s(:k0) = uten(:k0) vten_s(:k0) = vten(:k0) qrten_s(:k0) = qrten(:k0) qsten_s(:k0) = qsten(:k0) cush_s = cush cufrc_s(:k0) = cufrc(:k0) slflx_s(0:k0) = slflx(0:k0) qtflx_s(0:k0) = qtflx(0:k0) qcu_s(:k0) = qcu(:k0) qlu_s(:k0) = qlu(:k0) qiu_s(:k0) = qiu(:k0) fer_s(:k0) = fer(:k0) fdr_s(:k0) = fdr(:k0) xc_s(:k0) = xco(:k0) cin_s = cin cinlcl_s = cinlcl cbmf_s = cbmf qc_s(:k0) = qc(:k0) qldet_s(:k0) = qlten_det(:k0) qidet_s(:k0) = qiten_det(:k0) qlsub_s(:k0) = qlten_sink(:k0) qisub_s(:k0) = qiten_sink(:k0) cnt_s = cnt cnb_s = cnb qtten_s(:k0) = qtten(:k0) slten_s(:k0) = slten(:k0) ufrc_s(0:k0) = ufrc(0:k0) uflx_s(0:k0) = uflx(0:k0) vflx_s(0:k0) = vflx(0:k0) ufrcinvbase_s = ufrcinvbase ufrclcl_s = ufrclcl winvbase_s = winvbase wlcl_s = wlcl plcl_s = plcl pinv_s = pifc0(kinv-1) plfc_s = plfc prel_s = prel pbup_s = pifc0(kbup) ppen_s = pifc0(kpen-1) + ppen qtsrc_s = qtsrc thlsrc_s = thlsrc thvlsrc_s = thvlsrc emfkbup_s = emf(kbup) cbmflimit_s = cbmflimit tkeavg_s = tkeavg zinv_s = zifc0(kinv-1) rcwp_s = rcwp rlwp_s = rlwp riwp_s = riwp wu_s(0:k0) = wu(0:k0) qtu_s(0:k0) = qtu(0:k0) thlu_s(0:k0) = thlu(0:k0) thvu_s(0:k0) = thvu(0:k0) uu_s(0:k0) = uu(0:k0) vu_s(0:k0) = vu(0:k0) qtu_emf_s(0:k0) = qtu_emf(0:k0) thlu_emf_s(0:k0) = thlu_emf(0:k0) uu_emf_s(0:k0) = uu_emf(0:k0) vu_emf_s(0:k0) = vu_emf(0:k0) uemf_s(0:k0) = uemf(0:k0) dwten_s(:k0) = dwten(:k0) diten_s(:k0) = diten(:k0) excessu_arr_s(:k0) = excessu_arr(:k0) excess0_arr_s(:k0) = excess0_arr(:k0) xc_arr_s(:k0) = xc_arr(:k0) aquad_arr_s(:k0) = aquad_arr(:k0) bquad_arr_s(:k0) = bquad_arr(:k0) cquad_arr_s(:k0) = cquad_arr(:k0) bogbot_arr_s(:k0) = bogbot_arr(:k0) bogtop_arr_s(:k0) = bogtop_arr(:k0) do m = 1, ncnst trten_s(:k0,m) = trten(:k0,m) trflx_s(0:k0,m) = trflx(0:k0,m) tru_s(0:k0,m) = tru(0:k0,m) tru_emf_s(0:k0,m) = tru_emf(0:k0,m) enddo ! ----------------------------------------------------------------------------- ! ! Recalculate environmental variables for new cin calculation at "iter_cin = 2" ! ! using the updated state variables. Perform only for variables necessary for ! ! the new cin calculation. ! ! ----------------------------------------------------------------------------- ! qv0(:k0) = qv0_s(:k0) ql0(:k0) = ql0_s(:k0) qi0(:k0) = qi0_s(:k0) s0(:k0) = s0_s(:k0) t0(:k0) = t0_s(:k0) qt0(:k0) = (qv0(:k0) + ql0(:k0) + qi0(:k0)) thl0(:k0) = (t0(:k0) - xlv*ql0(:k0)/cp - xls*qi0(:k0)/cp)/exnmid0(:k0) thvl0(:k0) = (1._r8 + zvir*qt0(:k0))*thl0(:k0) ssthl0 = slope(k0,thl0,pmid0) ! Dimension of ssthl0(:k0) is implicit ssqt0 = slope(k0,qt0 ,pmid0) ssu0 = slope(k0,u0 ,pmid0) ssv0 = slope(k0,v0 ,pmid0) do m = 1, ncnst sstr0(:k0,m) = slope(k0,tr0(:k0,m),pmid0) enddo do k = 1, k0 thl0bot = thl0(k) + ssthl0(k) * ( pifc0(k-1) - pmid0(k) ) qt0bot = qt0(k) + ssqt0(k) * ( pifc0(k-1) - pmid0(k) ) call conden(pifc0(k-1),thl0bot,qt0bot,thj,qvj,qlj,qij,qse,id_check) if( id_check .eq. 1 ) then exit_conden(i) = 1. id_exit = .true. if (scverbose) then call write_parallel('------- UW ShCu: exit, conden') end if go to 333 end if thv0bot(k) = thj * ( 1._r8 + zvir*qvj - qlj - qij ) thvl0bot(k) = thl0bot * ( 1._r8 + zvir*qt0bot ) thl0top = thl0(k) + ssthl0(k) * ( pifc0(k) - pmid0(k) ) qt0top = qt0(k) + ssqt0(k) * ( pifc0(k) - pmid0(k) ) call conden(pifc0(k),thl0top,qt0top,thj,qvj,qlj,qij,qse,id_check) if( id_check .eq. 1 ) then exit_conden(i) = 1._r8 id_exit = .true. if (scverbose) then call write_parallel('------- UW ShCu: exit, conden') end if go to 333 end if thv0top(k) = thj * ( 1._r8 + zvir*qvj - qlj - qij ) thvl0top(k) = thl0top * ( 1._r8 + zvir*qt0top ) end do endif ! End of 'if(iter .ne. iter_cin)' if sentence. end do ! iter loop ! ----------------------- ! ! Update Output Variables ! ! ----------------------- ! cldhgt_out(i) = cldhgt umf_out(i,0:k0) = umf(0:k0) !the indices are not reversed, these variables go into compute_mcshallow_inv qvten_out(i,:k0) = qvten(:k0) qlten_out(i,:k0) = qlten(:k0) qiten_out(i,:k0) = qiten(:k0) sten_out(i,:k0) = sten(:k0) uten_out(i,:k0) = uten(:k0) vten_out(i,:k0) = vten(:k0) qrten_out(i,:k0) = qrten(:k0) qsten_out(i,:k0) = qsten(:k0) cufrc_out(i,:k0) = cufrc(:k0) qcu_out(i,:k0) = qcu(:k0) qlu_out(i,:k0) = qlu(:k0) qiu_out(i,:k0) = qiu(:k0) cush_inout(i) = cush cbmf_out(i) = cbmf qc_out(i,:k0) = qc(:k0) qldet_out(i,:k0) = qlten_det(:k0) qidet_out(i,:k0) = qiten_det(:k0) qlsub_out(i,:k0) = qlten_sink(:k0) qisub_out(i,:k0) = qiten_sink(:k0) ndrop_out(i,:k0) = qlten_det(:k0)/(4188.787_r8*rdrop**3) ! ndrop_out(i,:k0) = qlten_det(:k0)/(4.19e-12_r8) !(1.15e-11_r8) ! /drop mass nice_out(i,:k0) = qiten_det(:k0)/(3.0e-10_r8) ! /crystal mass cnt_out(i) = cnt cnb_out(i) = cnb do m = 1, ncnst trten_out(i,:k0,m) = trten(:k0,m) tr0_inout(i,:k0,m) = tr0_inout(i,:k0,m) + trten(:k0,m) * dt enddo ! ------------------------------------------------- ! ! Below are specific diagnostic output for detailed ! ! analysis of cumulus scheme ! ! ------------------------------------------------- ! fer_out(i,1:k0) = fer(:k0) fdr_out(i,1:k0) = fdr(:k0) xc_out(i,1:k0) = xco(:k0) cinh_out(i) = cin cinlclh_out(i) = cinlcl qtten_out(i,1:k0) = qtten(:k0) slten_out(i,1:k0) = slten(:k0) ufrc_out(i,0:k0) = ufrc(0:k0) ! uflx_out(i,0:k0) = uflx(0:k0) ! vflx_out(i,0:k0) = vflx(0:k0) ufrcinvbase_out(i) = ufrcinvbase ufrclcl_out(i) = ufrclcl winvbase_out(i) = winvbase wlcl_out(i) = wlcl plcl_out(i) = plcl pinv_out(i) = pifc0(kinv-1) plfc_out(i) = plfc prel_out(i) = prel pbup_out(i) = pifc0(kbup) ppen_out(i) = pifc0(kpen-1) + ppen qtsrc_out(i) = qtsrc thlsrc_out(i) = thlsrc thvlsrc_out(i) = thvlsrc emfkbup_out(i) = emf(kbup) cbmflimit_out(i) = cbmflimit tkeavg_out(i) = tkeavg zinv_out(i) = zifc0(kinv-1) rcwp_out(i) = rcwp rlwp_out(i) = rlwp riwp_out(i) = riwp wu_out(i,0:k0) = wu(0:k0) qtu_out(i,0:k0) = qtu(0:k0) thlu_out(i,0:k0) = thlu(0:k0) thvu_out(i,0:k0) = thvu(0:k0) uu_out(i,0:k0) = uu(0:k0) vu_out(i,0:k0) = vu(0:k0) qtu_emf_out(i,0:k0) = qtu_emf(0:k0) thlu_emf_out(i,0:k0) = thlu_emf(0:k0) uu_emf_out(i,0:k0) = uu_emf(0:k0) vu_emf_out(i,0:k0) = vu_emf(0:k0) uemf_out(i,0:k0) = uemf(0:k0) dwten_out(i,1:k0) = dwten(:k0) diten_out(i,1:k0) = diten(:k0) excessu_arr_out(i,1:k0) = excessu_arr(:k0) excess0_arr_out(i,1:k0) = excess0_arr(:k0) xc_arr_out(i,1:k0) = xc_arr(:k0) aquad_arr_out(i,1:k0) = aquad_arr(:k0) bquad_arr_out(i,1:k0) = bquad_arr(:k0) cquad_arr_out(i,1:k0) = cquad_arr(:k0) bogbot_arr_out(i,1:k0) = bogbot_arr(:k0) bogtop_arr_out(i,1:k0) = bogtop_arr(:k0) do m = 1, ncnst trflx_out(i,0:k0,m) = trflx(0:k0,m) tru_out(i,0:k0,m) = tru(0:k0,m) tru_emf_out(i,0:k0,m) = tru_emf(0:k0,m) enddo 333 if (id_exit) then exit_uwcu(i) = 1. if (scverbose) then call write_parallel('------- UW ShCu: Exited!') end if ! --------------------------------------------------------------------- ! ! Initialize output variables when cumulus convection was not performed.! ! --------------------------------------------------------------------- ! umf_out(i,0:k0) = 0._r8 qvten_out(i,:k0) = 0._r8 qlten_out(i,:k0) = 0._r8 qiten_out(i,:k0) = 0._r8 sten_out(i,:k0) = 0._r8 uten_out(i,:k0) = 0._r8 vten_out(i,:k0) = 0._r8 qrten_out(i,:k0) = 0._r8 qsten_out(i,:k0) = 0._r8 cufrc_out(i,:k0) = 0._r8 qcu_out(i,:k0) = 0._r8 qlu_out(i,:k0) = 0._r8 qiu_out(i,:k0) = 0._r8 cush_inout(i) = -1._r8 cbmf_out(i) = 0._r8 qc_out(i,:k0) = 0._r8 qldet_out(i,:k0) = 0._r8 qidet_out(i,:k0) = 0._r8 cnt_out(i) = 1._r8 cnb_out(i) = real(k0) fer_out(i,1:k0) = MAPL_UNDEF fdr_out(i,1:k0) = MAPL_UNDEF xc_out(i,1:k0) = MAPL_UNDEF cinh_out(i) = cin cinlclh_out(i) = cinlcl qtten_out(i,k0:1:-1) = 0._r8 slten_out(i,k0:1:-1) = 0._r8 ufrc_out(i,k0:0:-1) = 0._r8 ! uflx_out(i,k0:0:-1) = 0._r8 ! vflx_out(i,k0:0:-1) = 0._r8 ufrcinvbase_out(i) = 0._r8 ufrclcl_out(i) = 0._r8 winvbase_out(i) = 0._r8 wlcl_out(i) = MAPL_UNDEF plcl_out(i) = MAPL_UNDEF pinv_out(i) = MAPL_UNDEF prel_out(i) = MAPL_UNDEF plfc_out(i) = MAPL_UNDEF pbup_out(i) = MAPL_UNDEF ppen_out(i) = MAPL_UNDEF qtsrc_out(i) = MAPL_UNDEF thlsrc_out(i) = MAPL_UNDEF thvlsrc_out(i) = MAPL_UNDEF emfkbup_out(i) = 0._r8 cbmflimit_out(i) = 0._r8 tkeavg_out(i) = tkeavg zinv_out(i) = 0._r8 rcwp_out(i) = 0._r8 rlwp_out(i) = 0._r8 riwp_out(i) = 0._r8 wu_out(i,k0:0:-1) = MAPL_UNDEF qtu_out(i,k0:0:-1) = MAPL_UNDEF thlu_out(i,k0:0:-1) = MAPL_UNDEF thvu_out(i,k0:0:-1) = MAPL_UNDEF uu_out(i,k0:0:-1) = MAPL_UNDEF vu_out(i,k0:0:-1) = MAPL_UNDEF qtu_emf_out(i,k0:0:-1) = MAPL_UNDEF thlu_emf_out(i,k0:0:-1) = MAPL_UNDEF uu_emf_out(i,k0:0:-1) = MAPL_UNDEF vu_emf_out(i,k0:0:-1) = MAPL_UNDEF uemf_out(i,k0:0:-1) = MAPL_UNDEF dwten_out(i,k0:1:-1) = 0._r8 diten_out(i,k0:1:-1) = 0._r8 excessu_arr_out(i,k0:1:-1) = 0._r8 excess0_arr_out(i,k0:1:-1) = 0._r8 xc_arr_out(i,k0:1:-1) = 0._r8 aquad_arr_out(i,k0:1:-1) = 0._r8 bquad_arr_out(i,k0:1:-1) = 0._r8 cquad_arr_out(i,k0:1:-1) = 0._r8 bogbot_arr_out(i,k0:1:-1) = 0._r8 bogtop_arr_out(i,k0:1:-1) = 0._r8 do m = 1, ncnst trten_out(i,:k0,m) = 0._r8 trflx_out(i,k0:0:-1,m) = 0._r8 tru_out(i,k0:0:-1,m) = 0._r8 tru_emf_out(i,k0:0:-1,m) = 0._r8 enddo end if end do ! column i loop return end subroutine compute_uwshcu function slope(k0,field,p0) integer, intent(in):: k0 real*8 :: slope(k0) real*8, intent(in) :: field(k0) real*8, intent(in) :: p0(k0) real*8 :: below real*8 :: above integer :: k below = ( field(2) - field(1) ) / ( p0(2) - p0(1) ) do k = 2, k0 above = ( field(k) - field(k-1) ) / ( p0(k) - p0(k-1) ) if ( above .gt. 0._r8 ) then slope(k-1) = max(0._r8,min(above,below)) else slope(k-1) = min(0._r8,max(above,below)) end if below = above end do slope(k0) = slope(k0-1) return end function slope function qsinvert(qt,thl,ps_in) ! ----------------------------------------------------------------- ! ! Function calculating saturation pressure ps (or pLCL) from qt and ! ! thl ( liquid potential temperature, NOT liquid virtual potential ! ! temperature) by inverting Bolton formula. I should check later if ! ! current use of 'leff' instead of 'xlv' here is reasonable or not. ! ! ----------------------------------------------------------------- ! implicit none real*8 :: ps_in real*8 :: qsinvert real*8,intent(in) :: qt, thl real*8 :: Pis, err, dlnqsdT, dTdPis, dqsdT real*8 :: dPisdps, dlnqsdps, derrdps, dps real*8 :: rhi, TLCL, PiLCL, psmin, dpsmax real*8 :: Ti,qs,Ts,ps real :: Tgeos,Qgeos,Pgeos integer :: i real*8 :: es ! saturation vapor pressure real*8 :: gam ! (L/cp)*dqs/dT real*8 :: leff, nu psmin = 10000._r8 ! Default saturation pressure [Pa] if iteration does not converge dpsmax = 1._r8 ! Tolerance [Pa] for convergence of iteration ! ------------------------------------ ! ! Calculate best initial guess of pLCL ! ! ------------------------------------ ! Ti = thl*(ps_in/p00)**rovcp Tgeos = Ti Pgeos = ps_in qs = dble( GEOS_QSAT(Tgeos,Pgeos/100.) ) es = ps_in * qs / ( ep2 + (1._r8-ep2)*qs ) rhi = qt/qs if( rhi .le. 0.01_r8 ) then ! if (scverbose) then call write_parallel('Source air is too dry and pLCL is set to psmin in uwshcu.F90') ! end if qsinvert = psmin return end if ! print *,'Ti,Rhi,thl=',Ti,rhi,thl TLCL = 55._r8 + 1._r8/(1._r8/(Ti-55._r8)-log(rhi)/2840._r8) ! Bolton's formula. MWR.1980.Eq.(22) PiLCL = TLCL/thl ps = p00*(PiLCL)**(1._r8/rovcp) do i = 1, 10 Pis = (ps/p00)**rovcp ! Exner function Ts = thl*Pis Tgeos = Ts Pgeos = ps qs = dble( GEOS_QSAT(Tgeos,Pgeos/100.) ) Qgeos = qs dqsdT = dble( GEOS_DQSAT(Tgeos, Pgeos/100., QSAT=Qgeos ) ) gam = (xlv/cp)*dqsdT err = qt - qs nu = max(min((268._r8 - Ts)/20._r8,1.0_r8),0.0_r8) leff = (1._r8 - nu)*xlv + nu*xls dlnqsdT = gam*(cp/leff)/qs dTdPis = thl dPisdps = rovcp*Pis/ps dlnqsdps = -1._r8/(ps - (1. - ep2)*es) derrdps = -qs*(dlnqsdT * dTdPis * dPisdps + dlnqsdps) if (derrdps.eq.0._r8) then print *,"QSINVERT: derrdps=0 !!!" end if dps = -err/derrdps ps = ps + dps if( ps .lt. 0._r8 ) then ! if (scverbose) then ! call write_parallel('------- UW ShCu: pLCL iteration is negative and set to psmin') ! end if qsinvert = psmin return end if if( abs(dps) .le. dpsmax ) then qsinvert = ps return end if end do ! if (scverbose) then ! call write_parallel('------- UW ShCu: pLCL does not converge and is set to psmin') ! end if qsinvert = psmin return end function qsinvert real*8 function compute_alpha(del_CIN,ke) ! ------------------------------------------------ ! ! Subroutine to compute proportionality factor for ! ! implicit CIN calculation. ! ! ------------------------------------------------ ! real*8 :: del_CIN, ke real*8 :: x0, x1 integer :: iteration x0 = 0._r8 do iteration = 1, 10 x1 = x0 - (exp(-x0*ke*del_CIN) - x0)/(-ke*del_CIN*exp(-x0*ke*del_CIN) - 1.) x0 = x1 end do compute_alpha = x0 return end function compute_alpha real*8 function compute_mumin2(mulcl,rmaxfrac,mulow) ! --------------------------------------------------------- ! ! Subroutine to compute critical 'mu' (normalized CIN) such ! ! that updraft fraction at the LCL is equal to 'rmaxfrac'. ! ! --------------------------------------------------------- ! real*8 :: mulcl, rmaxfrac, mulow real*8 :: x0, x1, ex, ef, exf, f, fs integer :: iteration x0 = mulow do iteration = 1, 10 ex = exp(-x0**2) ef = erfc(x0) exf = ex/ef f = 0.5_r8*exf**2 - 0.5_r8*(ex/2._r8/rmaxfrac)**2 - (mulcl*2.5066_r8/2._r8)**2 fs = (2._r8*exf**2)*(exf/sqrt(MAPL_PI)-x0) + (0.5_r8*x0*ex**2)/(rmaxfrac**2) x1 = x0 - f/fs x0 = x1 end do compute_mumin2 = x0 return end function compute_mumin2 real*8 function compute_ppen(wtwb,D,bogbot,bogtop,rho0j,dpen) ! ----------------------------------------------------------- ! ! Subroutine to compute critical 'ppen[Pa]<0' ( pressure dis. ! ! from 'pifc0(kpen-1)' to the cumulus top where cumulus updraft ! ! vertical velocity is exactly zero ) by considering exact ! ! non-zero fer(kpen). ! ! ----------------------------------------------------------- ! real*8 :: wtwb, D, bogbot, bogtop, rho0j, dpen real*8 :: x0, x1, f, fs, SB, s00, aux integer :: iteration ! Buoyancy slope SB = ( bogtop - bogbot ) / dpen ! Sign of slope, 'f' at x = 0 ! If 's00>0', 'w' increases with height. s00 = bogbot / rho0j - D * wtwb if( D*dpen .lt. 1.e-4_r8 ) then if( s00 .ge. 0._r8 ) then x0 = dpen else x0 = max(0._r8,min(dpen,-0.5_r8*wtwb/s00)) endif else if( s00 .ge. 0._r8 ) then x0 = dpen else x0 = 0._r8 endif do iteration = 1, 5 aux = min(max(-2._r8*D*x0, -20.0), 20.0) f = exp(aux)*(wtwb-(bogbot-SB/(2._r8*D))/(D*rho0j)) + & (SB*x0+bogbot-SB/(2._r8*D))/(D*rho0j) fs = -2._r8*D*exp(aux)*(wtwb-(bogbot-SB/(2._r8*D))/(D*rho0j)) + & (SB)/(D*rho0j) ! if( fs .ge. 0._r8 ) then ! fs = max(fs, 1.e-8_r8) ! else ! fs = min(fs,-1.e-8_r8) ! endif x1 = x0 - f/fs x0 = x1 ! if (x0.gt.dpen*10._r8) then ! print *,'compute_ppen: x0 limit reached' ! x0 = 0._r8 ! exit ! end if end do endif compute_ppen = -max(0._r8,min(dpen,x0)) end function compute_ppen subroutine fluxbelowinv(cbmf,ps0,mkx,kinv,dt,xsrc,xmean,xtopin,xbotin,xflx) ! ------------------------------------------------------------------------- ! ! Subroutine to calculate turbulent fluxes at and below 'kinv-1' interfaces.! ! Check in the main program such that input 'cbmf' should not be zero. ! ! If the reconstructed inversion height does not go down below the 'kinv-1' ! ! interface, then turbulent flux at 'kinv-1' interface is simply a product ! ! of 'cmbf' and 'qtsrc-xbot' where 'xbot' is the value at the top interface ! ! of 'kinv-1' layer. This flux is linearly interpolated down to the surface ! ! assuming turbulent fluxes at surface are zero. If reconstructed inversion ! ! height goes down below the 'kinv-1' interface, subsidence warming &drying ! ! measured by 'xtop-xbot', where 'xtop' is the value at the base interface ! ! of 'kinv+1' layer, is added ONLY to the 'kinv-1' layer, using appropriate ! ! mass weighting ( rpinv and rcbmf, or rr = rpinv / rcbmf ) between current ! ! and next provisional time step. Also impose a limiter to enforce outliers ! ! of thermodynamic variables in 'kinv' layer to come back to normal values ! ! at the next step. ! ! ------------------------------------------------------------------------- ! integer, intent(in) :: mkx, kinv real*8, intent(in) :: cbmf, xsrc, xmean, xtopin, xbotin real*8, intent(in), dimension(0:mkx) :: ps0 real*8, intent(out), dimension(0:mkx) :: xflx integer k real, intent(in) :: dt real*8 rcbmf, rpeff, dp, rr, pinv_eff, xtop, xbot, pinv, xtop_ori, xbot_ori xflx(0:mkx) = 0._r8 dp = ps0(kinv-1) - ps0(kinv) xbot = xbotin xtop = xtopin ! -------------------------------------- ! ! Compute reconstructed inversion height ! ! -------------------------------------- ! xtop_ori = xtop xbot_ori = xbot rcbmf = ( cbmf * g * dt ) / dp ! Can be larger than 1 : 'OK' if( xbot .ge. xtop ) then rpeff = ( xmean - xtop ) / max( 1.e-20_r8, xbot - xtop ) else rpeff = ( xmean - xtop ) / min( -1.e-20_r8, xbot - xtop ) endif rpeff = min( max(0._r8,rpeff), 1. ) ! As of this, 0<= rpeff <= 1 if( rpeff .eq. 0._r8 .or. rpeff .eq. 1._r8 ) then xbot = xmean xtop = xmean endif ! Below two commented-out lines are the old code replacing the above 'if' block. ! if(rpeff.eq.1) xbot = xmean ! if(rpeff.eq.0) xtop = xmean rr = rpeff / rcbmf pinv = ps0(kinv-1) - rpeff * dp ! "pinv" before detraining mass pinv_eff = ps0(kinv-1) + ( rcbmf - rpeff ) * dp ! Effective "pinv" after detraining mass ! ----------------------------------------------------------------------- ! ! Compute turbulent fluxes. ! ! Below two cases exactly converges at 'kinv-1' interface when rr = 1. ! ! ----------------------------------------------------------------------- ! do k = 0, kinv - 1 xflx(k) = cbmf * ( xsrc - xbot ) * ( ps0(0) - ps0(k) ) / ( ps0(0) - pinv ) end do if( rr .le. 1. ) then xflx(kinv-1) = xflx(kinv-1) - ( 1. - rr ) * cbmf * ( xtop_ori - xbot_ori ) endif return end subroutine fluxbelowinv subroutine positive_moisture_single( xlv, xls, mkx, dt, qvmin, qlmin, qimin, dp, qv, ql, qi, s, qvten, qlten, qiten, sten ) ! ------------------------------------------------------------------------------- ! ! If any 'ql < qlmin, qi < qimin, qv < qvmin' are developed in any layer, ! ! force them to be larger than minimum value by (1) condensating water vapor ! ! into liquid or ice, and (2) by transporting water vapor from the very lower ! ! layer. '2._r8' is multiplied to the minimum values for safety. ! ! Update final state variables and tendencies associated with this correction. ! ! If any condensation happens, update (s,t) too. ! ! Note that (qv,ql,qi,s) are final state variables after applying corresponding ! ! input tendencies and corrective tendencies ! ! ------------------------------------------------------------------------------- ! implicit none integer, intent(in) :: mkx real*8, intent(in) :: xlv, xls real*8, intent(in) :: qvmin, qlmin, qimin real*8, intent(in) :: dp(mkx) real, intent(in) :: dt real*8, intent(inout) :: qv(mkx), ql(mkx), qi(mkx), s(mkx) real*8, intent(inout) :: qvten(mkx), qlten(mkx), qiten(mkx), sten(mkx) integer k real*8 dql, dqi, dqv, sum, aa, dum do k = mkx, 1, -1 ! From the top to the 1st (lowest) layer from the surface dql = max(0._r8,1._r8*qlmin-ql(k)) dqi = max(0._r8,1._r8*qimin-qi(k)) qlten(k) = qlten(k) + dql/dt qiten(k) = qiten(k) + dqi/dt qvten(k) = qvten(k) - (dql+dqi)/dt sten(k) = sten(k) + xlv * (dql/dt) + xls * (dqi/dt) ql(k) = ql(k) + dql qi(k) = qi(k) + dqi qv(k) = qv(k) - dql - dqi s(k) = s(k) + xlv * dql + xls * dqi dqv = max(0.,1.*qvmin-qv(k)) qvten(k) = qvten(k) + dqv/dt qv(k) = qv(k) + dqv if( k .ne. 1 ) then qv(k-1) = qv(k-1) - dqv*dp(k)/dp(k-1) qvten(k-1) = qvten(k-1) - dqv*dp(k)/dp(k-1)/dt endif qv(k) = max(qv(k),qvmin) ql(k) = max(ql(k),qlmin) qi(k) = max(qi(k),qimin) end do ! Extra moisture used to satisfy 'qv(i,1)=qvmin' is proportionally ! extracted from all the layers that has 'qv > 2*qvmin'. This fully ! preserves column moisture. if( dqv .gt. 1.e-20_r8 ) then sum = 0. do k = 1, mkx if( qv(k) .gt. 2._r8*qvmin ) sum = sum + qv(k)*dp(k) enddo aa = dqv*dp(1)/max(1.e-20_r8,sum) if( aa .lt. 0.5_r8 ) then do k = 1, mkx if( qv(k) .gt. 2._r8*qvmin ) then dum = aa*qv(k) qv(k) = qv(k) - dum qvten(k) = qvten(k) - dum/dt endif enddo else ! if (scverbose) then ! call write_parallel( 'Full positive_moisture is impossible in uwshcu') ! end if endif endif return end subroutine positive_moisture_single subroutine conden(p,thl,qt,th,qv,ql,qi,rvls,id_check) implicit none real*8, intent(in) :: p real*8, intent(in) :: thl real*8, intent(in) :: qt real*8, intent(out) :: th real*8, intent(out) :: qv real*8, intent(out) :: ql real*8, intent(out) :: qi real*8, intent(out) :: rvls integer, intent(out) :: id_check ! Local real*8 :: tc real :: temps,ps real*8 :: leff,nu,qc integer :: iteration real*8 :: qs ! Saturation specific humidity tc = thl*exnerfn(p) nu = max(min((268._r8-tc)/20._r8,1.0_r8),0.0_r8) ! ice fraction of condensate leff = (1._r8- nu)*xlv + nu*xls ! effective latent heat temps = tc ps = p qs = GEOS_QSAT(temps, ps/100.) rvls = qs if ( qs .ge. qt ) then ! no condensation id_check = 0 qv = qt qc = 0. ql = 0. qi = 0. th = thl !tc/exnerfn(p) else ! condensation do iteration = 1,10 temps = temps + ( (tc-temps)*cp/leff + qt - rvls )/( cp/leff + ep2*leff*rvls/(r*temps*temps) ) qs = GEOS_QSAT(temps, ps/100.) rvls = qs end do qc = max( qt-qs, 0._r8 ) qv = qt - qc ql = qc*(1._r8 - nu) qi = nu*qc th = temps/exnerfn(p) if ( abs((temps-(leff/cp)*qc)-tc) .ge. 1._r8 ) then id_check = 1 else id_check = 0 end if end if return end subroutine conden subroutine roots(a,b,c,r1,r2,status) ! --------------------------------------------------------- ! ! Subroutine to solve the second order polynomial equation. ! ! I should check this subroutine later. ! ! --------------------------------------------------------- ! real*8, intent(in) :: a real*8, intent(in) :: b real*8, intent(in) :: c real*8, intent(out) :: r1 real*8, intent(out) :: r2 integer , intent(out) :: status real*8 :: q status = 0 if( a .eq. 0._r8 ) then ! Form b*x + c = 0 if( b .eq. 0._r8 ) then ! Failure: c = 0 status = 1 else ! b*x + c = 0 r1 = -c/b endif r2 = r1 else if( b .eq. 0._r8 ) then ! Form a*x**2 + c = 0 if( a*c .gt. 0._r8 ) then ! Failure: x**2 = -c/a < 0 status = 2 else ! x**2 = -c/a r1 = sqrt(-c/a) endif r2 = -r1 else ! Form a*x**2 + b*x + c = 0 if( (b**2 - 4._r8*a*c) .lt. 0._r8 ) then ! Failure, no real roots status = 3 else q = -0.5_r8*(b + sign(1.0_r8,b)*sqrt(b**2 - 4._r8*a*c)) r1 = q/a r2 = c/q ! r1 = -0.5_r8*(b + sign(1.0_r8,b)*sqrt(b**2 - 4._r8*a*c))/a ! r2 = -0.5_r8*(b - sign(1.0_r8,b)*sqrt(b**2 - 4._r8*a*c))/a endif endif endif return end subroutine roots subroutine getbuoy(pbot,thv0bot,ptop,thv0top,thvubot,thvutop,plfc,cin) ! ----------------------------------------------------------- ! ! Subroutine to calculate integrated CIN [ J/kg = m2/s2 ] and ! ! 'cinlcl, plfc' if any. Assume 'thv' is linear in each layer ! ! both for cumulus and environment. Note that this subroutine ! ! only includes positive CIN in calculation - if there is any ! ! negative CIN, it is assumed to be zero. This is slightly ! ! different from 'single_cin' below, where both positive and ! ! negative CIN are included. ! ! ----------------------------------------------------------- ! real*8 pbot,thv0bot,ptop,thv0top,thvubot,thvutop,plfc,cin,frc if( thvubot .gt. thv0bot .and. thvutop .gt. thv0top ) then plfc = pbot return elseif( thvubot .le. thv0bot .and. thvutop .le. thv0top ) then cin = cin - ( (thvubot/thv0bot - 1._r8) + (thvutop/thv0top - 1._r8)) * (pbot - ptop) / & ( pbot/(r*thv0bot*exnerfn(pbot)) + ptop/(r*thv0top*exnerfn(ptop)) ) elseif( thvubot .gt. thv0bot .and. thvutop .le. thv0top ) then frc = ( thvutop/thv0top - 1._r8) / ( (thvutop/thv0top - 1._r8) - (thvubot/thv0bot - 1._r8) ) cin = cin - ( thvutop/thv0top - 1._r8) * ( (ptop + frc*(pbot - ptop)) - ptop ) / & ( pbot/(r*thv0bot*exnerfn(pbot)) + ptop/(r*thv0top*exnerfn(ptop)) ) else frc = ( thvubot/thv0bot - 1._r8) / ( (thvubot/thv0bot - 1._r8) - (thvutop/thv0top - 1._r8) ) plfc = pbot - frc * ( pbot - ptop ) cin = cin - ( thvubot/thv0bot - 1._r8)*(pbot - plfc)/ & ( pbot/(r*thv0bot*exnerfn(pbot)) + ptop/(r*thv0top * exnerfn(ptop))) endif return end subroutine getbuoy function single_cin(pbot,thv0bot,ptop,thv0top,thvubot,thvutop) ! ------------------------------------------------------- ! ! Function to calculate a single layer CIN by summing all ! ! positive and negative CIN. ! ! ------------------------------------------------------- ! real*8 :: single_cin real*8 :: pbot,thv0bot,ptop,thv0top,thvubot,thvutop single_cin = ( (1._r8- thvubot/thv0bot) + (1._r8- thvutop/thv0top)) * ( pbot - ptop ) / & ( pbot/(r*thv0bot*exnerfn(pbot)) + ptop/(r*thv0top*exnerfn(ptop)) ) return end function single_cin end module shallow