diff --git a/physics/cu_gf_deep.F90 b/physics/cu_gf_deep.F90 index a07523342..f025c4ec0 100644 --- a/physics/cu_gf_deep.F90 +++ b/physics/cu_gf_deep.F90 @@ -28,9 +28,9 @@ module cu_gf_deep real(kind=kind_phys), parameter :: pgcd = 0.1 ! !> aerosol awareness, do not user yet! - integer, parameter :: autoconv=1 - integer, parameter :: aeroevap=1 - real(kind=kind_phys), parameter :: ccnclean=250. + integer, parameter :: autoconv=2 !1 + integer, parameter :: aeroevap=3 !1 + real(kind=kind_phys), parameter :: scav_factor = 0.5 !> still 16 ensembles for clousres integer, parameter:: maxens3=16 @@ -56,6 +56,7 @@ subroutine cu_gf_deep_run( & ,ichoice & ! choice of closure, use "0" for ensemble average ,ipr & ! this flag can be used for debugging prints ,ccn & ! not well tested yet + ,ccnclean & ,dtime & ! dt over which forcing is applied ,imid & ! flag to turn on mid level convection ,kpbl & ! level of boundary layer height @@ -176,15 +177,15 @@ subroutine cu_gf_deep_run( & q,qo,zuo,zdo,zdm real(kind=kind_phys), dimension (its:ite) & ,intent (in ) :: & - dx,ccn,z1,psur,xland + dx,z1,psur,xland real(kind=kind_phys), dimension (its:ite) & ,intent (inout ) :: & - mconv + mconv,ccn real(kind=kind_phys) & ,intent (in ) :: & - dtime + dtime,ccnclean ! @@ -291,7 +292,7 @@ subroutine cu_gf_deep_run( & real(kind=kind_phys), dimension (its:ite) :: & edt,edto,edtm,aa1,aa0,xaa0,hkb, & hkbo,xhkb, & - xmb,pwavo, & + xmb,pwavo,ccnloss, & pwevo,bu,bud,cap_max, & cap_max_increment,closure_n,psum,psumh,sig,sigd real(kind=kind_phys), dimension (its:ite) :: & @@ -305,7 +306,7 @@ subroutine cu_gf_deep_run( & integer :: & iloop,nens3,ki,kk,i,k real(kind=kind_phys) :: & - dz,dzo,mbdt,radius, & + dz,dzo,mbdt,radius,pefc, & zcutdown,depth_min,zkbmax,z_detr,zktop, & dh,cap_maxs,trash,trash2,frh,sig_thresh real(kind=kind_phys) entdo,dp,subin,detdo,entup, & @@ -504,8 +505,8 @@ subroutine cu_gf_deep_run( & ! !--- minimum depth (m), clouds must have ! - depth_min=1000. - if(imid.eq.1)depth_min=500. + depth_min=3000. + if(imid.eq.1)depth_min=2500. ! !--- maximum depth (mb) of capping !--- inversion (larger cap = no convection) @@ -844,8 +845,8 @@ subroutine cu_gf_deep_run( & exit endif enddo - ktop(i)=ktopkeep(i) - if(ierr(i).eq.0)ktop(i)=ktopkeep(i) + !ktop(i)=ktopkeep(i) + !if(ierr(i).eq.0)ktop(i)=ktopkeep(i) enddo do 37 i=its,itf kzdown(i)=0 @@ -947,14 +948,14 @@ subroutine cu_gf_deep_run( & call cup_up_moisture('mid',ierr,zo_cup,qco,qrco,pwo,pwavo, & p_cup,kbcon,ktop,dbyo,clw_all,xland1, & qo,gammao_cup,zuo,qeso_cup,k22,qo_cup,c0, & - zqexec,ccn,rho,c1d,tn_cup,up_massentr,up_massdetr,psum,psumh, & + zqexec,ccn,ccnclean,rho,c1d,tn_cup,autoconv,up_massentr,up_massdetr,psum,psumh, & 1,itf,ktf, & its,ite, kts,kte) else call cup_up_moisture('deep',ierr,zo_cup,qco,qrco,pwo,pwavo, & p_cup,kbcon,ktop,dbyo,clw_all,xland1, & qo,gammao_cup,zuo,qeso_cup,k22,qo_cup,c0, & - zqexec,ccn,rho,c1d,tn_cup,up_massentr,up_massdetr,psum,psumh, & + zqexec,ccn,ccnclean,rho,c1d,tn_cup,autoconv,up_massentr,up_massdetr,psum,psumh, & 1,itf,ktf, & its,ite, kts,kte) endif @@ -1022,8 +1023,8 @@ subroutine cu_gf_deep_run( & exit endif enddo - ktop(i)=ktopkeep(i) - if(ierr(i).eq.0)ktop(i)=ktopkeep(i) + !ktop(i)=ktopkeep(i) + !if(ierr(i).eq.0)ktop(i)=ktopkeep(i) enddo 41 continue do i=its,itf @@ -1478,8 +1479,8 @@ subroutine cu_gf_deep_run( & !> - Call cup_dd_edt() to determine downdraft strength in terms of windshear ! call cup_dd_edt(ierr,us,vs,zo,ktop,kbcon,edt,po,pwavo, & - pwo,ccn,pwevo,edtmax,edtmin,edtc,psum,psumh, & - rho,aeroevap,itf,ktf, & + pwo,ccn,ccnclean,pwevo,edtmax,edtmin,edtc,psum,psumh, & + rho,aeroevap,pefc,itf,ktf, & its,ite, kts,kte) do i=its,itf if(ierr(i)/=0)cycle @@ -1715,6 +1716,14 @@ subroutine cu_gf_deep_run( & xt(i,k)= dellat(i,k)*mbdt+tn(i,k) xt(i,k)=max(190.,xt(i,k)) enddo + + ! Smooth dellas (HCB) + do k=kts+1,ktf + xt(i,k)=tn(i,k)+0.25*(dellat(i,k-1) + 2.*dellat(i,k) + dellat(i,k+1)) * mbdt + xt(i,k)=max(190.,xt(i,k)) + xq(i,k)=max(1.e-16, qo(i,k)+0.25*(dellaq(i,k-1) + 2.*dellaq(i,k) + dellaq(i,k+1)) * mbdt) + xhe(i,k)=heo(i,k)+0.25*(dellah(i,k-1) + 2.*dellah(i,k) + dellah(i,k+1)) * mbdt + enddo endif enddo do i=its,itf @@ -2019,6 +2028,16 @@ subroutine cu_gf_deep_run( & endif enddo endif + + do i=its,itf + if(ierr(i).eq.0) then + if(aeroevap.gt.1)then + ! aerosol scavagening + ccnloss(i)=ccn(i)*pefc*xmb(i) ! HCB + ccn(i) = ccn(i) - ccnloss(i)*scav_factor + endif + endif + enddo ! ! since kinetic energy is being dissipated, add heating accordingly (from ecmwf) ! @@ -2317,8 +2336,8 @@ end subroutine rain_evap_below_cloudbase subroutine cup_dd_edt(ierr,us,vs,z,ktop,kbcon,edt,p,pwav, & - pw,ccn,pwev,edtmax,edtmin,edtc,psum2,psumh, & - rho,aeroevap,itf,ktf, & + pw,ccn,ccnclean,pwev,edtmax,edtmin,edtc,psum2,psumh, & + rho,aeroevap,pefc,itf,ktf, & its,ite, kts,kte ) implicit none @@ -2336,15 +2355,22 @@ subroutine cup_dd_edt(ierr,us,vs,z,ktop,kbcon,edt,p,pwav, & real(kind=kind_phys), dimension (its:ite,1) & ,intent (out ) :: & edtc + real(kind=kind_phys), intent (out ) :: & + pefc real(kind=kind_phys), dimension (its:ite) & ,intent (out ) :: & edt real(kind=kind_phys), dimension (its:ite) & ,intent (in ) :: & - pwav,pwev,ccn,psum2,psumh,edtmax,edtmin + pwav,pwev,psum2,psumh,edtmax,edtmin integer, dimension (its:ite) & ,intent (in ) :: & ktop,kbcon + real(kind=kind_phys), intent (in ) :: & !HCB + ccnclean + real(kind=kind_phys), dimension (its:ite) & + ,intent (inout ) :: & + ccn integer, dimension (its:ite) & ,intent (inout) :: & ierr @@ -2356,11 +2382,13 @@ subroutine cup_dd_edt(ierr,us,vs,z,ktop,kbcon,edt,p,pwav, & real(kind=kind_phys) einc,pef,pefb,prezk,zkbc real(kind=kind_phys), dimension (its:ite) :: & vshear,sdp,vws - real(kind=kind_phys) :: prop_c,pefc,aeroadd,alpha3,beta3 - prop_c=8. !10.386 - alpha3 = 1.9 - beta3 = -1.13 + real(kind=kind_phys) :: prop_c,aeroadd,alpha3,beta3 + prop_c=0. !10.386 + alpha3 = 0.75 + beta3 = -0.15 pefc=0. + pefb=0. + pef=0. ! !--- determine downdraft strength in terms of windshear @@ -2410,18 +2438,23 @@ subroutine cup_dd_edt(ierr,us,vs,z,ktop,kbcon,edt,p,pwav, & pefb=1./(1.+prezk) if(pefb.gt.0.9)pefb=0.9 if(pefb.lt.0.1)pefb=0.1 + pefb=pef + edt(i)=1.-.5*(pefb+pef) if(aeroevap.gt.1)then - aeroadd=(ccnclean**beta3)*((psumh(i))**(alpha3-1)) !*1.e6 -! prop_c=.9/aeroadd + aeroadd=0. + if((psumh(i)>0.).and.(psum2(i)>0.))then + aeroadd=((1.e-2*ccnclean)**beta3)*((psumh(i)*1.e0)**(alpha3-1)) prop_c=.5*(pefb+pef)/aeroadd - aeroadd=(ccn(i)**beta3)*((psum2(i))**(alpha3-1)) !*1.e6 + aeroadd=((1.e-2*ccn(i))**beta3)*((psum2(i)*1.e0)**(alpha3-1)) aeroadd=prop_c*aeroadd pefc=aeroadd + if(pefc.gt.0.9)pefc=0.9 if(pefc.lt.0.1)pefc=0.1 edt(i)=1.-pefc if(aeroevap.eq.2)edt(i)=1.-.25*(pefb+pef+2.*pefc) + endif endif @@ -3105,12 +3138,12 @@ subroutine cup_forcing_ens_3d(closure_n,xland,aa0,aa1,xaa0,mbdt,dtime,ierr,ierr2 xf_ens(i,5)=max(0.,xff_ens3(5)) xf_ens(i,6)=max(0.,xff_ens3(6)) xf_ens(i,14)=max(0.,xff_ens3(14)) - a1=max(1.e-5,pr_ens(i,7)) + a1=max(1.e-3,pr_ens(i,7)) xf_ens(i,7)=max(0.,xff_ens3(7)/a1) - a1=max(1.e-5,pr_ens(i,8)) + a1=max(1.e-3,pr_ens(i,8)) xf_ens(i,8)=max(0.,xff_ens3(8)/a1) ! forcing(i,7)=xf_ens(i,8) - a1=max(1.e-5,pr_ens(i,9)) + a1=max(1.e-3,pr_ens(i,9)) xf_ens(i,9)=max(0.,xff_ens3(9)/a1) a1=max(1.e-3,pr_ens(i,15)) xf_ens(i,15)=max(0.,xff_ens3(15)/a1) @@ -3875,7 +3908,7 @@ end subroutine cup_output_ens_3d subroutine cup_up_moisture(name,ierr,z_cup,qc,qrc,pw,pwav, & p_cup,kbcon,ktop,dby,clw_all,xland1, & q,gamma_cup,zu,qes_cup,k22,qe_cup,c0, & - zqexec,ccn,rho,c1d,t, & + zqexec,ccn,ccnclean,rho,c1d,t,autoconv, & up_massentr,up_massdetr,psum,psumh, & itest,itf,ktf, & its,ite, kts,kte ) @@ -3891,6 +3924,7 @@ subroutine cup_up_moisture(name,ierr,z_cup,qc,qrc,pw,pwav, & integer & ,intent (in ) :: & + autoconv, & itest,itf,ktf, & its,ite, kts,kte ! cd= detrainment function @@ -3914,7 +3948,7 @@ subroutine cup_up_moisture(name,ierr,z_cup,qc,qrc,pw,pwav, & ,intent (in ) :: & kbcon,ktop,k22,xland1 real(kind=kind_phys), intent (in ) :: & ! HCB - c0 + c0,ccnclean ! ! input and output ! @@ -3937,9 +3971,9 @@ subroutine cup_up_moisture(name,ierr,z_cup,qc,qrc,pw,pwav, & ,intent (out ) :: & qc,qrc,pw,clw_all real(kind=kind_phys), dimension (its:ite,kts:kte) :: & - qch,qrcb,pwh,clw_allh,c1d,t + qch,qrcb,pwh,clw_allh,c1d,c1d_b,t real(kind=kind_phys), dimension (its:ite) :: & - pwavh + pwavh,kklev real(kind=kind_phys), dimension (its:ite) & ,intent (out ) :: & pwav,psum,psumh @@ -3963,7 +3997,7 @@ subroutine cup_up_moisture(name,ierr,z_cup,qc,qrc,pw,pwav, & ! prop_b(kts:kte)=0 iall=0 - clwdet=50. + c1d_b=c1d bdsp=bdispm ! !--- no precip for small clouds @@ -4016,11 +4050,12 @@ subroutine cup_up_moisture(name,ierr,z_cup,qc,qrc,pw,pwav, & ! ! if(name == "deep" )then do k=k22(i)+1,kbcon(i) - if(t(i,k) > 273.16) then - c0t = c0 - else - c0t = c0 * exp(0.07 * (t(i,k) - 273.16)) - endif + c0t = c0 + !if(t(i,k) > 273.16) then + ! c0t = c0 + !else + ! c0t = c0 * exp(0.07 * (t(i,k) - 273.16)) + !endif qc(i,k)= (qc(i,k-1)*zu(i,k-1)-.5*up_massdetr(i,k-1)* qc(i,k-1)+ & up_massentr(i,k-1)*q(i,k-1)) / & (zu(i,k-1)-.5*up_massdetr(i,k-1)+up_massentr(i,k-1)) @@ -4041,13 +4076,12 @@ subroutine cup_up_moisture(name,ierr,z_cup,qc,qrc,pw,pwav, & !now do the rest ! do k=kbcon(i)+1,ktop(i) - !c0=.004 HCB tuning - !if(t(i,k).lt.270.)c0=.002 HCB tuning - if(t(i,k) > 273.16) then - c0t = c0 - else - c0t = c0 * exp(0.07 * (t(i,k) - 273.16)) - endif + c0t = c0 + !if(t(i,k) > 273.16) then + ! c0t = c0 + !else + ! c0t = c0 * exp(0.07 * (t(i,k) - 273.16)) + !endif denom=zu(i,k-1)-.5*up_massdetr(i,k-1)+up_massentr(i,k-1) if(denom.lt.1.e-16)then ierr(i)=51 @@ -4076,21 +4110,27 @@ subroutine cup_up_moisture(name,ierr,z_cup,qc,qrc,pw,pwav, & (zu(i,k-1)-.5*up_massdetr(i,k-1)+up_massentr(i,k-1)) if(qc(i,k).le.qrch)then - qc(i,k)=qrch + qc(i,k)=qrch+1e-8 endif if(qch(i,k).le.qrch)then - qch(i,k)=qrch + qch(i,k)=qrch+1e-8 endif ! !------- total condensed water before rainout ! + if(name == "deep" )then + clwdet=0.1 ! 05/11/2021 + kklev(i)=maxloc(zu(i,:),1) ! 05/05/2021 + if(k.lt.kklev(i)) clwdet=0. ! 05/05/2021 + else + clwdet=0.1 ! 05/05/2021 + endif + if(k.gt.kbcon(i)+1)c1d(i,k)=clwdet*up_massdetr(i,k-1) + if(k.gt.kbcon(i)+1)c1d_b(i,k)=clwdet*up_massdetr(i,k-1) clw_all(i,k)=max(0.,qc(i,k)-qrch) - qrc(i,k)=max(0.,(qc(i,k)-qrch)) ! /(1.+c0*dz*zu(i,k)) clw_allh(i,k)=max(0.,qch(i,k)-qrch) - qrcb(i,k)=max(0.,(qch(i,k)-qrch)) ! /(1.+c0*dz*zu(i,k)) - if(autoconv.eq.2) then - + if(autoconv.eq.2) then ! ! normalized berry ! @@ -4098,41 +4138,38 @@ subroutine cup_up_moisture(name,ierr,z_cup,qc,qrc,pw,pwav, & ! this will also determine proportionality constant prop_b, which, if applied, ! would give the same results as c0 under these conditions ! - q1=1.e3*rhoc*qrcb(i,k) ! g/m^3 ! g[h2o]/cm^3 - berryc0=q1*q1/(60.0*(5.0 + 0.0366*ccnclean/ & + q1=1.e3*rhoc*clw_allh(i,k) ! g/m^3 ! g[h2o]/cm^3 + berryc0=q1*q1/(60.0*(5.0 + 0.0366*ccnclean/ & ( q1 * bdsp) ) ) !/( - qrcb_h=((qch(i,k)-qrch)*zu(i,k)-qrcb(i,k-1)*(.5*up_massdetr(i,k-1)))/ & - (zu(i,k)+.5*up_massdetr(i,k-1)+c0t*dz*zu(i,k)) - prop_b(k)=c0t*qrcb_h*zu(i,k)/(1.e-3*berryc0) + qrcb_h=(qch(i,k)-qrch)/(1.+(c1d_b(i,k)+c0t)*dz) + prop_b(k)=(c0t*qrcb_h)/max(1.e-8,(1.e-3*berryc0)) + if(prop_b(k)>5.) prop_b(k)=5. pwh(i,k)=zu(i,k)*1.e-3*berryc0*dz*prop_b(k) ! 2. - berryc=qrcb(i,k) - qrcb(i,k)=((qch(i,k)-qrch)*zu(i,k)-pwh(i,k)-qrcb(i,k-1)*(.5*up_massdetr(i,k-1)))/ & - (zu(i,k)+.5*up_massdetr(i,k-1)) + qrcb(i,k)=(max(0.,(qch(i,k)-qrch))*zu(i,k)-pwh(i,k))/(zu(i,k)*(1+c1d_b(i,k)*dz)) if(qrcb(i,k).lt.0.)then - berryc0=(qrcb(i,k-1)*(.5*up_massdetr(i,k-1))-(qch(i,k)-qrch)*zu(i,k))/zu(i,k)*1.e-3*dz*prop_b(k) + berryc0=max(0.,(qch(i,k)-qrch))/(1.e-3*dz*prop_b(k)) pwh(i,k)=zu(i,k)*1.e-3*berryc0*dz*prop_b(k) qrcb(i,k)=0. endif qch(i,k)=qrcb(i,k)+qrch pwavh(i)=pwavh(i)+pwh(i,k) - psumh(i)=psumh(i)+clw_allh(i,k)*zu(i,k) *dz + psumh(i)=psumh(i)+pwh(i,k) ! HCB + !psumh(i)=psumh(i)+clw_allh(i,k)*zu(i,k) *dz ! ! then the real berry ! - q1=1.e3*rhoc*qrc(i,k) ! g/m^3 ! g[h2o]/cm^3 - berryc0=q1*q1/(60.0*(5.0 + 0.0366*ccn(i)/ & + q1=1.e3*rhoc*clw_all(i,k) ! g/m^3 ! g[h2o]/cm^3 + berryc0=q1*q1/(60.0*(5.0 + 0.0366*ccn(i)/ & ( q1 * bdsp) ) ) !/( berryc0=1.e-3*berryc0*dz*prop_b(k) ! 2. - berryc=qrc(i,k) - qrc(i,k)=((qc(i,k)-qrch)*zu(i,k)-zu(i,k)*berryc0-qrc(i,k-1)*(.5*up_massdetr(i,k-1)))/ & - (zu(i,k)+.5*up_massdetr(i,k-1)) + qrc(i,k)=(max(0.,(qc(i,k)-qrch))*zu(i,k)-zu(i,k)*berryc0)/(zu(i,k)*(1+c1d(i,k)*dz)) if(qrc(i,k).lt.0.)then - berryc0=((qc(i,k)-qrch)*zu(i,k)-qrc(i,k-1)*(.5*up_massdetr(i,k-1)))/zu(i,k) + berryc0=max(0.,(qc(i,k)-qrch))/(1.e-3*dz*prop_b(k)) qrc(i,k)=0. endif pw(i,k)=berryc0*zu(i,k) qc(i,k)=qrc(i,k)+qrch -! + ! if not running with berry at all, do the following ! else !c0=.002 @@ -4149,7 +4186,8 @@ subroutine cup_up_moisture(name,ierr,z_cup,qc,qrc,pw,pwav, & if(qrc(i,k).lt.0.)then ! hli new test 02/12/19 qrc(i,k)=0. endif - pw(i,k)=c0t*dz*qrc(i,k)*zu(i,k) + pw(i,k)=c0t*dz*qrc(i,k)*zu(i,k) + !-----srf-08aug2017-----begin ! pw(i,k)=(c1d(i,k)+c0)*dz*max(0.,qrc(i,k) -qrc_crit)! units kg[rain]/kg[air] !-----srf-08aug2017-----end @@ -4161,7 +4199,7 @@ subroutine cup_up_moisture(name,ierr,z_cup,qc,qrc,pw,pwav, & qc(i,k)=qrc(i,k)+qrch endif !autoconv pwav(i)=pwav(i)+pw(i,k) - psum(i)=psum(i)+clw_all(i,k)*zu(i,k) *dz + psum(i)=psum(i)+pw(i,k) ! HCB enddo ! k=kbcon,ktop ! do not include liquid/ice in qc do k=k22(i)+1,ktop(i) @@ -4304,6 +4342,7 @@ subroutine rates_up_pdf(rand_vmas,ipr,name,ktop,ierr,p_cup,entr_rate_2d,hkbo,heo kfinalzu=ktf-2 ktop(i)=kfinalzu 412 continue + ktop(i)=ktopdby(i) ! HCB kklev=min(kklev+3,ktop(i)-2) ! ! at least overshoot by one level diff --git a/physics/cu_gf_driver.F90 b/physics/cu_gf_driver.F90 index 157247f6a..b2eb7660b 100644 --- a/physics/cu_gf_driver.F90 +++ b/physics/cu_gf_driver.F90 @@ -7,7 +7,7 @@ module cu_gf_driver ! DH* TODO: replace constants with arguments to cu_gf_driver_run use physcons , g => con_g, cp => con_cp, xlv => con_hvap, r_v => con_rv use machine , only: kind_phys - use cu_gf_deep, only: cu_gf_deep_run,neg_check,autoconv,aeroevap,fct1d3 + use cu_gf_deep, only: cu_gf_deep_run,neg_check,fct1d3 use cu_gf_sh , only: cu_gf_sh_run implicit none @@ -75,14 +75,14 @@ end subroutine cu_gf_driver_finalize !! !>\section gen_gf_driver GSD GF Cumulus Scheme General Algorithm !> @{ - subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,cactiv, & - forcet,forceqv_spechum,phil,raincv,qv_spechum,t,cld1d, & - us,vs,t2di,w,qv2di_spechum,p2di,psuri, & - hbot,htop,kcnv,xland,hfx2,qfx2,cliw,clcw, & + subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,flag_init, & + cactiv,cactiv_m,forcet,forceqv_spechum,phil,raincv,qv_spechum, & + t,cld1d,us,vs,t2di,w,qv2di_spechum,p2di,psuri, & + hbot,htop,kcnv,xland,hfx2,qfx2,aod_gf,cliw,clcw, & pbl,ud_mf,dd_mf,dt_mf,cnvw_moist,cnvc,imfshalcnv, & flag_for_scnv_generic_tend,flag_for_dcnv_generic_tend, & dtend,dtidx,ntqv,ntiw,ntcw,index_of_temperature,index_of_x_wind, & - index_of_y_wind,index_of_process_scnv,index_of_process_dcnv, & + index_of_y_wind,index_of_process_scnv,index_of_process_dcnv, & ldiag3d,qci_conv,errmsg,errflg) !------------------------------------------------------------- implicit none @@ -97,7 +97,9 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,cactiv, & !integer, parameter :: ichoicem=5 ! 0 2 5 13 integer, parameter :: ichoicem=13 ! 0 2 5 13 integer, parameter :: ichoice_s=3 ! 0 1 2 3 - real(kind=kind_phys), parameter :: aodccn=0.1 + + real(kind=kind_phys), parameter :: aodc0=0.14 + real(kind=kind_phys), parameter :: aodreturn=30. real(kind=kind_phys) :: dts,fpi,fp integer, parameter :: dicycle=0 ! diurnal cycle flag integer, parameter :: dicycle_m=0 !- diurnal cycle flag @@ -105,14 +107,14 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,cactiv, & !------------------------------------------------------------- integer :: its,ite, jts,jte, kts,kte integer, intent(in ) :: im,km,ntracer - logical, intent(in ) :: flag_for_scnv_generic_tend,flag_for_dcnv_generic_tend + logical, intent(in ) :: flag_for_scnv_generic_tend,flag_for_dcnv_generic_tend,flag_init logical, intent(in ) :: ldiag3d real(kind=kind_phys), optional, intent(inout) :: dtend(:,:,:) integer, intent(in) :: dtidx(:,:), & index_of_x_wind, index_of_y_wind, index_of_temperature, & index_of_process_scnv, index_of_process_dcnv, ntqv, ntcw, ntiw - + real(kind=kind_phys), dimension( : , : ), intent(in ) :: forcet,forceqv_spechum,w,phil real(kind=kind_phys), dimension( : , : ), intent(inout ) :: t,us,vs real(kind=kind_phys), dimension( : , : ), intent(inout ) :: qci_conv @@ -133,6 +135,7 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,cactiv, & ! Specific humidity from FV3 real(kind=kind_phys), dimension (:,:), intent(in) :: qv2di_spechum real(kind=kind_phys), dimension (:,:), intent(inout) :: qv_spechum + real(kind=kind_phys), dimension (:), intent(inout) :: aod_gf ! Local water vapor mixing ratios and cloud water mixing ratios real(kind=kind_phys), dimension (im,km) :: qv2di, qv, forceqv, cnvw ! @@ -140,7 +143,7 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,cactiv, & real(kind=kind_phys), intent(in ) :: dt integer, intent(in ) :: imfshalcnv - integer, dimension(:), intent(inout) :: cactiv + integer, dimension(:), intent(inout) :: cactiv,cactiv_m character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -151,6 +154,8 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,cactiv, & real(kind=kind_phys), dimension (im,4) :: rand_clos real(kind=kind_phys), dimension (im,km,11) :: gdc,gdc2 real(kind=kind_phys), dimension (im) :: ht + real(kind=kind_phys), dimension (im) :: ccn_gf,ccn_m + real(kind=kind_phys) :: ccnclean real(kind=kind_phys), dimension (im) :: dx real(kind=kind_phys), dimension (im,km) :: outt,outq,outqc,phh,subm,cupclw,cupclws real(kind=kind_phys), dimension (im,km) :: dhdt,zu,zus,zd,phf,zum,zdm,outum,outvm @@ -179,9 +184,9 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,cactiv, & ! omega (omeg), windspeed (us,vs), and a flag (ierr) to turn off ! convection for this call only and at that particular gridpoint ! - real(kind=kind_phys), dimension (im,km) :: qcheck,zo,t2d,q2d,po,p2d,rhoi + real(kind=kind_phys), dimension (im,km) :: qcheck,zo,t2d,q2d,po,p2d,rhoi,clw_ten real(kind=kind_phys), dimension (im,km) :: tn,qo,tshall,qshall,dz8w,omeg - real(kind=kind_phys), dimension (im) :: ccn,z1,psur,cuten,cutens,cutenm + real(kind=kind_phys), dimension (im) :: z1,psur,cuten,cutens,cutenm real(kind=kind_phys), dimension (im) :: umean,vmean,pmean real(kind=kind_phys), dimension (im) :: xmbs,xmbs2,xmb,xmbm,xmb_dumm,mconv @@ -190,7 +195,7 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,cactiv, & integer :: high_resolution real(kind=kind_phys) :: clwtot,clwtot1,excess,tcrit,tscl_kf,dp,dq,sub_spread,subcenter real(kind=kind_phys) :: dsubclw,dsubclws,dsubclwm,dtime_max,ztm,ztq,hfm,qfm,rkbcon,rktop - real(kind=kind_phys), dimension(km) :: massflx,trcflx_in1,clw_in1,clw_ten1,po_cup + real(kind=kind_phys), dimension(km) :: massflx,trcflx_in1,clw_in1,po_cup ! real(kind=kind_phys), dimension(km) :: trcflx_in2,clw_in2,clw_ten2 real(kind=kind_phys), dimension (im) :: flux_tun,tun_rad_mid,tun_rad_shall,tun_rad_deep character*50 :: ierrc(im),ierrcm(im) @@ -280,7 +285,6 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,cactiv, & ! dx for scale awareness ! dx=40075000./float(lonf) ! tscl_kf=dx/25000. - ccn(its:ite)=150. if (imfshalcnv == 3) then ishallow_g3 = 1 @@ -335,7 +339,24 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,cactiv, & do i= its,itf forcing(i,:)=0. forcing2(i,:)=0. - ccn(i)=100. + ccn_gf(i) = 0. + ccn_m(i) = 0. + + ! set aod and ccn + if (flag_init) then + aod_gf(i)=aodc0 + else + if((cactiv(i).eq.0) .and. (cactiv_m(i).eq.0))then + if(aodc0>aod_gf(i)) aod_gf(i)=aod_gf(i)+((aodc0-aod_gf(i))*(dt/(aodreturn*60))) + if(aod_gf(i)>aodc0) aod_gf(i)=aodc0 + endif + endif + + ccn_gf(i)=max(5., (aod_gf(i)/0.0027)**(1/0.640)) + ccn_m(i)=ccn_gf(i) + + ccnclean=max(5., (aodc0/0.0027)**(1/0.640)) + hbot(i) =kte htop(i) =kts raincv(i)=0. @@ -558,7 +579,8 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,cactiv, & ,dicycle_m & ,ichoicem & ,ipr & - ,ccn & + ,ccn_m & + ,ccnclean & ,dt & ,imid_gf & ,kpbli & @@ -638,7 +660,8 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,cactiv, & ,dicycle & ,ichoice & ,ipr & - ,ccn & + ,ccn_gf & + ,ccnclean & ,dt & ,0 & @@ -761,7 +784,9 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,cactiv, & massflx(:)=0. trcflx_in1(:)=0. clw_in1(:)=0. - clw_ten1(:)=0. + do k=kts,ktf + clw_ten(i, k)=0. + enddo po_cup(:)=0. kstop=kts if(ktopm(i).gt.kts .or. ktop(i).gt.kts)kstop=max(ktopm(i),ktop(i)) @@ -851,20 +876,22 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,cactiv, & massflx (1)=0. trcflx_in1(1)=0. call fct1d3 (kstop,kte,dtime_max,po_cup, & - clw_in1,massflx,trcflx_in1,clw_ten1,g) + clw_in1,massflx,trcflx_in1,clw_ten(i,:),g) do k=1,kstop tem = dt*(outqcs(i,k)*cutens(i)+outqc(i,k)*cuten(i) & +outqcm(i,k)*cutenm(i) & - +clw_ten1(k) & + +clw_ten(i,k) & ) - tem1 = max(0.0, min(1.0, (tcr-t(i,k))*tcrf)) - if (clcw(i,k) .gt. -999.0) then - cliw(i,k) = max(0.,cliw(i,k) + tem * tem1) ! ice - clcw(i,k) = max(0.,clcw(i,k) + tem *(1.0-tem1)) ! water - else - cliw(i,k) = max(0.,cliw(i,k) + tem) - endif + !tem1 = max(0.0, min(1.0, (tcr-t(i,k))*tcrf)) + !if (clcw(i,k) .gt. -999.0) then + ! cliw(i,k) = max(0.,cliw(i,k) + tem * tem1) ! ice + ! clcw(i,k) = max(0.,clcw(i,k) + tem *(1.0-tem1)) ! water + !else + ! cliw(i,k) = max(0.,cliw(i,k) + tem) + !endif + if(t(i,k).le.270.) cliw(i,k) = max(0.,cliw(i,k) + tem) ! HCB + if(t(i,k).gt.270) clcw(i,k) = max(0.,clcw(i,k) + tem) enddo @@ -893,6 +920,29 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,cactiv, & cactiv(i)=0 if(pretm(i).gt.0)raincv(i)=.001*cutenm(i)*pretm(i)*dt endif ! pret > 0 + + if(pretm(i).gt.0)then + cactiv_m(i)=1 + else + cactiv_m(i)=0 + endif + + ! Unify ccn + if(ccn_m(i).lt.ccn_gf(i))then + ccn_gf(i)=ccn_m(i) + endif + + if(ccn_gf(i)<0) ccn_gf(i)=0 + + ! Convert ccn back to aod + aod_gf(i)=0.0027*(ccn_gf(i)**0.64) + if(aod_gf(i)<0.007)then + aod_gf(i)=0.007 + ccn_gf(i)=(aod_gf(i)/0.0027)**(1/0.640) + elseif(aod_gf(i)>aodc0)then + aod_gf(i)=aodc0 + ccn_gf(i)=(aod_gf(i)/0.0027)**(1/0.640) + endif enddo 100 continue ! @@ -958,7 +1008,7 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,cactiv, & if(qidx>=1) then do k=kts,ktf do i=its,itf - tem = (cuten(i)*outq(i,k) + cutenm(i)*outqm(i,k))* dt + tem = (cuten(i)*outq(i,k) + cutenm(i)*outqm(i,k))* dt tem = tem/(1.0_kind_phys+tem) dtend(i,k,qidx) = dtend(i,k,qidx) + tem enddo @@ -969,14 +1019,14 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,cactiv, & do k=kts,ktf do i=its,itf tem_shal = dt*(outqcs(i,k)*cutens(i)+outqcm(i,k)*cutenm(i)) - tem_deep = dt*(outqc(i,k)*cuten(i)+clw_ten1(k)) + tem_deep = dt*(outqc(i,k)*cuten(i)+clw_ten(i,k)) tem = tem_shal+tem_deep tem1 = max(0.0, min(1.0, (tcr-t(i,k))*tcrf)) weight_sum = abs(tem_shal)+abs(tem_deep) if(weight_sum<1e-12) then cycle endif - + if (clcw_save(i,k) .gt. -999.0) then cliw_both = max(0.,cliw_save(i,k) + tem * tem1) - cliw_save(i,k) clcw_both = max(0.,clcw_save(i,k) + tem) - clcw_save(i,k) diff --git a/physics/cu_gf_driver.meta b/physics/cu_gf_driver.meta index 84db197bc..cb7ceabd9 100644 --- a/physics/cu_gf_driver.meta +++ b/physics/cu_gf_driver.meta @@ -119,6 +119,14 @@ kind = kind_phys intent = in optional = F +[flag_init] + standard_name = flag_for_first_time_step + long_name = flag signaling first time step for time integration loop + units = flag + dimensions = () + type = logical + intent = in + optional = F [cactiv] standard_name = conv_activity_counter long_name = convective activity memory @@ -127,6 +135,14 @@ type = integer intent = inout optional = F +[cactiv_m] + standard_name = mid_conv_activity_counter + long_name = mid-level cloud convective activity memory + units = none + dimensions = (horizontal_loop_extent) + type = integer + intent = inout + optional = F [forcet] standard_name = temperature_tendency_due_to_dynamics long_name = temperature tendency due to dynamics only @@ -303,6 +319,15 @@ kind = kind_phys intent = in optional = F +[aod_gf] + standard_name = aod_gf_deep + long_name = aerosol optical depth used in Grell-Freitas Convective Parameterization + units = none + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout + optional = F [cliw] standard_name = ice_water_mixing_ratio_convective_transport_tracer long_name = ratio of mass of ice water to mass of dry air plus vapor (without condensates) in the convectively transported tracer array diff --git a/physics/cu_gf_driver_post.F90 b/physics/cu_gf_driver_post.F90 index 4e172ed5a..eab5eefd6 100644 --- a/physics/cu_gf_driver_post.F90 +++ b/physics/cu_gf_driver_post.F90 @@ -20,7 +20,7 @@ end subroutine cu_gf_driver_post_finalize !> \section arg_table_cu_gf_driver_post_run Argument Table !! \htmlinclude cu_gf_driver_post_run.html !! - subroutine cu_gf_driver_post_run (im, t, q, prevst, prevsq, cactiv, conv_act, errmsg, errflg) + subroutine cu_gf_driver_post_run (im, t, q, prevst, prevsq, cactiv, cactiv_m, conv_act, conv_act_m, errmsg, errflg) use machine, only: kind_phys @@ -33,7 +33,9 @@ subroutine cu_gf_driver_post_run (im, t, q, prevst, prevsq, cactiv, conv_act, er real(kind_phys), intent(out) :: prevst(:,:) real(kind_phys), intent(out) :: prevsq(:,:) integer, intent(in) :: cactiv(:) + integer, intent(in) :: cactiv_m(:) real(kind_phys), intent(out) :: conv_act(:) + real(kind_phys), intent(out) :: conv_act_m(:) character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -53,6 +55,11 @@ subroutine cu_gf_driver_post_run (im, t, q, prevst, prevsq, cactiv, conv_act, er else conv_act(i)=0.0 endif + if (cactiv_m(i).gt.0) then + conv_act_m(i) = conv_act_m(i)+1.0 + else + conv_act_m(i)=0.0 + endif enddo end subroutine cu_gf_driver_post_run diff --git a/physics/cu_gf_driver_post.meta b/physics/cu_gf_driver_post.meta index 152409fbd..6e68d62f9 100644 --- a/physics/cu_gf_driver_post.meta +++ b/physics/cu_gf_driver_post.meta @@ -59,6 +59,14 @@ type = integer intent = in optional = F +[cactiv_m] + standard_name = conv_mid_activity_counter + long_name = midlevel convective activity memory + units = none + dimensions = (horizontal_loop_extent) + type = integer + intent = in + optional = F [conv_act] standard_name = gf_memory_counter long_name = Memory counter for GF @@ -68,6 +76,15 @@ kind = kind_phys intent = out optional = F +[conv_act_m] + standard_name = gf_mid_memory_counter + long_name = Memory counter for GF midlevel + units = none + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out + optional = F [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP diff --git a/physics/cu_gf_driver_pre.F90 b/physics/cu_gf_driver_pre.F90 index 3512f65f9..4d4ae9162 100644 --- a/physics/cu_gf_driver_pre.F90 +++ b/physics/cu_gf_driver_pre.F90 @@ -21,7 +21,8 @@ end subroutine cu_gf_driver_pre_finalize !! \htmlinclude cu_gf_driver_pre_run.html !! subroutine cu_gf_driver_pre_run (flag_init, flag_restart, kdt, fhour, dtp, t, q, prevst, prevsq, & - forcet, forceq, cactiv, conv_act, errmsg, errflg) + forcet, forceq, cactiv, cactiv_m, conv_act, conv_act_m, & + errmsg, errflg) use machine, only: kind_phys @@ -39,7 +40,9 @@ subroutine cu_gf_driver_pre_run (flag_init, flag_restart, kdt, fhour, dtp, t, q, real(kind_phys), intent(out) :: forcet(:,:) real(kind_phys), intent(out) :: forceq(:,:) integer, intent(out) :: cactiv(:) + integer, intent(out) :: cactiv_m(:) real(kind_phys), intent(in) :: conv_act(:) + real(kind_phys), intent(in) :: conv_act_m(:) character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -68,6 +71,7 @@ subroutine cu_gf_driver_pre_run (flag_init, flag_restart, kdt, fhour, dtp, t, q, endif cactiv(:)=nint(conv_act(:)) + cactiv_m(:)=nint(conv_act_m(:)) end subroutine cu_gf_driver_pre_run diff --git a/physics/cu_gf_driver_pre.meta b/physics/cu_gf_driver_pre.meta index 1a7fbe4a3..ee64fb5a9 100644 --- a/physics/cu_gf_driver_pre.meta +++ b/physics/cu_gf_driver_pre.meta @@ -111,6 +111,14 @@ type = integer intent = out optional = F +[cactiv_m] + standard_name = conv_mid_activity_counter + long_name = midlevel convective activity memory + units = none + dimensions = (horizontal_loop_extent) + type = integer + intent = out + optional = F [conv_act] standard_name = gf_memory_counter long_name = Memory counter for GF @@ -120,6 +128,15 @@ kind = kind_phys intent = in optional = F +[conv_act_m] + standard_name = gf_mid_memory_counter + long_name = Memory counter for GF midlevel + units = none + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in + optional = F [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP