module module_cu_camzm use shr_kind_mod, only: r8 => shr_kind_r8 use module_cam_support, only: masterproc, pcols, pver, pverp use cldwat, only: cldwat_fice use physconst, only: cpair, epsilo, gravit, latice, latvap, tmelt, rair, & cpwv, cpliq, rh2o use module_cam_support,only: endrun, iulog implicit none save private public zmconv_readnl public zm_convi public zm_convr public zm_conv_evap public convtran public momtran real(r8), parameter :: unset_r8 = huge(1.0_r8) real(r8) :: zmconv_c0_lnd = unset_r8 real(r8) :: zmconv_c0_ocn = unset_r8 real(r8) :: zmconv_ke = unset_r8 real(r8) rl real(r8) cpres real(r8), parameter :: capelmt = 70._r8 real(r8) :: ke real(r8) :: c0_lnd real(r8) :: c0_ocn real(r8) tau real(r8),parameter :: a = 21.656_r8 real(r8),parameter :: b = 5418._r8 real(r8),parameter :: c1 = 6.112_r8 real(r8),parameter :: c2 = 17.67_r8 real(r8),parameter :: c3 = 243.5_r8 real(r8) :: tfreez real(r8) :: eps1 logical :: no_deep_pbl real(r8) :: rgrav real(r8) :: rgas real(r8) :: grav real(r8) :: cp integer limcnv real(r8),parameter :: tiedke_add = 0.5_r8 contains subroutine zmconv_readnl(nlfile) character(len=*), intent(in) :: nlfile c0_lnd = 0.0059D0 c0_ocn = 0.0450D0 ke = 1.0E-6 end subroutine zmconv_readnl subroutine zm_convi(DT,DX,limcnv_in, no_deep_pbl_in) integer, intent(in) :: limcnv_in logical, intent(in), optional :: no_deep_pbl_in REAL , INTENT(IN) :: DT, DX real(r8) :: DTT character(len=32) :: hgrid real(r8) :: deltax real(r8) :: ref_dx real(r8) :: taumin real(r8) :: taumax limcnv = limcnv_in tfreez = tmelt eps1 = epsilo rl = latvap cpres = cpair rgrav = 1.0_r8/gravit rgas = rair grav = gravit cp = cpres ref_dx = 275000._r8 taumin = 600._r8 taumax = 3600._r8 deltax=dx if ( present(no_deep_pbl_in) ) then no_deep_pbl = no_deep_pbl_in else no_deep_pbl = .true. endif DTT=DT tau = max(taumin, taumax*min(1._r8,deltax/ref_dx)) if ( masterproc ) then write(iulog,*) 'delta X =',deltax call wrf_debug(1,iulog) write(iulog,*) 'Convective relaxation time scale (tau) is a tunable parameter in CAM and is a function of spatial resolution.' write(iulog,*) 'Users are encouraged to consult with the PNNL WRF-CAM5 development team if they want to change tau.' call wrf_debug(1,iulog) write(iulog,*) 'tuning parameters zm_convi: tau',tau call wrf_message(iulog) write(iulog,*) 'Standard 2-deg CAM5.1 sets tau=3600s and is reduced for the mesoscale model WRF' call wrf_debug(1,iulog) write(iulog,*) 'tuning parameters zm_convi: c0_lnd',c0_lnd, ', c0_ocn', c0_ocn call wrf_debug(1,iulog) write(iulog,*) 'tuning parameters zm_convi: ke',ke call wrf_debug(1,iulog) write(iulog,*) 'tuning parameters zm_convi: no_deep_pbl',no_deep_pbl call wrf_debug(1,iulog) endif if (masterproc) write(iulog,*)'**** ZM: DILUTE Buoyancy Calculation ****' call wrf_debug(1,iulog) end subroutine zm_convi subroutine zm_convr(lchnk ,ncol , & t ,qh ,prec ,jctop ,jcbot , & pblh ,zm ,geos ,zi ,qtnd , & heat ,pap ,paph ,dpp , & delt ,mcon ,cme ,cape , & tpert ,dlf ,pflx ,zdu ,rprd , & mu ,md ,du ,eu ,ed , & dp ,dsubcld ,jt ,maxg ,ideep , & lengath ,ql ,rliq ,landfrac) use module_cam_support, only: pcnst =>pcnst_runtime integer, intent(in) :: lchnk integer, intent(in) :: ncol real(r8), intent(in) :: t(pcols,pver) real(r8), intent(in) :: qh(pcols,pver,pcnst) real(r8), intent(in) :: pap(pcols,pver) real(r8), intent(in) :: paph(pcols,pver+1) real(r8), intent(in) :: dpp(pcols,pver) real(r8), intent(in) :: zm(pcols,pver) real(r8), intent(in) :: geos(pcols) real(r8), intent(in) :: zi(pcols,pver+1) real(r8), intent(in) :: pblh(pcols) real(r8), intent(in) :: tpert(pcols) real(r8), intent(in) :: landfrac(pcols) real(r8), intent(out) :: qtnd(pcols,pver) real(r8), intent(out) :: heat(pcols,pver) real(r8), intent(out) :: mcon(pcols,pverp) real(r8), intent(out) :: dlf(pcols,pver) real(r8), intent(out) :: pflx(pcols,pverp) real(r8), intent(out) :: cme(pcols,pver) real(r8), intent(out) :: cape(pcols) real(r8), intent(out) :: zdu(pcols,pver) real(r8), intent(out) :: rprd(pcols,pver) real(r8), intent(out) :: mu(pcols,pver) real(r8), intent(out) :: eu(pcols,pver) real(r8), intent(out) :: du(pcols,pver) real(r8), intent(out) :: md(pcols,pver) real(r8), intent(out) :: ed(pcols,pver) real(r8), intent(out) :: dp(pcols,pver) real(r8), intent(out) :: dsubcld(pcols) real(r8), intent(out) :: jctop(pcols) real(r8), intent(out) :: jcbot(pcols) real(r8), intent(out) :: prec(pcols) real(r8), intent(out) :: rliq(pcols) real(r8) zs(pcols) real(r8) dlg(pcols,pver) real(r8) pflxg(pcols,pverp) real(r8) cug(pcols,pver) real(r8) evpg(pcols,pver) real(r8) mumax(pcols) integer jt(pcols) integer maxg(pcols) integer ideep(pcols) integer lengath real(r8) ql(pcols,pver) real(r8) pblt(pcols) real(r8) q(pcols,pver) real(r8) p(pcols,pver) real(r8) z(pcols,pver) real(r8) s(pcols,pver) real(r8) tp(pcols,pver) real(r8) zf(pcols,pver+1) real(r8) pf(pcols,pver+1) real(r8) qstp(pcols,pver) real(r8) tl(pcols) integer lcl(pcols) integer lel(pcols) integer lon(pcols) integer maxi(pcols) integer index(pcols) real(r8) precip real(r8) qg(pcols,pver) real(r8) tg(pcols,pver) real(r8) pg(pcols,pver) real(r8) zg(pcols,pver) real(r8) sg(pcols,pver) real(r8) tpg(pcols,pver) real(r8) zfg(pcols,pver+1) real(r8) qstpg(pcols,pver) real(r8) ug(pcols,pver) real(r8) vg(pcols,pver) real(r8) cmeg(pcols,pver) real(r8) rprdg(pcols,pver) real(r8) capeg(pcols) real(r8) tlg(pcols) real(r8) landfracg(pcols) integer lclg(pcols) integer lelg(pcols) real(r8) dqdt(pcols,pver) real(r8) dsdt(pcols,pver) real(r8) sd(pcols,pver) real(r8) qd(pcols,pver) real(r8) mc(pcols,pver) real(r8) qhat(pcols,pver) real(r8) qu(pcols,pver) real(r8) su(pcols,pver) real(r8) qs(pcols,pver) real(r8) shat(pcols,pver) real(r8) hmn(pcols,pver) real(r8) hsat(pcols,pver) real(r8) qlg(pcols,pver) real(r8) dudt(pcols,pver) real(r8) dvdt(pcols,pver) real(r8) mb(pcols) integer jlcl(pcols) integer j0(pcols) integer jd(pcols) real(r8) delt integer i integer ii integer k integer msg real(r8) qdifr real(r8) sdifr msg = limcnv - 1 qtnd(:,:) = 0._r8 heat(:,:) = 0._r8 mcon(:,:) = 0._r8 rliq(:ncol) = 0._r8 prec(:ncol) = 0._r8 do k = 1,pver do i = 1,ncol dqdt(i,k) = 0._r8 dsdt(i,k) = 0._r8 dudt(i,k) = 0._r8 dvdt(i,k) = 0._r8 pflx(i,k) = 0._r8 pflxg(i,k) = 0._r8 cme(i,k) = 0._r8 rprd(i,k) = 0._r8 zdu(i,k) = 0._r8 ql(i,k) = 0._r8 qlg(i,k) = 0._r8 dlf(i,k) = 0._r8 dlg(i,k) = 0._r8 end do end do do i = 1,ncol pflx(i,pverp) = 0 pflxg(i,pverp) = 0 end do do i = 1,ncol pblt(i) = pver dsubcld(i) = 0._r8 jctop(i) = pver jcbot(i) = 1 end do do i = 1,ncol zs(i) = geos(i)*rgrav pf(i,pver+1) = paph(i,pver+1)*0.01_r8 zf(i,pver+1) = zi(i,pver+1) + zs(i) end do do k = 1,pver do i = 1,ncol p(i,k) = pap(i,k)*0.01_r8 pf(i,k) = paph(i,k)*0.01_r8 z(i,k) = zm(i,k) + zs(i) zf(i,k) = zi(i,k) + zs(i) end do end do do k = pver - 1,msg + 1,-1 do i = 1,ncol if (abs(z(i,k)-zs(i)-pblh(i)) < (zf(i,k)-zf(i,k+1))*0.5_r8) pblt(i) = k end do end do do k = 1,pver do i = 1,ncol q(i,k) = qh(i,k,1) s(i,k) = t(i,k) + (grav/cpres)*z(i,k) tp(i,k)=0.0_r8 shat(i,k) = s(i,k) qhat(i,k) = q(i,k) end do end do do i = 1,ncol capeg(i) = 0._r8 lclg(i) = 1 lelg(i) = pver maxg(i) = 1 tlg(i) = 400._r8 dsubcld(i) = 0._r8 end do call buoyan_dilute(lchnk ,ncol , & q ,t ,p ,z ,pf , & tp ,qstp ,tl ,rl ,cape , & pblt ,lcl ,lel ,lon ,maxi , & rgas ,grav ,cpres ,msg , & tpert ) lengath = 0 do i=1,ncol if (cape(i) > capelmt) then lengath = lengath + 1 index(lengath) = i end if end do if (lengath.eq.0) return do ii=1,lengath i=index(ii) ideep(ii)=i end do do k = 1,pver do i = 1,lengath dp(i,k) = 0.01_r8*dpp(ideep(i),k) qg(i,k) = q(ideep(i),k) tg(i,k) = t(ideep(i),k) pg(i,k) = p(ideep(i),k) zg(i,k) = z(ideep(i),k) sg(i,k) = s(ideep(i),k) tpg(i,k) = tp(ideep(i),k) zfg(i,k) = zf(ideep(i),k) qstpg(i,k) = qstp(ideep(i),k) ug(i,k) = 0._r8 vg(i,k) = 0._r8 end do end do do i = 1,lengath zfg(i,pver+1) = zf(ideep(i),pver+1) end do do i = 1,lengath capeg(i) = cape(ideep(i)) lclg(i) = lcl(ideep(i)) lelg(i) = lel(ideep(i)) maxg(i) = maxi(ideep(i)) tlg(i) = tl(ideep(i)) landfracg(i) = landfrac(ideep(i)) end do do k = msg + 1,pver do i = 1,lengath if (k >= maxg(i)) then dsubcld(i) = dsubcld(i) + dp(i,k) end if end do end do do k = msg + 2,pver do i = 1,lengath sdifr = 0._r8 qdifr = 0._r8 if (sg(i,k) > 0._r8 .or. sg(i,k-1) > 0._r8) & sdifr = abs((sg(i,k)-sg(i,k-1))/max(sg(i,k-1),sg(i,k))) if (qg(i,k) > 0._r8 .or. qg(i,k-1) > 0._r8) & qdifr = abs((qg(i,k)-qg(i,k-1))/max(qg(i,k-1),qg(i,k))) if (sdifr > 1.E-6_r8) then shat(i,k) = log(sg(i,k-1)/sg(i,k))*sg(i,k-1)*sg(i,k)/(sg(i,k-1)-sg(i,k)) else shat(i,k) = 0.5_r8* (sg(i,k)+sg(i,k-1)) end if if (qdifr > 1.E-6_r8) then qhat(i,k) = log(qg(i,k-1)/qg(i,k))*qg(i,k-1)*qg(i,k)/(qg(i,k-1)-qg(i,k)) else qhat(i,k) = 0.5_r8* (qg(i,k)+qg(i,k-1)) end if end do end do call cldprp(lchnk , & qg ,tg ,ug ,vg ,pg , & zg ,sg ,mu ,eu ,du , & md ,ed ,sd ,qd ,mc , & qu ,su ,zfg ,qs ,hmn , & hsat ,shat ,qlg , & cmeg ,maxg ,lelg ,jt ,jlcl , & maxg ,j0 ,jd ,rl ,lengath , & rgas ,grav ,cpres ,msg , & pflxg ,evpg ,cug ,rprdg ,limcnv ,landfracg) do k = msg + 1,pver do i = 1,lengath du (i,k) = du (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k) eu (i,k) = eu (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k) ed (i,k) = ed (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k) cug (i,k) = cug (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k) cmeg (i,k) = cmeg (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k) rprdg(i,k) = rprdg(i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k) evpg (i,k) = evpg (i,k)* (zfg(i,k)-zfg(i,k+1))/dp(i,k) end do end do call closure(lchnk , & qg ,tg ,pg ,zg ,sg , & tpg ,qs ,qu ,su ,mc , & du ,mu ,md ,qd ,sd , & qhat ,shat ,dp ,qstpg ,zfg , & qlg ,dsubcld ,mb ,capeg ,tlg , & lclg ,lelg ,jt ,maxg ,1 , & lengath ,rgas ,grav ,cpres ,rl , & msg ,capelmt ) do i=1,lengath mumax(i) = 0 end do do k=msg + 2,pver do i=1,lengath mumax(i) = max(mumax(i), mu(i,k)/dp(i,k)) end do end do do i=1,lengath if (mumax(i) > 0._r8) then mb(i) = min(mb(i),0.5_r8/(delt*mumax(i))) else mb(i) = 0._r8 endif end do if (no_deep_pbl) then do i=1,lengath if (zm(ideep(i),jt(i)) < pblh(ideep(i))) mb(i) = 0 end do end if do k=msg+1,pver do i=1,lengath mu (i,k) = mu (i,k)*mb(i) md (i,k) = md (i,k)*mb(i) mc (i,k) = mc (i,k)*mb(i) du (i,k) = du (i,k)*mb(i) eu (i,k) = eu (i,k)*mb(i) ed (i,k) = ed (i,k)*mb(i) cmeg (i,k) = cmeg (i,k)*mb(i) rprdg(i,k) = rprdg(i,k)*mb(i) cug (i,k) = cug (i,k)*mb(i) evpg (i,k) = evpg (i,k)*mb(i) pflxg(i,k+1)= pflxg(i,k+1)*mb(i)*100._r8/grav end do end do call q1q2_pjr(lchnk , & dqdt ,dsdt ,qg ,qs ,qu , & su ,du ,qhat ,shat ,dp , & mu ,md ,sd ,qd ,qlg , & dsubcld ,jt ,maxg ,1 ,lengath , & cpres ,rl ,msg , & dlg ,evpg ,cug ) do k = msg + 1,pver !DIR$ CONCURRENT do i = 1,lengath q(ideep(i),k) = qh(ideep(i),k,1) + 2._r8*delt*dqdt(i,k) qtnd(ideep(i),k) = dqdt (i,k) cme (ideep(i),k) = cmeg (i,k) rprd(ideep(i),k) = rprdg(i,k) zdu (ideep(i),k) = du (i,k) mcon(ideep(i),k) = mc (i,k) heat(ideep(i),k) = dsdt (i,k)*cpres dlf (ideep(i),k) = dlg (i,k) pflx(ideep(i),k) = pflxg(i,k) ql (ideep(i),k) = qlg (i,k) end do end do !DIR$ CONCURRENT do i = 1,lengath jctop(ideep(i)) = jt(i) jcbot(ideep(i)) = maxg(i) pflx(ideep(i),pverp) = pflxg(i,pverp) end do do k = pver,msg + 1,-1 do i = 1,ncol prec(i) = prec(i) - dpp(i,k)* (q(i,k)-qh(i,k,1)) - dpp(i,k)*dlf(i,k)*2._r8*delt end do end do do i = 1,ncol prec(i) = rgrav*max(prec(i),0._r8)/ (2._r8*delt)/1000._r8 end do do k = 1, pver do i = 1, ncol rliq(i) = rliq(i) + dlf(i,k)*dpp(i,k)/gravit end do end do rliq(:ncol) = rliq(:ncol) /1000._r8 return end subroutine zm_convr subroutine zm_conv_evap(ncol,lchnk, & t,pmid,pdel,q, & tend_s, tend_s_snwprd, tend_s_snwevmlt, tend_q, & prdprec, cldfrc, deltat, & prec, snow, ntprprd, ntsnprd, flxprec, flxsnow ) use wv_saturation, only: aqsat integer,intent(in) :: ncol, lchnk real(r8),intent(in), dimension(pcols,pver) :: t real(r8),intent(in), dimension(pcols,pver) :: pmid real(r8),intent(in), dimension(pcols,pver) :: pdel real(r8),intent(in), dimension(pcols,pver) :: q real(r8),intent(inout), dimension(pcols,pver) :: tend_s real(r8),intent(inout), dimension(pcols,pver) :: tend_q real(r8),intent(out ), dimension(pcols,pver) :: tend_s_snwprd real(r8),intent(out ), dimension(pcols,pver) :: tend_s_snwevmlt real(r8), intent(in ) :: prdprec(pcols,pver) real(r8), intent(in ) :: cldfrc(pcols,pver) real(r8), intent(in ) :: deltat real(r8), intent(inout) :: prec(pcols) real(r8), intent(out) :: snow(pcols) real(r8) :: est (pcols,pver) real(r8) :: fice (pcols,pver) real(r8) :: fsnow_conv(pcols,pver) real(r8) :: qsat (pcols,pver) real(r8),intent(out) :: flxprec(pcols,pverp) real(r8),intent(out) :: flxsnow(pcols,pverp) real(r8),intent(out) :: ntprprd(pcols,pver) real(r8),intent(out) :: ntsnprd(pcols,pver) real(r8) :: work1 real(r8) :: work2 real(r8) :: evpvint(pcols) real(r8) :: evpprec(pcols) real(r8) :: evpsnow(pcols) real(r8) :: snowmlt(pcols) real(r8) :: flxsntm(pcols) real(r8) :: evplimit real(r8) :: rlat(pcols) integer :: i,k prec(:ncol) = prec(:ncol)*1000._r8 call aqsat (t ,pmid ,est ,qsat ,pcols , & ncol ,pver ,1 ,pver ) call cldwat_fice(ncol, t, fice, fsnow_conv) flxprec(:ncol,1) = 0._r8 flxsnow(:ncol,1) = 0._r8 evpvint(:ncol) = 0._r8 do k = 1, pver do i = 1, ncol if (t(i,k) > tmelt) then flxsntm(i) = 0._r8 snowmlt(i) = flxsnow(i,k) * gravit/ pdel(i,k) else flxsntm(i) = flxsnow(i,k) snowmlt(i) = 0._r8 end if evplimit = max(1._r8 - q(i,k)/qsat(i,k), 0._r8) evpprec(i) = ke * (1._r8 - cldfrc(i,k)) * evplimit * sqrt(flxprec(i,k)) evplimit = max(0._r8, (qsat(i,k)-q(i,k)) / deltat) evplimit = min(evplimit, flxprec(i,k) * gravit / pdel(i,k)) evplimit = min(evplimit, (prec(i) - evpvint(i)) * gravit / pdel(i,k)) evpprec(i) = min(evplimit, evpprec(i)) if (flxprec(i,k) > 0._r8) then work1 = min(max(0._r8,flxsntm(i)/flxprec(i,k)),1._r8) evpsnow(i) = evpprec(i) * work1 else evpsnow(i) = 0._r8 end if evpvint(i) = evpvint(i) + evpprec(i) * pdel(i,k)/gravit ntprprd(i,k) = prdprec(i,k) - evpprec(i) if (flxprec(i,k).gt.0._r8) then work1 = min(max(0._r8,flxsnow(i,k)/flxprec(i,k)),1._r8) else work1 = 0._r8 endif work2 = max(fsnow_conv(i,k), work1) if (snowmlt(i).gt.0._r8) work2 = 0._r8 ntsnprd(i,k) = prdprec(i,k)*work2 - evpsnow(i) - snowmlt(i) tend_s_snwprd (i,k) = prdprec(i,k)*work2*latice tend_s_snwevmlt(i,k) = - ( evpsnow(i) + snowmlt(i) )*latice flxprec(i,k+1) = flxprec(i,k) + ntprprd(i,k) * pdel(i,k)/gravit flxsnow(i,k+1) = flxsnow(i,k) + ntsnprd(i,k) * pdel(i,k)/gravit flxprec(i,k+1) = max(flxprec(i,k+1), 0._r8) flxsnow(i,k+1) = max(flxsnow(i,k+1), 0._r8) tend_s(i,k) =-evpprec(i)*latvap + ntsnprd(i,k)*latice tend_q(i,k) = evpprec(i) end do end do prec(:ncol) = flxprec(:ncol,pver+1) / 1000._r8 snow(:ncol) = flxsnow(:ncol,pver+1) / 1000._r8 end subroutine zm_conv_evap subroutine convtran(lchnk , & doconvtran,q ,ncnst ,mu ,md , & du ,eu ,ed ,dp ,dsubcld , & jt ,mx ,ideep ,il1g ,il2g , & nstep ,fracis ,dqdt ,dpdry ) use shr_kind_mod, only: r8 => shr_kind_r8 use constituents, only: cnst_get_type_byind implicit none integer, intent(in) :: lchnk integer, intent(in) :: ncnst logical, intent(in) :: doconvtran(ncnst) real(r8), intent(in) :: q(pcols,pver,ncnst) real(r8), intent(in) :: mu(pcols,pver) real(r8), intent(in) :: md(pcols,pver) real(r8), intent(in) :: du(pcols,pver) real(r8), intent(in) :: eu(pcols,pver) real(r8), intent(in) :: ed(pcols,pver) real(r8), intent(in) :: dp(pcols,pver) real(r8), intent(in) :: dsubcld(pcols) real(r8), intent(in) :: fracis(pcols,pver,ncnst) integer, intent(in) :: jt(pcols) integer, intent(in) :: mx(pcols) integer, intent(in) :: ideep(pcols) integer, intent(in) :: il1g integer, intent(in) :: il2g integer, intent(in) :: nstep real(r8), intent(in) :: dpdry(pcols,pver) real(r8), intent(out) :: dqdt(pcols,pver,ncnst) integer i integer k integer kbm integer kk integer kkp1 integer km1 integer kp1 integer ktm integer m real(r8) cabv real(r8) cbel real(r8) cdifr real(r8) chat(pcols,pver) real(r8) cond(pcols,pver) real(r8) const(pcols,pver) real(r8) fisg(pcols,pver) real(r8) conu(pcols,pver) real(r8) dcondt(pcols,pver) real(r8) small real(r8) mbsth real(r8) mupdudp real(r8) minc real(r8) maxc real(r8) fluxin real(r8) fluxout real(r8) netflux real(r8) dutmp(pcols,pver) real(r8) eutmp(pcols,pver) real(r8) edtmp(pcols,pver) real(r8) dptmp(pcols,pver) small = 1.e-36_r8 mbsth = 1.e-15_r8 ktm = pver kbm = pver do i = il1g, il2g ktm = min(ktm,jt(i)) kbm = min(kbm,mx(i)) end do do m = 2, ncnst if (doconvtran(m)) then if (cnst_get_type_byind(m).eq.'dry') then do k = 1,pver do i =il1g,il2g dptmp(i,k) = dpdry(i,k) dutmp(i,k) = du(i,k)*dp(i,k)/dpdry(i,k) eutmp(i,k) = eu(i,k)*dp(i,k)/dpdry(i,k) edtmp(i,k) = ed(i,k)*dp(i,k)/dpdry(i,k) end do end do else do k = 1,pver do i =il1g,il2g dptmp(i,k) = dp(i,k) dutmp(i,k) = du(i,k) eutmp(i,k) = eu(i,k) edtmp(i,k) = ed(i,k) end do end do endif do k = 1,pver do i =il1g,il2g const(i,k) = q(ideep(i),k,m) fisg(i,k) = fracis(ideep(i),k,m) end do end do do k = 1,pver km1 = max(1,k-1) do i = il1g, il2g minc = min(const(i,km1),const(i,k)) maxc = max(const(i,km1),const(i,k)) if (minc < 0) then cdifr = 0._r8 else cdifr = abs(const(i,k)-const(i,km1))/max(maxc,small) endif if (cdifr > 1.E-6_r8) then cabv = max(const(i,km1),maxc*1.e-12_r8) cbel = max(const(i,k),maxc*1.e-12_r8) chat(i,k) = log(cabv/cbel)/(cabv-cbel)*cabv*cbel else chat(i,k) = 0.5_r8* (const(i,k)+const(i,km1)) end if conu(i,k) = chat(i,k) cond(i,k) = chat(i,k) dcondt(i,k) = 0._r8 end do end do k = 2 km1 = 1 kk = pver do i = il1g,il2g mupdudp = mu(i,kk) + dutmp(i,kk)*dptmp(i,kk) if (mupdudp > mbsth) then conu(i,kk) = (+eutmp(i,kk)*fisg(i,kk)*const(i,kk)*dptmp(i,kk))/mupdudp endif if (md(i,k) < -mbsth) then cond(i,k) = (-edtmp(i,km1)*fisg(i,km1)*const(i,km1)*dptmp(i,km1))/md(i,k) endif end do do kk = pver-1,1,-1 kkp1 = min(pver,kk+1) do i = il1g,il2g mupdudp = mu(i,kk) + dutmp(i,kk)*dptmp(i,kk) if (mupdudp > mbsth) then conu(i,kk) = ( mu(i,kkp1)*conu(i,kkp1)+eutmp(i,kk)*fisg(i,kk)* & const(i,kk)*dptmp(i,kk) )/mupdudp endif end do end do do k = 3,pver km1 = max(1,k-1) do i = il1g,il2g if (md(i,k) < -mbsth) then cond(i,k) = ( md(i,km1)*cond(i,km1)-edtmp(i,km1)*fisg(i,km1)*const(i,km1) & *dptmp(i,km1) )/md(i,k) endif end do end do do k = ktm,pver km1 = max(1,k-1) kp1 = min(pver,k+1) do i = il1g,il2g fluxin = mu(i,kp1)*conu(i,kp1)+ mu(i,k)*min(chat(i,k),const(i,km1)) & -(md(i,k) *cond(i,k) + md(i,kp1)*min(chat(i,kp1),const(i,kp1))) fluxout = mu(i,k)*conu(i,k) + mu(i,kp1)*min(chat(i,kp1),const(i,k)) & -(md(i,kp1)*cond(i,kp1) + md(i,k)*min(chat(i,k),const(i,k))) netflux = fluxin - fluxout if (abs(netflux) < max(fluxin,fluxout)*1.e-12_r8) then netflux = 0._r8 endif dcondt(i,k) = netflux/dptmp(i,k) end do end do !DIR$ NOINTERCHANGE do k = kbm,pver km1 = max(1,k-1) do i = il1g,il2g if (k == mx(i)) then fluxin = mu(i,k)*min(chat(i,k),const(i,km1)) - md(i,k)*cond(i,k) fluxout = mu(i,k)*conu(i,k) - md(i,k)*min(chat(i,k),const(i,k)) netflux = fluxin - fluxout if (abs(netflux) < max(fluxin,fluxout)*1.e-12_r8) then netflux = 0._r8 endif dcondt(i,k) = netflux/dptmp(i,k) else if (k > mx(i)) then dcondt(i,k) = 0._r8 end if end do end do dqdt(:,:,m) = 0._r8 do k = 1,pver kp1 = min(pver,k+1) !DIR$ CONCURRENT do i = il1g,il2g dqdt(ideep(i),k,m) = dcondt(i,k) end do end do end if end do return end subroutine convtran subroutine momtran(lchnk, ncol, & domomtran,q ,ncnst ,mu ,md , & du ,eu ,ed ,dp ,dsubcld , & jt ,mx ,ideep ,il1g ,il2g , & nstep ,dqdt ,pguall ,pgdall, icwu, icwd, dt, seten ) use shr_kind_mod, only: r8 => shr_kind_r8 implicit none integer, intent(in) :: lchnk integer, intent(in) :: ncol integer, intent(in) :: ncnst logical, intent(in) :: domomtran(ncnst) real(r8), intent(in) :: q(pcols,pver,ncnst) real(r8), intent(in) :: mu(pcols,pver) real(r8), intent(in) :: md(pcols,pver) real(r8), intent(in) :: du(pcols,pver) real(r8), intent(in) :: eu(pcols,pver) real(r8), intent(in) :: ed(pcols,pver) real(r8), intent(in) :: dp(pcols,pver) real(r8), intent(in) :: dsubcld(pcols) real(r8), intent(in) :: dt integer, intent(in) :: jt(pcols) integer, intent(in) :: mx(pcols) integer, intent(in) :: ideep(pcols) integer, intent(in) :: il1g integer, intent(in) :: il2g integer, intent(in) :: nstep real(r8), intent(out) :: dqdt(pcols,pver,ncnst) integer i integer k integer kbm integer kk integer kkp1 integer kkm1 integer km1 integer kp1 integer ktm integer m integer ii real(r8) cabv real(r8) cbel real(r8) cdifr real(r8) chat(pcols,pver) real(r8) cond(pcols,pver) real(r8) const(pcols,pver) real(r8) conu(pcols,pver) real(r8) dcondt(pcols,pver) real(r8) small real(r8) mbsth real(r8) mupdudp real(r8) minc real(r8) maxc real(r8) fluxin real(r8) fluxout real(r8) netflux real(r8) momcu real(r8) momcd real(r8) sum real(r8) sum2 real(r8) mududp(pcols,pver) real(r8) mddudp(pcols,pver) real(r8) pgu(pcols,pver) real(r8) pgd(pcols,pver) real(r8),intent(out) :: pguall(pcols,pver,ncnst) real(r8),intent(out) :: pgdall(pcols,pver,ncnst) real(r8),intent(out) :: icwu(pcols,pver,ncnst) real(r8),intent(out) :: icwd(pcols,pver,ncnst) real(r8),intent(out) :: seten(pcols,pver) real(r8) gseten(pcols,pver) real(r8) mflux(pcols,pverp,ncnst) real(r8) wind0(pcols,pver,ncnst) real(r8) windf(pcols,pver,ncnst) real(r8) fkeb, fket, ketend_cons, ketend, utop, ubot, vtop, vbot, gset2 pguall(:,:,:) = 0.0_r8 pgdall(:,:,:) = 0.0_r8 icwu(:ncol,:,:) = q(:ncol,:,:) icwd(:ncol,:,:) = q(:ncol,:,:) mflux(:,:,:) = 0.0_r8 wind0(:,:,:) = 0.0_r8 windf(:,:,:) = 0.0_r8 seten(:,:) = 0.0_r8 gseten(:,:) = 0.0_r8 momcu = 0.4_r8 momcd = 0.4_r8 small = 1.e-36_r8 mbsth = 1.e-15_r8 ktm = pver kbm = pver do i = il1g, il2g ktm = min(ktm,jt(i)) kbm = min(kbm,mx(i)) end do do m = 1, ncnst if (domomtran(m)) then do k = 1,pver do i =il1g,il2g const(i,k) = q(ideep(i),k,m) wind0(i,k,m) = const(i,k) end do end do do k = 1,pver km1 = max(1,k-1) do i = il1g, il2g chat(i,k) = 0.5_r8* (const(i,k)+const(i,km1)) conu(i,k) = chat(i,k) cond(i,k) = chat(i,k) dcondt(i,k) = 0._r8 end do end do k=1 pgu(:il2g,k) = 0.0_r8 pgd(:il2g,k) = 0.0_r8 do k=2,pver-1 km1 = max(1,k-1) kp1 = min(pver,k+1) do i = il1g,il2g mududp(i,k) = ( mu(i,k) * (const(i,k)- const(i,km1))/dp(i,km1) & + mu(i,kp1) * (const(i,kp1) - const(i,k))/dp(i,k)) pgu(i,k) = - momcu * 0.5_r8 * mududp(i,k) mddudp(i,k) = ( md(i,k) * (const(i,k)- const(i,km1))/dp(i,km1) & + md(i,kp1) * (const(i,kp1) - const(i,k))/dp(i,k)) pgd(i,k) = - momcd * 0.5_r8 * mddudp(i,k) end do end do k = pver km1 = max(1,k-1) do i=il1g,il2g mududp(i,k) = mu(i,k) * (const(i,k)- const(i,km1))/dp(i,km1) pgu(i,k) = - momcu * mududp(i,k) mddudp(i,k) = md(i,k) * (const(i,k)- const(i,km1))/dp(i,km1) pgd(i,k) = - momcd * mddudp(i,k) end do k = 2 km1 = 1 kk = pver kkm1 = max(1,kk-1) do i = il1g,il2g mupdudp = mu(i,kk) + du(i,kk)*dp(i,kk) if (mupdudp > mbsth) then conu(i,kk) = (+eu(i,kk)*const(i,kk)*dp(i,kk)+pgu(i,kk)*dp(i,kk))/mupdudp endif if (md(i,k) < -mbsth) then cond(i,k) = (-ed(i,km1)*const(i,km1)*dp(i,km1))-pgd(i,km1)*dp(i,km1)/md(i,k) endif end do do kk = pver-1,1,-1 kkm1 = max(1,kk-1) kkp1 = min(pver,kk+1) do i = il1g,il2g mupdudp = mu(i,kk) + du(i,kk)*dp(i,kk) if (mupdudp > mbsth) then conu(i,kk) = ( mu(i,kkp1)*conu(i,kkp1)+eu(i,kk)* & const(i,kk)*dp(i,kk)+pgu(i,kk)*dp(i,kk))/mupdudp endif end do end do do k = 3,pver km1 = max(1,k-1) do i = il1g,il2g if (md(i,k) < -mbsth) then cond(i,k) = ( md(i,km1)*cond(i,km1)-ed(i,km1)*const(i,km1) & *dp(i,km1)-pgd(i,km1)*dp(i,km1) )/md(i,k) endif end do end do sum = 0._r8 sum2 = 0._r8 do k = ktm,pver km1 = max(1,k-1) kp1 = min(pver,k+1) do i = il1g,il2g ii = ideep(i) dcondt(i,k) = & +(mu(i,kp1)* (conu(i,kp1)-chat(i,kp1)) & -mu(i,k)* (conu(i,k)-chat(i,k)) & +md(i,kp1)* (cond(i,kp1)-chat(i,kp1)) & -md(i,k)* (cond(i,k)-chat(i,k)) & )/dp(i,k) end do end do !DIR$ NOINTERCHANGE do k = kbm,pver km1 = max(1,k-1) do i = il1g,il2g if (k == mx(i)) then dcondt(i,k) = (1./dp(i,k))* & (-mu(i,k)*(conu(i,k)-chat(i,k)) & -md(i,k)*(cond(i,k)-chat(i,k)) & ) end if end do end do dqdt(:,:,m) = 0._r8 do k = 1,pver do i = il1g,il2g ii = ideep(i) dqdt(ii,k,m) = dcondt(i,k) pguall(ii,k,m) = -pgu(i,k) pgdall(ii,k,m) = -pgd(i,k) icwu(ii,k,m) = conu(i,k) icwd(ii,k,m) = cond(i,k) end do end do do k = ktm,pver do i = il1g,il2g ii = ideep(i) mflux(i,k,m) = & -mu(i,k)* (conu(i,k)-chat(i,k)) & -md(i,k)* (cond(i,k)-chat(i,k)) end do end do do k = ktm,pver do i = il1g,il2g ii = ideep(i) km1 = max(1,k-1) kp1 = k+1 windf(i,k,m) = const(i,k) - (mflux(i,kp1,m) - mflux(i,k,m)) * dt /dp(i,k) end do end do end if end do do k = ktm,pver km1 = max(1,k-1) kp1 = min(pver,k+1) do i = il1g,il2g ii = ideep(i) utop = (wind0(i,k,1)+wind0(i,km1,1))/2. vtop = (wind0(i,k,2)+wind0(i,km1,2))/2. ubot = (wind0(i,kp1,1)+wind0(i,k,1))/2. vbot = (wind0(i,kp1,2)+wind0(i,k,2))/2. fket = utop*mflux(i,k,1) + vtop*mflux(i,k,2) fkeb = ubot*mflux(i,k+1,1) + vbot*mflux(i,k+1,2) ketend_cons = (fket-fkeb)/dp(i,k) ketend = ((windf(i,k,1)**2 + windf(i,k,2)**2) - (wind0(i,k,1)**2 + wind0(i,k,2)**2))*0.5/dt gset2 = ketend_cons - ketend gseten(i,k) = gset2 end do end do do k = 1,pver do i = il1g,il2g ii = ideep(i) seten(ii,k) = gseten(i,k) end do end do return end subroutine momtran subroutine buoyan(lchnk ,ncol , & q ,t ,p ,z ,pf , & tp ,qstp ,tl ,rl ,cape , & pblt ,lcl ,lel ,lon ,mx , & rd ,grav ,cp ,msg , & tpert ) implicit none integer, intent(in) :: lchnk integer, intent(in) :: ncol real(r8), intent(in) :: q(pcols,pver) real(r8), intent(in) :: t(pcols,pver) real(r8), intent(in) :: p(pcols,pver) real(r8), intent(in) :: z(pcols,pver) real(r8), intent(in) :: pf(pcols,pver+1) real(r8), intent(in) :: pblt(pcols) real(r8), intent(in) :: tpert(pcols) real(r8), intent(out) :: tp(pcols,pver) real(r8), intent(out) :: qstp(pcols,pver) real(r8), intent(out) :: tl(pcols) real(r8), intent(out) :: cape(pcols) integer lcl(pcols) integer lel(pcols) integer lon(pcols) integer mx(pcols) real(r8) capeten(pcols,5) real(r8) tv(pcols,pver) real(r8) tpv(pcols,pver) real(r8) buoy(pcols,pver) real(r8) a1(pcols) real(r8) a2(pcols) real(r8) estp(pcols) real(r8) pl(pcols) real(r8) plexp(pcols) real(r8) hmax(pcols) real(r8) hmn(pcols) real(r8) y(pcols) logical plge600(pcols) integer knt(pcols) integer lelten(pcols,5) real(r8) cp real(r8) e real(r8) grav integer i integer k integer msg integer n real(r8) rd real(r8) rl do n = 1,5 do i = 1,ncol lelten(i,n) = pver capeten(i,n) = 0._r8 end do end do do i = 1,ncol lon(i) = pver knt(i) = 0 lel(i) = pver mx(i) = lon(i) cape(i) = 0._r8 hmax(i) = 0._r8 end do tp(:ncol,:) = t(:ncol,:) qstp(:ncol,:) = q(:ncol,:) tv(:ncol,:) = t(:ncol,:) *(1._r8+1.608*q(:ncol,:))/ (1._r8+q(:ncol,:)) tpv(:ncol,:) = tv(:ncol,:) buoy(:ncol,:) = 0._r8 do k = pver,msg + 1,-1 do i = 1,ncol hmn(i) = cp*t(i,k) + grav*z(i,k) + rl*q(i,k) if (k >= nint(pblt(i)) .and. k <= lon(i) .and. hmn(i) > hmax(i)) then hmax(i) = hmn(i) mx(i) = k end if end do end do do i = 1,ncol lcl(i) = mx(i) e = p(i,mx(i))*q(i,mx(i))/ (eps1+q(i,mx(i))) tl(i) = 2840._r8/ (3.5_r8*log(t(i,mx(i)))-log(e)-4.805_r8) + 55._r8 if (tl(i) < t(i,mx(i))) then plexp(i) = (1._r8/ (0.2854_r8* (1._r8-0.28_r8*q(i,mx(i))))) pl(i) = p(i,mx(i))* (tl(i)/t(i,mx(i)))**plexp(i) else tl(i) = t(i,mx(i)) pl(i) = p(i,mx(i)) end if end do do k = pver,msg + 2,-1 do i = 1,ncol if (k <= mx(i) .and. (p(i,k) > pl(i) .and. p(i,k-1) <= pl(i))) then lcl(i) = k - 1 end if end do end do do i = 1,ncol plge600(i) = pl(i).ge.600._r8 end do do k = pver,msg + 1,-1 do i=1,ncol if (k > lcl(i) .and. k <= mx(i) .and. plge600(i)) then tv(i,k) = t(i,k)* (1._r8+1.608_r8*q(i,k))/ (1._r8+q(i,k)) qstp(i,k) = q(i,mx(i)) tp(i,k) = t(i,mx(i))* (p(i,k)/p(i,mx(i)))**(0.2854_r8* (1._r8-0.28_r8*q(i,mx(i)))) tpv(i,k) = (tp(i,k)+tpert(i))*(1._r8+1.608_r8*q(i,mx(i)))/ (1._r8+q(i,mx(i))) buoy(i,k) = tpv(i,k) - tv(i,k) + tiedke_add end if end do end do do k = pver,msg + 1,-1 do i=1,ncol if (k == lcl(i) .and. plge600(i)) then tv(i,k) = t(i,k)* (1._r8+1.608_r8*q(i,k))/ (1._r8+q(i,k)) qstp(i,k) = q(i,mx(i)) tp(i,k) = tl(i)* (p(i,k)/pl(i))**(0.2854_r8* (1._r8-0.28_r8*qstp(i,k))) estp(i) = c1*exp((c2* (tp(i,k)-tfreez))/((tp(i,k)-tfreez)+c3)) qstp(i,k) = eps1*estp(i)/ (p(i,k)-estp(i)) a1(i) = cp / rl + qstp(i,k) * (1._r8+ qstp(i,k) / eps1) * rl * eps1 / & (rd * tp(i,k) ** 2) a2(i) = .5_r8* (qstp(i,k)* (1._r8+2._r8/eps1*qstp(i,k))* & (1._r8+qstp(i,k)/eps1)*eps1**2*rl*rl/ & (rd**2*tp(i,k)**4)-qstp(i,k)* & (1._r8+qstp(i,k)/eps1)*2._r8*eps1*rl/ & (rd*tp(i,k)**3)) a1(i) = 1._r8/a1(i) a2(i) = -a2(i)*a1(i)**3 y(i) = q(i,mx(i)) - qstp(i,k) tp(i,k) = tp(i,k) + a1(i)*y(i) + a2(i)*y(i)**2 estp(i) = c1*exp((c2* (tp(i,k)-tfreez))/ ((tp(i,k)-tfreez)+c3)) qstp(i,k) = eps1*estp(i) / (p(i,k)-estp(i)) tpv(i,k) = (tp(i,k)+tpert(i))* (1._r8+1.608_r8*qstp(i,k)) / (1._r8+q(i,mx(i))) buoy(i,k) = tpv(i,k) - tv(i,k) + tiedke_add end if end do end do do k = pver - 1,msg + 1,-1 do i=1,ncol if (k < lcl(i) .and. plge600(i)) then tv(i,k) = t(i,k)* (1._r8+1.608_r8*q(i,k))/ (1._r8+q(i,k)) qstp(i,k) = qstp(i,k+1) tp(i,k) = tp(i,k+1)* (p(i,k)/p(i,k+1))**(0.2854_r8* (1._r8-0.28_r8*qstp(i,k))) estp(i) = c1*exp((c2* (tp(i,k)-tfreez))/((tp(i,k)-tfreez)+c3)) qstp(i,k) = eps1*estp(i)/ (p(i,k)-estp(i)) a1(i) = cp/rl + qstp(i,k)* (1._r8+qstp(i,k)/eps1)*rl*eps1/ (rd*tp(i,k)**2) a2(i) = .5_r8* (qstp(i,k)* (1._r8+2._r8/eps1*qstp(i,k))* & (1._r8+qstp(i,k)/eps1)*eps1**2*rl*rl/ & (rd**2*tp(i,k)**4)-qstp(i,k)* & (1._r8+qstp(i,k)/eps1)*2._r8*eps1*rl/ & (rd*tp(i,k)**3)) a1(i) = 1._r8/a1(i) a2(i) = -a2(i)*a1(i)**3 y(i) = qstp(i,k+1) - qstp(i,k) tp(i,k) = tp(i,k) + a1(i)*y(i) + a2(i)*y(i)**2 estp(i) = c1*exp((c2* (tp(i,k)-tfreez))/ ((tp(i,k)-tfreez)+c3)) qstp(i,k) = eps1*estp(i)/ (p(i,k)-estp(i)) tpv(i,k) = (tp(i,k)+tpert(i))* (1._r8+1.608_r8*qstp(i,k))/(1._r8+q(i,mx(i))) buoy(i,k) = tpv(i,k) - tv(i,k) + tiedke_add end if end do end do do k = msg + 2,pver do i = 1,ncol if (k < lcl(i) .and. plge600(i)) then if (buoy(i,k+1) > 0._r8 .and. buoy(i,k) <= 0._r8) then knt(i) = min(5,knt(i) + 1) lelten(i,knt(i)) = k end if end if end do end do do n = 1,5 do k = msg + 1,pver do i = 1,ncol if (plge600(i) .and. k <= mx(i) .and. k > lelten(i,n)) then capeten(i,n) = capeten(i,n) + rd*buoy(i,k)*log(pf(i,k+1)/pf(i,k)) end if end do end do end do do n = 1,5 do i = 1,ncol if (capeten(i,n) > cape(i)) then cape(i) = capeten(i,n) lel(i) = lelten(i,n) end if end do end do do i = 1,ncol cape(i) = max(cape(i), 0._r8) end do return end subroutine buoyan subroutine cldprp(lchnk , & q ,t ,u ,v ,p , & z ,s ,mu ,eu ,du , & md ,ed ,sd ,qd ,mc , & qu ,su ,zf ,qst ,hmn , & hsat ,shat ,ql , & cmeg ,jb ,lel ,jt ,jlcl , & mx ,j0 ,jd ,rl ,il2g , & rd ,grav ,cp ,msg , & pflx ,evp ,cu ,rprd ,limcnv ,landfrac) implicit none integer, intent(in) :: lchnk real(r8), intent(in) :: q(pcols,pver) real(r8), intent(in) :: t(pcols,pver) real(r8), intent(in) :: p(pcols,pver) real(r8), intent(in) :: z(pcols,pver) real(r8), intent(in) :: s(pcols,pver) real(r8), intent(in) :: zf(pcols,pverp) real(r8), intent(in) :: u(pcols,pver) real(r8), intent(in) :: v(pcols,pver) real(r8), intent(in) :: landfrac(pcols) integer, intent(in) :: jb(pcols) integer, intent(in) :: lel(pcols) integer, intent(out) :: jt(pcols) integer, intent(out) :: jlcl(pcols) integer, intent(in) :: mx(pcols) integer, intent(out) :: j0(pcols) integer, intent(out) :: jd(pcols) integer, intent(in) :: limcnv integer, intent(in) :: il2g integer, intent(in) :: msg real(r8), intent(in) :: rl real(r8), intent(in) :: shat(pcols,pver) real(r8), intent(out) :: rprd(pcols,pver) real(r8), intent(out) :: du(pcols,pver) real(r8), intent(out) :: ed(pcols,pver) real(r8), intent(out) :: eu(pcols,pver) real(r8), intent(out) :: hmn(pcols,pver) real(r8), intent(out) :: hsat(pcols,pver) real(r8), intent(out) :: mc(pcols,pver) real(r8), intent(out) :: md(pcols,pver) real(r8), intent(out) :: mu(pcols,pver) real(r8), intent(out) :: pflx(pcols,pverp) real(r8), intent(out) :: qd(pcols,pver) real(r8), intent(out) :: ql(pcols,pver) real(r8), intent(out) :: qst(pcols,pver) real(r8), intent(out) :: qu(pcols,pver) real(r8), intent(out) :: sd(pcols,pver) real(r8), intent(out) :: su(pcols,pver) real(r8) rd real(r8) grav real(r8) cp real(r8) gamma(pcols,pver) real(r8) dz(pcols,pver) real(r8) iprm(pcols,pver) real(r8) hu(pcols,pver) real(r8) hd(pcols,pver) real(r8) eps(pcols,pver) real(r8) f(pcols,pver) real(r8) k1(pcols,pver) real(r8) i2(pcols,pver) real(r8) ihat(pcols,pver) real(r8) i3(pcols,pver) real(r8) idag(pcols,pver) real(r8) i4(pcols,pver) real(r8) qsthat(pcols,pver) real(r8) hsthat(pcols,pver) real(r8) gamhat(pcols,pver) real(r8) cu(pcols,pver) real(r8) evp(pcols,pver) real(r8) cmeg(pcols,pver) real(r8) qds(pcols,pver) real(r8) c0mask(pcols) real(r8) hmin(pcols) real(r8) expdif(pcols) real(r8) expnum(pcols) real(r8) ftemp(pcols) real(r8) eps0(pcols) real(r8) rmue(pcols) real(r8) zuef(pcols) real(r8) zdef(pcols) real(r8) epsm(pcols) real(r8) ratmjb(pcols) real(r8) est(pcols) real(r8) totpcp(pcols) real(r8) totevp(pcols) real(r8) alfa(pcols) real(r8) ql1 real(r8) tu real(r8) estu real(r8) qstu real(r8) small real(r8) mdt integer khighest integer klowest integer kount integer i,k logical doit(pcols) logical done(pcols) do i = 1,il2g ftemp(i) = 0._r8 expnum(i) = 0._r8 expdif(i) = 0._r8 c0mask(i) = c0_ocn * (1._r8-landfrac(i)) + c0_lnd * landfrac(i) end do do k = 1,pver do i = 1,il2g dz(i,k) = zf(i,k) - zf(i,k+1) end do end do pflx(:il2g,1) = 0 do k = 1,pver do i = 1,il2g k1(i,k) = 0._r8 i2(i,k) = 0._r8 i3(i,k) = 0._r8 i4(i,k) = 0._r8 mu(i,k) = 0._r8 f(i,k) = 0._r8 eps(i,k) = 0._r8 eu(i,k) = 0._r8 du(i,k) = 0._r8 ql(i,k) = 0._r8 cu(i,k) = 0._r8 evp(i,k) = 0._r8 cmeg(i,k) = 0._r8 qds(i,k) = q(i,k) md(i,k) = 0._r8 ed(i,k) = 0._r8 sd(i,k) = s(i,k) qd(i,k) = q(i,k) mc(i,k) = 0._r8 qu(i,k) = q(i,k) su(i,k) = s(i,k) est(i) = c1*exp((c2* (t(i,k)-tfreez))/((t(i,k)-tfreez)+c3)) if ( p(i,k)-est(i) > 0._r8 ) then qst(i,k) = eps1*est(i)/ (p(i,k)-est(i)) else qst(i,k) = 1.0_r8 end if gamma(i,k) = qst(i,k)*(1._r8 + qst(i,k)/eps1)*eps1*rl/(rd*t(i,k)**2)*rl/cp hmn(i,k) = cp*t(i,k) + grav*z(i,k) + rl*q(i,k) hsat(i,k) = cp*t(i,k) + grav*z(i,k) + rl*qst(i,k) hu(i,k) = hmn(i,k) hd(i,k) = hmn(i,k) rprd(i,k) = 0._r8 end do end do do k=1,msg do i=1,il2g rprd(i,k) = 0._r8 end do end do do i = 1,il2g hsthat(i,msg+1) = hsat(i,msg+1) qsthat(i,msg+1) = qst(i,msg+1) gamhat(i,msg+1) = gamma(i,msg+1) totpcp(i) = 0._r8 totevp(i) = 0._r8 end do do k = msg + 2,pver do i = 1,il2g if (abs(qst(i,k-1)-qst(i,k)) > 1.E-6_r8) then qsthat(i,k) = log(qst(i,k-1)/qst(i,k))*qst(i,k-1)*qst(i,k)/ (qst(i,k-1)-qst(i,k)) else qsthat(i,k) = qst(i,k) end if hsthat(i,k) = cp*shat(i,k) + rl*qsthat(i,k) if (abs(gamma(i,k-1)-gamma(i,k)) > 1.E-6_r8) then gamhat(i,k) = log(gamma(i,k-1)/gamma(i,k))*gamma(i,k-1)*gamma(i,k)/ & (gamma(i,k-1)-gamma(i,k)) else gamhat(i,k) = gamma(i,k) end if end do end do jt(:) = pver do i = 1,il2g jt(i) = max(lel(i),limcnv+1) jt(i) = min(jt(i),pver) jd(i) = pver jlcl(i) = lel(i) hmin(i) = 1.E6_r8 end do do k = msg + 1,pver do i = 1,il2g if (hsat(i,k) <= hmin(i) .and. k >= jt(i) .and. k <= jb(i)) then hmin(i) = hsat(i,k) j0(i) = k end if end do end do do i = 1,il2g j0(i) = min(j0(i),jb(i)-2) j0(i) = max(j0(i),jt(i)+2) j0(i) = min(j0(i),pver) end do do k = msg + 1,pver do i = 1,il2g if (k >= jt(i) .and. k <= jb(i)) then hu(i,k) = hmn(i,mx(i)) + cp*tiedke_add su(i,k) = s(i,mx(i)) + tiedke_add end if end do end do do k = pver - 1,msg + 1,-1 do i = 1,il2g if (k < jb(i) .and. k >= jt(i)) then k1(i,k) = k1(i,k+1) + (hmn(i,mx(i))-hmn(i,k))*dz(i,k) ihat(i,k) = 0.5_r8* (k1(i,k+1)+k1(i,k)) i2(i,k) = i2(i,k+1) + ihat(i,k)*dz(i,k) idag(i,k) = 0.5_r8* (i2(i,k+1)+i2(i,k)) i3(i,k) = i3(i,k+1) + idag(i,k)*dz(i,k) iprm(i,k) = 0.5_r8* (i3(i,k+1)+i3(i,k)) i4(i,k) = i4(i,k+1) + iprm(i,k)*dz(i,k) end if end do end do do i = 1,il2g hmin(i) = 1.E6_r8 end do do k = msg + 1,pver do i = 1,il2g if (k >= j0(i) .and. k <= jb(i) .and. hmn(i,k) <= hmin(i)) then hmin(i) = hmn(i,k) expdif(i) = hmn(i,mx(i)) - hmin(i) end if end do end do do k = msg + 2,pver do i = 1,il2g expnum(i) = 0._r8 ftemp(i) = 0._r8 if (k < jt(i) .or. k >= jb(i)) then k1(i,k) = 0._r8 expnum(i) = 0._r8 else expnum(i) = hmn(i,mx(i)) - (hsat(i,k-1)*(zf(i,k)-z(i,k)) + & hsat(i,k)* (z(i,k-1)-zf(i,k)))/(z(i,k-1)-z(i,k)) end if if ((expdif(i) > 100._r8 .and. expnum(i) > 0._r8) .and. & k1(i,k) > expnum(i)*dz(i,k)) then ftemp(i) = expnum(i)/k1(i,k) f(i,k) = ftemp(i) + i2(i,k)/k1(i,k)*ftemp(i)**2 + & (2._r8*i2(i,k)**2-k1(i,k)*i3(i,k))/k1(i,k)**2* & ftemp(i)**3 + (-5._r8*k1(i,k)*i2(i,k)*i3(i,k)+ & 5._r8*i2(i,k)**3+k1(i,k)**2*i4(i,k))/ & k1(i,k)**3*ftemp(i)**4 f(i,k) = max(f(i,k),0._r8) f(i,k) = min(f(i,k),0.0002_r8) end if end do end do do i = 1,il2g if (j0(i) < jb(i)) then if (f(i,j0(i)) < 1.E-6_r8 .and. f(i,j0(i)+1) > f(i,j0(i))) j0(i) = j0(i) + 1 end if end do do k = msg + 2,pver do i = 1,il2g if (k >= jt(i) .and. k <= j0(i)) then f(i,k) = max(f(i,k),f(i,k-1)) end if end do end do do i = 1,il2g eps0(i) = f(i,j0(i)) eps(i,jb(i)) = eps0(i) end do do k = pver,msg + 1,-1 do i = 1,il2g if (k >= j0(i) .and. k <= jb(i)) then eps(i,k) = f(i,j0(i)) end if end do end do do k = pver,msg + 1,-1 do i = 1,il2g if (k < j0(i) .and. k >= jt(i)) eps(i,k) = f(i,k) end do end do do i = 1,il2g if (eps0(i) > 0._r8) then mu(i,jb(i)) = 1._r8 eu(i,jb(i)) = mu(i,jb(i))/dz(i,jb(i)) end if end do do k = pver,msg + 1,-1 do i = 1,il2g if (eps0(i) > 0._r8 .and. (k >= jt(i) .and. k < jb(i))) then zuef(i) = zf(i,k) - zf(i,jb(i)) rmue(i) = (1._r8/eps0(i))* (exp(eps(i,k+1)*zuef(i))-1._r8)/zuef(i) mu(i,k) = (1._r8/eps0(i))* (exp(eps(i,k )*zuef(i))-1._r8)/zuef(i) eu(i,k) = (rmue(i)-mu(i,k+1))/dz(i,k) du(i,k) = (rmue(i)-mu(i,k))/dz(i,k) end if end do end do khighest = pverp klowest = 1 do i=1,il2g khighest = min(khighest,lel(i)) klowest = max(klowest,jb(i)) end do do k = klowest-1,khighest,-1 !cdir$ ivdep do i = 1,il2g if (k <= jb(i)-1 .and. k >= lel(i) .and. eps0(i) > 0._r8) then if (mu(i,k) < 0.01_r8) then hu(i,k) = hu(i,jb(i)) mu(i,k) = 0._r8 eu(i,k) = 0._r8 du(i,k) = mu(i,k+1)/dz(i,k) else hu(i,k) = mu(i,k+1)/mu(i,k)*hu(i,k+1) + & dz(i,k)/mu(i,k)* (eu(i,k)*hmn(i,k)- du(i,k)*hsat(i,k)) end if end if end do end do do i=1,il2g doit(i) = .true. end do do k=klowest-2,khighest-1,-1 do i=1,il2g if (doit(i) .and. k <= jb(i)-2 .and. k >= lel(i)-1) then if (hu(i,k) <= hsthat(i,k) .and. hu(i,k+1) > hsthat(i,k+1) & .and. mu(i,k) >= 0.02_r8) then if (hu(i,k)-hsthat(i,k) < -2000._r8) then jt(i) = k + 1 doit(i) = .false. else jt(i) = k if (eps0(i) <= 0._r8) doit(i) = .false. end if else if (hu(i,k) > hu(i,jb(i)) .or. mu(i,k) < 0.01_r8) then jt(i) = k + 1 doit(i) = .false. end if end if end do end do do k = pver,msg + 1,-1 !cdir$ ivdep do i = 1,il2g if (k >= lel(i) .and. k <= jt(i) .and. eps0(i) > 0._r8) then mu(i,k) = 0._r8 eu(i,k) = 0._r8 du(i,k) = 0._r8 hu(i,k) = hu(i,jb(i)) end if if (k == jt(i) .and. eps0(i) > 0._r8) then du(i,k) = mu(i,k+1)/dz(i,k) eu(i,k) = 0._r8 mu(i,k) = 0._r8 end if end do end do do i = 1,il2g alfa(i) = 0.1_r8 jt(i) = min(jt(i),jb(i)-1) jd(i) = max(j0(i),jt(i)+1) jd(i) = min(jd(i),jb(i)) hd(i,jd(i)) = hmn(i,jd(i)-1) if (jd(i) < jb(i) .and. eps0(i) > 0._r8) then epsm(i) = eps0(i) md(i,jd(i)) = -alfa(i)*epsm(i)/eps0(i) end if end do do k = msg + 1,pver do i = 1,il2g if ((k > jd(i) .and. k <= jb(i)) .and. eps0(i) > 0._r8) then zdef(i) = zf(i,jd(i)) - zf(i,k) md(i,k) = -alfa(i)/ (2._r8*eps0(i))*(exp(2._r8*epsm(i)*zdef(i))-1._r8)/zdef(i) end if end do end do do k = msg + 1,pver do i = 1,il2g if ((k >= jt(i) .and. k <= jb(i)) .and. eps0(i) > 0._r8 .and. jd(i) < jb(i)) then ratmjb(i) = min(abs(mu(i,jb(i))/md(i,jb(i))),1._r8) md(i,k) = md(i,k)*ratmjb(i) end if end do end do small = 1.e-20_r8 do k = msg + 1,pver do i = 1,il2g if ((k >= jt(i) .and. k <= pver) .and. eps0(i) > 0._r8) then ed(i,k-1) = (md(i,k-1)-md(i,k))/dz(i,k-1) mdt = min(md(i,k),-small) hd(i,k) = (md(i,k-1)*hd(i,k-1) - dz(i,k-1)*ed(i,k-1)*hmn(i,k-1))/mdt end if end do end do do k = msg + 2,pver do i = 1,il2g if ((k >= jd(i) .and. k <= jb(i)) .and. eps0(i) > 0._r8 .and. jd(i) < jb(i)) then qds(i,k) = qsthat(i,k) + gamhat(i,k)*(hd(i,k)-hsthat(i,k))/ & (rl*(1._r8 + gamhat(i,k))) end if end do end do do i = 1,il2g done(i) = .false. end do kount = 0 do k = pver,msg + 2,-1 do i = 1,il2g if (( .not. done(i) .and. k > jt(i) .and. k < jb(i)) .and. eps0(i) > 0._r8) then su(i,k) = mu(i,k+1)/mu(i,k)*su(i,k+1) + & dz(i,k)/mu(i,k)* (eu(i,k)-du(i,k))*s(i,k) qu(i,k) = mu(i,k+1)/mu(i,k)*qu(i,k+1) + dz(i,k)/mu(i,k)* (eu(i,k)*q(i,k)- & du(i,k)*qst(i,k)) tu = su(i,k) - grav/cp*zf(i,k) estu = c1*exp((c2* (tu-tfreez))/ ((tu-tfreez)+c3)) qstu = eps1*estu/ ((p(i,k)+p(i,k-1))/2._r8-estu) if (qu(i,k) >= qstu) then jlcl(i) = k kount = kount + 1 done(i) = .true. end if end if end do if (kount >= il2g) goto 690 end do 690 continue do k = msg + 2,pver do i = 1,il2g if (k == jb(i) .and. eps0(i) > 0._r8) then qu(i,k) = q(i,mx(i)) su(i,k) = (hu(i,k)-rl*qu(i,k))/cp end if if ((k > jt(i) .and. k <= jlcl(i)) .and. eps0(i) > 0._r8) then su(i,k) = shat(i,k) + (hu(i,k)-hsthat(i,k))/(cp* (1._r8+gamhat(i,k))) qu(i,k) = qsthat(i,k) + gamhat(i,k)*(hu(i,k)-hsthat(i,k))/ & (rl* (1._r8+gamhat(i,k))) end if end do end do do k = pver,msg + 2,-1 do i = 1,il2g if (k >= jt(i) .and. k < jb(i) .and. eps0(i) > 0._r8) then cu(i,k) = ((mu(i,k)*su(i,k)-mu(i,k+1)*su(i,k+1))/ & dz(i,k)- (eu(i,k)-du(i,k))*s(i,k))/(rl/cp) if (k == jt(i)) cu(i,k) = 0._r8 cu(i,k) = max(0._r8,cu(i,k)) end if end do end do do k = pver,msg + 2,-1 do i = 1,il2g rprd(i,k) = 0._r8 if (k >= jt(i) .and. k < jb(i) .and. eps0(i) > 0._r8 .and. mu(i,k) >= 0.0_r8) then if (mu(i,k) > 0._r8) then ql1 = 1._r8/mu(i,k)* (mu(i,k+1)*ql(i,k+1)- & dz(i,k)*du(i,k)*ql(i,k+1)+dz(i,k)*cu(i,k)) ql(i,k) = ql1/ (1._r8+dz(i,k)*c0mask(i)) else ql(i,k) = 0._r8 end if totpcp(i) = totpcp(i) + dz(i,k)*(cu(i,k)-du(i,k)*ql(i,k+1)) rprd(i,k) = c0mask(i)*mu(i,k)*ql(i,k) end if end do end do do i = 1,il2g qd(i,jd(i)) = qds(i,jd(i)) sd(i,jd(i)) = (hd(i,jd(i)) - rl*qd(i,jd(i)))/cp end do do k = msg + 2,pver do i = 1,il2g if (k >= jd(i) .and. k < jb(i) .and. eps0(i) > 0._r8) then qd(i,k+1) = qds(i,k+1) evp(i,k) = -ed(i,k)*q(i,k) + (md(i,k)*qd(i,k)-md(i,k+1)*qd(i,k+1))/dz(i,k) evp(i,k) = max(evp(i,k),0._r8) mdt = min(md(i,k+1),-small) sd(i,k+1) = ((rl/cp*evp(i,k)-ed(i,k)*s(i,k))*dz(i,k) + md(i,k)*sd(i,k))/mdt totevp(i) = totevp(i) - dz(i,k)*ed(i,k)*q(i,k) end if end do end do do i = 1,il2g totevp(i) = totevp(i) + md(i,jd(i))*qd(i,jd(i)) - md(i,jb(i))*qd(i,jb(i)) end do if (.false.) then do i = 1,il2g k = jb(i) if (eps0(i) > 0._r8) then evp(i,k) = -ed(i,k)*q(i,k) + (md(i,k)*qd(i,k))/dz(i,k) evp(i,k) = max(evp(i,k),0._r8) totevp(i) = totevp(i) - dz(i,k)*ed(i,k)*q(i,k) end if end do endif do i = 1,il2g totpcp(i) = max(totpcp(i),0._r8) totevp(i) = max(totevp(i),0._r8) end do do k = msg + 2,pver do i = 1,il2g if (totevp(i) > 0._r8 .and. totpcp(i) > 0._r8) then md(i,k) = md (i,k)*min(1._r8, totpcp(i)/(totevp(i)+totpcp(i))) ed(i,k) = ed (i,k)*min(1._r8, totpcp(i)/(totevp(i)+totpcp(i))) evp(i,k) = evp(i,k)*min(1._r8, totpcp(i)/(totevp(i)+totpcp(i))) else md(i,k) = 0._r8 ed(i,k) = 0._r8 evp(i,k) = 0._r8 end if cmeg(i,k) = cu(i,k) - evp(i,k) rprd(i,k) = rprd(i,k)-evp(i,k) end do end do pflx(:il2g,1) = 0._r8 do k = 2,pverp do i = 1,il2g pflx(i,k) = pflx(i,k-1) + rprd(i,k-1)*dz(i,k-1) end do end do do k = msg + 1,pver do i = 1,il2g mc(i,k) = mu(i,k) + md(i,k) end do end do return end subroutine cldprp subroutine closure(lchnk , & q ,t ,p ,z ,s , & tp ,qs ,qu ,su ,mc , & du ,mu ,md ,qd ,sd , & qhat ,shat ,dp ,qstp ,zf , & ql ,dsubcld ,mb ,cape ,tl , & lcl ,lel ,jt ,mx ,il1g , & il2g ,rd ,grav ,cp ,rl , & msg ,capelmt ) implicit none integer, intent(in) :: lchnk real(r8), intent(inout) :: q(pcols,pver) real(r8), intent(inout) :: t(pcols,pver) real(r8), intent(inout) :: p(pcols,pver) real(r8), intent(inout) :: mb(pcols) real(r8), intent(in) :: z(pcols,pver) real(r8), intent(in) :: s(pcols,pver) real(r8), intent(in) :: tp(pcols,pver) real(r8), intent(in) :: qs(pcols,pver) real(r8), intent(in) :: qu(pcols,pver) real(r8), intent(in) :: su(pcols,pver) real(r8), intent(in) :: mc(pcols,pver) real(r8), intent(in) :: du(pcols,pver) real(r8), intent(in) :: mu(pcols,pver) real(r8), intent(in) :: md(pcols,pver) real(r8), intent(in) :: qd(pcols,pver) real(r8), intent(in) :: sd(pcols,pver) real(r8), intent(in) :: qhat(pcols,pver) real(r8), intent(in) :: shat(pcols,pver) real(r8), intent(in) :: dp(pcols,pver) real(r8), intent(in) :: qstp(pcols,pver) real(r8), intent(in) :: zf(pcols,pver+1) real(r8), intent(in) :: ql(pcols,pver) real(r8), intent(in) :: cape(pcols) real(r8), intent(in) :: tl(pcols) real(r8), intent(in) :: dsubcld(pcols) integer, intent(in) :: lcl(pcols) integer, intent(in) :: lel(pcols) integer, intent(in) :: jt(pcols) integer, intent(in) :: mx(pcols) real(r8) dtpdt(pcols,pver) real(r8) dqsdtp(pcols,pver) real(r8) dtmdt(pcols,pver) real(r8) dqmdt(pcols,pver) real(r8) dboydt(pcols,pver) real(r8) thetavp(pcols,pver) real(r8) thetavm(pcols,pver) real(r8) dtbdt(pcols),dqbdt(pcols),dtldt(pcols) real(r8) beta real(r8) capelmt real(r8) cp real(r8) dadt(pcols) real(r8) debdt real(r8) dltaa real(r8) eb real(r8) grav integer i integer il1g integer il2g integer k, kmin, kmax integer msg real(r8) rd real(r8) rl do i = il1g,il2g mb(i) = 0._r8 eb = p(i,mx(i))*q(i,mx(i))/ (eps1+q(i,mx(i))) dtbdt(i) = (1._r8/dsubcld(i))* (mu(i,mx(i))*(shat(i,mx(i))-su(i,mx(i)))+ & md(i,mx(i))* (shat(i,mx(i))-sd(i,mx(i)))) dqbdt(i) = (1._r8/dsubcld(i))* (mu(i,mx(i))*(qhat(i,mx(i))-qu(i,mx(i)))+ & md(i,mx(i))* (qhat(i,mx(i))-qd(i,mx(i)))) debdt = eps1*p(i,mx(i))/ (eps1+q(i,mx(i)))**2*dqbdt(i) dtldt(i) = -2840._r8* (3.5_r8/t(i,mx(i))*dtbdt(i)-debdt/eb)/ & (3.5_r8*log(t(i,mx(i)))-log(eb)-4.805_r8)**2 end do do k = msg + 1,pver do i = il1g,il2g dtmdt(i,k) = 0._r8 dqmdt(i,k) = 0._r8 end do end do do k = msg + 1,pver - 1 do i = il1g,il2g if (k == jt(i)) then dtmdt(i,k) = (1._r8/dp(i,k))*(mu(i,k+1)* (su(i,k+1)-shat(i,k+1)- & rl/cp*ql(i,k+1))+md(i,k+1)* (sd(i,k+1)-shat(i,k+1))) dqmdt(i,k) = (1._r8/dp(i,k))*(mu(i,k+1)* (qu(i,k+1)- & qhat(i,k+1)+ql(i,k+1))+md(i,k+1)*(qd(i,k+1)-qhat(i,k+1))) end if end do end do beta = 0._r8 do k = msg + 1,pver - 1 do i = il1g,il2g if (k > jt(i) .and. k < mx(i)) then dtmdt(i,k) = (mc(i,k)* (shat(i,k)-s(i,k))+mc(i,k+1)* (s(i,k)-shat(i,k+1)))/ & dp(i,k) - rl/cp*du(i,k)*(beta*ql(i,k)+ (1-beta)*ql(i,k+1)) dqmdt(i,k) = (mu(i,k+1)* (qu(i,k+1)-qhat(i,k+1)+cp/rl* (su(i,k+1)-s(i,k)))- & mu(i,k)* (qu(i,k)-qhat(i,k)+cp/rl*(su(i,k)-s(i,k)))+md(i,k+1)* & (qd(i,k+1)-qhat(i,k+1)+cp/rl*(sd(i,k+1)-s(i,k)))-md(i,k)* & (qd(i,k)-qhat(i,k)+cp/rl*(sd(i,k)-s(i,k))))/dp(i,k) + & du(i,k)* (beta*ql(i,k)+(1-beta)*ql(i,k+1)) end if end do end do do k = msg + 1,pver do i = il1g,il2g if (k >= lel(i) .and. k <= lcl(i)) then thetavp(i,k) = tp(i,k)* (1000._r8/p(i,k))** (rd/cp)*(1._r8+1.608_r8*qstp(i,k)-q(i,mx(i))) thetavm(i,k) = t(i,k)* (1000._r8/p(i,k))** (rd/cp)*(1._r8+0.608_r8*q(i,k)) dqsdtp(i,k) = qstp(i,k)* (1._r8+qstp(i,k)/eps1)*eps1*rl/(rd*tp(i,k)**2) dtpdt(i,k) = tp(i,k)/ (1._r8+rl/cp* (dqsdtp(i,k)-qstp(i,k)/tp(i,k)))* & (dtbdt(i)/t(i,mx(i))+rl/cp* (dqbdt(i)/tl(i)-q(i,mx(i))/ & tl(i)**2*dtldt(i))) dboydt(i,k) = ((dtpdt(i,k)/tp(i,k)+1._r8/(1._r8+1.608_r8*qstp(i,k)-q(i,mx(i)))* & (1.608_r8 * dqsdtp(i,k) * dtpdt(i,k) -dqbdt(i))) - (dtmdt(i,k)/t(i,k)+0.608_r8/ & (1._r8+0.608_r8*q(i,k))*dqmdt(i,k)))*grav*thetavp(i,k)/thetavm(i,k) end if end do end do do k = msg + 1,pver do i = il1g,il2g if (k > lcl(i) .and. k < mx(i)) then thetavp(i,k) = tp(i,k)* (1000._r8/p(i,k))** (rd/cp)*(1._r8+0.608_r8*q(i,mx(i))) thetavm(i,k) = t(i,k)* (1000._r8/p(i,k))** (rd/cp)*(1._r8+0.608_r8*q(i,k)) dboydt(i,k) = (dtbdt(i)/t(i,mx(i))+0.608_r8/ (1._r8+0.608_r8*q(i,mx(i)))*dqbdt(i)- & dtmdt(i,k)/t(i,k)-0.608_r8/ (1._r8+0.608_r8*q(i,k))*dqmdt(i,k))* & grav*thetavp(i,k)/thetavm(i,k) end if end do end do dadt(il1g:il2g) = 0._r8 kmin = minval(lel(il1g:il2g)) kmax = maxval(mx(il1g:il2g)) - 1 do k = kmin, kmax do i = il1g,il2g if ( k >= lel(i) .and. k <= mx(i) - 1) then dadt(i) = dadt(i) + dboydt(i,k)* (zf(i,k)-zf(i,k+1)) endif end do end do do i = il1g,il2g dltaa = -1._r8* (cape(i)-capelmt) if (dadt(i) /= 0._r8) mb(i) = max(dltaa/tau/dadt(i),0._r8) end do return end subroutine closure subroutine q1q2_pjr(lchnk , & dqdt ,dsdt ,q ,qs ,qu , & su ,du ,qhat ,shat ,dp , & mu ,md ,sd ,qd ,ql , & dsubcld ,jt ,mx ,il1g ,il2g , & cp ,rl ,msg , & dl ,evp ,cu ) implicit none real(r8), intent(in) :: cp integer, intent(in) :: lchnk integer, intent(in) :: il1g integer, intent(in) :: il2g integer, intent(in) :: msg real(r8), intent(in) :: q(pcols,pver) real(r8), intent(in) :: qs(pcols,pver) real(r8), intent(in) :: qu(pcols,pver) real(r8), intent(in) :: su(pcols,pver) real(r8), intent(in) :: du(pcols,pver) real(r8), intent(in) :: qhat(pcols,pver) real(r8), intent(in) :: shat(pcols,pver) real(r8), intent(in) :: dp(pcols,pver) real(r8), intent(in) :: mu(pcols,pver) real(r8), intent(in) :: md(pcols,pver) real(r8), intent(in) :: sd(pcols,pver) real(r8), intent(in) :: qd(pcols,pver) real(r8), intent(in) :: ql(pcols,pver) real(r8), intent(in) :: evp(pcols,pver) real(r8), intent(in) :: cu(pcols,pver) real(r8), intent(in) :: dsubcld(pcols) real(r8),intent(out) :: dqdt(pcols,pver),dsdt(pcols,pver) real(r8),intent(out) :: dl(pcols,pver) integer kbm integer ktm integer jt(pcols) integer mx(pcols) integer i integer k real(r8) emc real(r8) rl do k = msg + 1,pver do i = il1g,il2g dsdt(i,k) = 0._r8 dqdt(i,k) = 0._r8 dl(i,k) = 0._r8 end do end do ktm = pver kbm = pver do i = il1g, il2g ktm = min(ktm,jt(i)) kbm = min(kbm,mx(i)) end do do k = ktm,pver-1 do i = il1g,il2g emc = -cu (i,k) & +evp(i,k) dsdt(i,k) = -rl/cp*emc & + (+mu(i,k+1)* (su(i,k+1)-shat(i,k+1)) & -mu(i,k)* (su(i,k)-shat(i,k)) & +md(i,k+1)* (sd(i,k+1)-shat(i,k+1)) & -md(i,k)* (sd(i,k)-shat(i,k)) & )/dp(i,k) dqdt(i,k) = emc + & (+mu(i,k+1)* (qu(i,k+1)-qhat(i,k+1)) & -mu(i,k)* (qu(i,k)-qhat(i,k)) & +md(i,k+1)* (qd(i,k+1)-qhat(i,k+1)) & -md(i,k)* (qd(i,k)-qhat(i,k)) & )/dp(i,k) dl(i,k) = du(i,k)*ql(i,k+1) end do end do !DIR$ NOINTERCHANGE! do k = kbm,pver do i = il1g,il2g if (k == mx(i)) then dsdt(i,k) = (1._r8/dsubcld(i))* & (-mu(i,k)* (su(i,k)-shat(i,k)) & -md(i,k)* (sd(i,k)-shat(i,k)) & ) dqdt(i,k) = (1._r8/dsubcld(i))* & (-mu(i,k)*(qu(i,k)-qhat(i,k)) & -md(i,k)*(qd(i,k)-qhat(i,k)) & ) else if (k > mx(i)) then dsdt(i,k) = dsdt(i,k-1) dqdt(i,k) = dqdt(i,k-1) end if end do end do return end subroutine q1q2_pjr subroutine buoyan_dilute(lchnk ,ncol , & q ,t ,p ,z ,pf , & tp ,qstp ,tl ,rl ,cape , & pblt ,lcl ,lel ,lon ,mx , & rd ,grav ,cp ,msg , & tpert ) implicit none integer, intent(in) :: lchnk integer, intent(in) :: ncol real(r8), intent(in) :: q(pcols,pver) real(r8), intent(in) :: t(pcols,pver) real(r8), intent(in) :: p(pcols,pver) real(r8), intent(in) :: z(pcols,pver) real(r8), intent(in) :: pf(pcols,pver+1) real(r8), intent(in) :: pblt(pcols) real(r8), intent(in) :: tpert(pcols) real(r8), intent(out) :: tp(pcols,pver) real(r8), intent(out) :: qstp(pcols,pver) real(r8), intent(out) :: tl(pcols) real(r8), intent(out) :: cape(pcols) integer lcl(pcols) integer lel(pcols) integer lon(pcols) integer mx(pcols) real(r8) capeten(pcols,5) real(r8) tv(pcols,pver) real(r8) tpv(pcols,pver) real(r8) buoy(pcols,pver) real(r8) a1(pcols) real(r8) a2(pcols) real(r8) estp(pcols) real(r8) pl(pcols) real(r8) plexp(pcols) real(r8) hmax(pcols) real(r8) hmn(pcols) real(r8) y(pcols) logical plge600(pcols) integer knt(pcols) integer lelten(pcols,5) real(r8) cp real(r8) e real(r8) grav integer i integer k integer msg integer n real(r8) rd real(r8) rl do n = 1,5 do i = 1,ncol lelten(i,n) = pver capeten(i,n) = 0._r8 end do end do do i = 1,ncol lon(i) = pver knt(i) = 0 lel(i) = pver mx(i) = lon(i) cape(i) = 0._r8 hmax(i) = 0._r8 end do tp(:ncol,:) = t(:ncol,:) qstp(:ncol,:) = q(:ncol,:) tv(:ncol,:) = t(:ncol,:) *(1._r8+1.608_r8*q(:ncol,:))/ (1._r8+q(:ncol,:)) tpv(:ncol,:) = tv(:ncol,:) buoy(:ncol,:) = 0._r8 do k = pver,msg + 1,-1 do i = 1,ncol hmn(i) = cp*t(i,k) + grav*z(i,k) + rl*q(i,k) if (k >= nint(pblt(i)) .and. k <= lon(i) .and. hmn(i) > hmax(i)) then hmax(i) = hmn(i) mx(i) = k end if end do end do do i = 1,ncol lcl(i) = mx(i) tl(i) = t(i,mx(i)) pl(i) = p(i,mx(i)) end do call parcel_dilute(lchnk, ncol, msg, mx, p, t, q, tpert, tp, tpv, qstp, pl, tl, lcl) do i = 1,ncol plge600(i) = pl(i).ge.600._r8 end do do k = pver,msg + 1,-1 do i=1,ncol if (k <= mx(i) .and. plge600(i)) then tv(i,k) = t(i,k)* (1._r8+1.608_r8*q(i,k))/ (1._r8+q(i,k)) buoy(i,k) = tpv(i,k) - tv(i,k) + tiedke_add else qstp(i,k) = q(i,k) tp(i,k) = t(i,k) tpv(i,k) = tv(i,k) endif end do end do do k = msg + 2,pver do i = 1,ncol if (k < lcl(i) .and. plge600(i)) then if (buoy(i,k+1) > 0. .and. buoy(i,k) <= 0._r8) then knt(i) = min(5,knt(i) + 1) lelten(i,knt(i)) = k end if end if end do end do do n = 1,5 do k = msg + 1,pver do i = 1,ncol if (plge600(i) .and. k <= mx(i) .and. k > lelten(i,n)) then capeten(i,n) = capeten(i,n) + rd*buoy(i,k)*log(pf(i,k+1)/pf(i,k)) end if end do end do end do do n = 1,5 do i = 1,ncol if (capeten(i,n) > cape(i)) then cape(i) = capeten(i,n) lel(i) = lelten(i,n) end if end do end do do i = 1,ncol cape(i) = max(cape(i), 0._r8) end do return end subroutine buoyan_dilute subroutine parcel_dilute (lchnk, ncol, msg, klaunch, p, t, q, tpert, tp, tpv, qstp, pl, tl, lcl) implicit none integer, intent(in) :: lchnk integer, intent(in) :: ncol integer, intent(in) :: msg integer, intent(in), dimension(pcols) :: klaunch(pcols) real(r8), intent(in), dimension(pcols,pver) :: p real(r8), intent(in), dimension(pcols,pver) :: t real(r8), intent(in), dimension(pcols,pver) :: q real(r8), intent(in), dimension(pcols) :: tpert real(r8), intent(inout), dimension(pcols,pver) :: tp real(r8), intent(inout), dimension(pcols,pver) :: qstp real(r8), intent(inout), dimension(pcols) :: tl real(r8), intent(inout), dimension(pcols) :: pl integer, intent(inout), dimension(pcols) :: lcl real(r8), intent(out), dimension(pcols,pver) :: tpv real(r8) tmix(pcols,pver) real(r8) qtmix(pcols,pver) real(r8) qsmix(pcols,pver) real(r8) smix(pcols,pver) real(r8) xsh2o(pcols,pver) real(r8) ds_xsh2o(pcols,pver) real(r8) ds_freeze(pcols,pver) real(r8) mp(pcols) real(r8) qtp(pcols) real(r8) sp(pcols) real(r8) sp0(pcols) real(r8) qtp0(pcols) real(r8) mp0(pcols) real(r8) lwmax real(r8) dmpdp real(r8) dmpdz real(r8) dpdz,dzdp real(r8) senv real(r8) qtenv real(r8) penv real(r8) tenv real(r8) new_s real(r8) new_q real(r8) dp real(r8) tfguess real(r8) tscool real(r8) qxsk, qxskp1 real(r8) dsdp, dqtdp, dqxsdp real(r8) slcl,qtlcl,qslcl integer rcall integer nit_lheat integer i,k,ii nit_lheat = 2 dmpdz=-1.e-3_r8 lwmax = 1.e-3_r8 tscool = 0.0_r8 qtmix=0._r8 smix=0._r8 qtenv = 0._r8 senv = 0._r8 tenv = 0._r8 penv = 0._r8 qtp0 = 0._r8 sp0 = 0._r8 mp0 = 0._r8 qtp = 0._r8 sp = 0._r8 mp = 0._r8 new_q = 0._r8 new_s = 0._r8 do k = pver, msg+1, -1 do i=1,ncol if (k == klaunch(i)) then qtp0(i) = q(i,k) sp0(i) = entropy(t(i,k),p(i,k),qtp0(i)) mp0(i) = 1._r8 smix(i,k) = sp0(i) qtmix(i,k) = qtp0(i) tfguess = t(i,k) rcall = 1 call ientropy (rcall,i,lchnk,smix(i,k),p(i,k),qtmix(i,k),tmix(i,k),qsmix(i,k),tfguess) end if if (k < klaunch(i)) then dp = (p(i,k)-p(i,k+1)) qtenv = 0.5_r8*(q(i,k)+q(i,k+1)) tenv = 0.5_r8*(t(i,k)+t(i,k+1)) penv = 0.5_r8*(p(i,k)+p(i,k+1)) senv = entropy(tenv,penv,qtenv) dpdz = -(penv*grav)/(rgas*tenv) dzdp = 1._r8/dpdz dmpdp = dmpdz*dzdp sp(i) = sp(i) - dmpdp*dp*senv qtp(i) = qtp(i) - dmpdp*dp*qtenv mp(i) = mp(i) - dmpdp*dp smix(i,k) = (sp0(i) + sp(i)) / (mp0(i) + mp(i)) qtmix(i,k) = (qtp0(i) + qtp(i)) / (mp0(i) + mp(i)) tfguess = tmix(i,k+1) rcall = 2 call ientropy(rcall,i,lchnk,smix(i,k),p(i,k),qtmix(i,k),tmix(i,k),qsmix(i,k),tfguess) if (qsmix(i,k) <= qtmix(i,k) .and. qsmix(i,k+1) > qtmix(i,k+1)) then lcl(i) = k qxsk = qtmix(i,k) - qsmix(i,k) qxskp1 = qtmix(i,k+1) - qsmix(i,k+1) dqxsdp = (qxsk - qxskp1)/dp pl(i) = p(i,k+1) - qxskp1/dqxsdp dsdp = (smix(i,k) - smix(i,k+1))/dp dqtdp = (qtmix(i,k) - qtmix(i,k+1))/dp slcl = smix(i,k+1) + dsdp* (pl(i)-p(i,k+1)) qtlcl = qtmix(i,k+1) + dqtdp*(pl(i)-p(i,k+1)) tfguess = tmix(i,k) rcall = 3 call ientropy (rcall,i,lchnk,slcl,pl(i),qtlcl,tl(i),qslcl,tfguess) endif end if end do end do xsh2o = 0._r8 ds_xsh2o = 0._r8 ds_freeze = 0._r8 do k = pver, msg+1, -1 do i=1,ncol if (k == klaunch(i)) then tp(i,k) = tmix(i,k) qstp(i,k) = q(i,k) tpv(i,k) = (tp(i,k) + tpert(i)) * (1._r8+1.608_r8*qstp(i,k)) / (1._r8+qstp(i,k)) end if if (k < klaunch(i)) then do ii=0,nit_lheat-1 xsh2o(i,k) = max (0._r8, qtmix(i,k) - qsmix(i,k) - lwmax) ds_xsh2o(i,k) = ds_xsh2o(i,k+1) - cpliq * log (tmix(i,k)/tfreez) * max(0._r8,(xsh2o(i,k)-xsh2o(i,k+1))) if (tmix(i,k) <= tfreez+tscool .and. ds_freeze(i,k+1) == 0._r8) then ds_freeze(i,k) = (latice/tmix(i,k)) * max(0._r8,qtmix(i,k)-qsmix(i,k)-xsh2o(i,k)) end if if (tmix(i,k) <= tfreez+tscool .and. ds_freeze(i,k+1) /= 0._r8) then ds_freeze(i,k) = ds_freeze(i,k+1)+(latice/tmix(i,k)) * max(0._r8,(qsmix(i,k+1)-qsmix(i,k))) end if new_s = smix(i,k) + ds_xsh2o(i,k) + ds_freeze(i,k) new_q = qtmix(i,k) - xsh2o(i,k) tfguess = tmix(i,k) rcall =4 call ientropy (rcall,i,lchnk,new_s, p(i,k), new_q, tmix(i,k), qsmix(i,k), tfguess) end do tp(i,k) = tmix(i,k) if (new_q > qsmix(i,k)) then qstp(i,k) = qsmix(i,k) else qstp(i,k) = new_q end if tpv(i,k) = (tp(i,k)+tpert(i))* (1._r8+1.608_r8*qstp(i,k)) / (1._r8+ new_q) end if end do end do return end subroutine parcel_dilute real(r8) function entropy(TK,p,qtot) real(r8), intent(in) :: p,qtot,TK real(r8) :: qv,qsat,e,esat,L,eref,pref pref = 1000.0_r8 eref = 6.106_r8 L = rl - (cpliq - cpwv)*(TK-tfreez) esat = c1*exp(c2*(TK-tfreez)/(c3+TK-tfreez)) qsat=eps1*esat/(p-esat) qv = min(qtot,qsat) e = qv*p / (eps1 +qv) entropy = (cpres + qtot*cpliq)*log( TK/tfreez) - rgas*log( (p-e)/pref ) + & L*qv/TK - qv*rh2o*log(qv/qsat) return end FUNCTION entropy SUBROUTINE ientropy (rcall,icol,lchnk,s,p,qt,T,qsat,Tfg) integer, intent(in) :: icol, lchnk, rcall real(r8), intent(in) :: s, p, Tfg, qt real(r8), intent(out) :: qsat, T real(r8) :: qv,Ts,dTs,fs1,fs2,esat real(r8) :: pref,eref,L,e real(r8) :: this_lat,this_lon integer :: LOOPMAX,i LOOPMAX = 100 pref = 1000.0_r8 eref = 6.106_r8 Ts = Tfg converge: do i=0, LOOPMAX L = rl - (cpliq - cpwv)*(Ts-tfreez) esat = c1*exp(c2*(Ts-tfreez)/(c3+Ts-tfreez)) qsat = eps1*esat/(p-esat) qv = min(qt,qsat) e = qv*p / (eps1 +qv) fs1 = (cpres + qt*cpliq)*log( Ts/tfreez ) - rgas*log( (p-e)/pref ) + & L*qv/Ts - qv*rh2o*log(qv/qsat) - s L = rl - (cpliq - cpwv)*(Ts-1._r8-tfreez) esat = c1*exp(c2*(Ts-1._r8-tfreez)/(c3+Ts-1._r8-tfreez)) qsat = eps1*esat/(p-esat) qv = min(qt,qsat) e = qv*p / (eps1 +qv) fs2 = (cpres + qt*cpliq)*log( (Ts-1._r8)/tfreez ) - rgas*log( (p-e)/pref ) + & L*qv/(Ts-1._r8) - qv*rh2o*log(qv/qsat) - s dTs = fs1/(fs2 - fs1) Ts = Ts+dTs if (abs(dTs).lt.0.001_r8) exit converge if (i .eq. LOOPMAX - 1) then this_lat = 0. this_lon = 0. write(iulog,*) '*** ZM_CONV: IENTROPY: Failed and about to exit, info follows ****' call wrf_debug(1,iulog) write(iulog,100) 'ZM_CONV: IENTROPY. Details: call#,lchnk,icol= ',rcall,lchnk,icol, & ' lat: ',this_lat,' lon: ',this_lon, & ' P(mb)= ', p, ' Tfg(K)= ', Tfg, ' qt(g/kg) = ', 1000._r8*qt, & ' qsat(g/kg) = ', 1000._r8*qsat,', s(J/kg) = ',s call wrf_debug(1,iulog) write(iulog,*) '*** Please report this crash to Po-Lun.Ma@pnnl.gov ***' call wrf_debug(1,iulog) call endrun('**** ZM_CONV IENTROPY: Tmix did not converge ****') end if enddo converge esat = c1*exp(c2*(Ts-tfreez)/(c3+Ts-tfreez)) qsat=eps1*esat/(p-esat) qv = min(qt,qsat) T = Ts 100 format (A,I1,I4,I4,7(A,F6.2)) return end SUBROUTINE ientropy end module module_cu_camzm