From 5ca808ea8242d3fbfdf28828b89647cf781a679e Mon Sep 17 00:00:00 2001 From: "Man.Zhang" Date: Fri, 13 Dec 2019 16:48:14 -0700 Subject: [PATCH 01/14] initialize HWRF sasas scheme using preprocessor directives controlled --- physics/samfdeepcnv.f | 254 +++++++++++++++++++++++++++++++++--------- 1 file changed, 199 insertions(+), 55 deletions(-) diff --git a/physics/samfdeepcnv.f b/physics/samfdeepcnv.f index bb5d5deb1..abd1700c9 100644 --- a/physics/samfdeepcnv.f +++ b/physics/samfdeepcnv.f @@ -137,7 +137,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & & dh, dhh, dp, & dq, dqsdp, dqsdt, dt, & dt2, dtmax, dtmin, - & dxcrtas, dxcrtuf, + & dxcrtas, dxcrtuf, dxcrtuf_hwrf, & dv1h, dv2h, dv3h, & dv1q, dv2q, dv3q, & dz, dz1, e1, edtmax, @@ -196,13 +196,13 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & ! & bb1, bb2, wucb ! c physical parameters -! parameter(grav=grav,asolfac=0.958) + parameter(grav=grav,asolfac=0.89) !HWRF ! parameter(grav=grav) ! parameter(elocp=hvap/cp,el2orc=hvap*hvap/(rv*cp)) ! parameter(c0s=.002,c1=.002,d0=.01) ! parameter(d0=.01) parameter(d0=.001) -! parameter(c0l=c0s*asolfac) + parameter(c0l=c0s*asolfac) ! ! asolfac: aerosol-aware parameter based on Lim (2011) ! asolfac= cx / c0s(=.002) @@ -221,7 +221,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & ! parameter(cinacrmx=-120.,cinacrmn=-120.) parameter(cinacrmx=-120.,cinacrmn=-80.) parameter(bet1=1.875,cd1=.506,f1=2.0,gam1=.5) - parameter(betaw=.03,dxcrtas=8.e3,dxcrtuf=15.e3) + parameter(betaw=.03,dxcrtas=8.e3,dxcrtuf=15.e3,dxcrtuf_hwrf=25.e3) ! ! local variables and arrays real(kind=kind_phys) pfld(im,km), to(im,km), qo(im,km), @@ -267,10 +267,53 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & real(kind=kind_phys) tf, tcr, tcrf parameter (tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf)) +#if HWRF==1 + real*8 :: gasdev,ran1 !zhang + real :: rr !zhang + logical,save :: pert_sas_local !zhang + integer,save :: ens_random_seed_local,env_pp_local !zhang + integer :: ensda_physics_pert !zhang + real,save :: ens_sasamp_local !zhang + data ens_random_seed_local/0/ + data env_pp_local/0/ + CHARACTER(len=3) :: env_memb,env_pp +#endif + ! Initialize CCPP error handling variables errmsg = '' errflg = 0 +#if HWRF==1 + if ( ens_random_seed_local .eq. 0 ) then + CALL nl_get_ensda_physics_pert(1,ensda_physics_pert) + ens_random_seed_local=ens_random_seed + env_pp_local=ensda_physics_pert + pert_sas_local=.false. + ens_sasamp_local=0.0 +! env_pp=1: do physics perturbations for ensda members, ens_random_seed must be 99 + if ( env_pp_local .eq. 1 ) then + if ( ens_random_seed .ne. 99 ) then + pert_sas_local=.true. + ens_sasamp_local=ens_sasamp + else +! ens_random_seed=99 do physics perturbation for ensemble forecasts, env_pp must be zero + ens_random_seed_local=ens_random_seed + pert_sas_local=pert_sas + ens_sasamp_local=ens_sasamp + endif + else + ens_random_seed_local=ens_random_seed + pert_sas_local=pert_sas + ens_sasamp_local=ens_sasamp + endif + print*, "DESAS ==", ens_random_seed_local,pert_sas_local,ens_sasamp_local,ensda_physics_pert + endif +#endif + + +#endif + +#ifndef HWRF_SCALESAS elocp = hvap/cp el2orc = hvap*hvap/(rv*cp) @@ -281,6 +324,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & !> ## Determine whether to perform aerosol transport do_aerosols = (itc > 0) .and. (ntc > 0) .and. (ntr > 0) if (do_aerosols) do_aerosols = (ntr >= itc + ntc - 3) +#endif ! c----------------------------------------------------------------------- !> ## Compute preliminary quantities needed for static, dynamic, and feedback control portions of the algorithm. @@ -328,12 +372,22 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & xpwev(i)= 0. vshear(i) = 0. gdx(i) = sqrt(garea(i)) + +#ifdef HWRF_SCALESAS + scaldfunc(i)=-1.0 ! initialized wang + sigmagfm(i)=-1.0 + sigmuout(i)=-1.0 +#endif enddo ! !> - determine aerosol-aware rain conversion parameter over land do i=1,im if(islimsk(i) == 1) then +#ifdef HWRF_SCALESAS + c0(i) = c0l +#else c0(i) = c0s*asolfac +#endif else c0(i) = c0s endif @@ -366,6 +420,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & dt_mf(i,k) = 0. enddo enddo + if(mp_phys == mp_phys_mg) then do k = 1, km do i = 1, im @@ -398,8 +453,15 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & ! model tunable parameters are all here edtmaxl = .3 edtmaxs = .3 +#ifdef HWRF_SCALESAS + clam = .1 + aafac = .1 + betal = .05 + betas = .05 + evfact = 0.3 + evfactl = 0.3 +#else ! clam = .1 -! aafac = .1 aafac = .05 ! betal = .15 ! betas = .15 @@ -408,12 +470,17 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & ! evef = 0.07 ! evfact = 0.3 ! evfactl = 0.3 +#endif ! +#ifdef HWRF_SCALESAS + crtlamu = 1.0e-4 + cxlamu = 1.0e-3 +#else crtlame = 1.0e-4 - crtlamd = 1.0e-4 -! -! cxlame = 1.0e-3 cxlame = 1.0e-4 +#endif + + crtlamd = 1.0e-4 cxlamd = 1.0e-4 xlamde = 1.0e-4 xlamdd = 1.0e-4 @@ -467,6 +534,9 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & do k = 1, km1 do i=1,im zi(i,k) = 0.5*(zo(i,k)+zo(i,k+1)) +#ifdef HWRF_SCALESAS + xlamue(i,k) = clam / zi(i,k) +#endif enddo enddo c @@ -514,6 +584,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & ! ! initialize tracer variables ! +#ifndef HWRF_SCALESAS do n = 3, ntr+2 kk = n-2 do k = 1, km @@ -527,6 +598,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & enddo enddo enddo +#endif ! !> - Calculate saturation specific humidity and enforce minimum moisture values. do k = 1, km @@ -623,6 +695,8 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo + +#ifndef HWRF_SCALESAS do n = 1, ntr do k = 1, km1 do i=1,im @@ -632,6 +706,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & enddo enddo enddo +#endif c c look for the level of free convection as cloud base c @@ -701,6 +776,18 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & ptem1= .5*(cinpcrmx-cinpcrmn) cinpcr = cinpcrmx - ptem * ptem1 tem1 = pfld(i,kb(i)) - pfld(i,kbcon(i)) +#if HWRF==1 +! randomly perturb the convection trigger +!zzz if( pert_sas_local .and. ens_random_seed_local .gt. 0 ) then + if( pert_sas_local ) then +!zz print*,"ens_random_seed==",ens_random_seed,ens_random_seed_local + ens_random_seed_local=ran1(-ens_random_seed_local)*1000 + rr=2.0*ens_sasamp_local*ran1(-ens_random_seed_local)-ens_sasamp_local +!zz print*, "zhang inde sas=a", cinpcr,ens_sasamp_local,ens_random_seed_local,cinpcr + cinpcr=cinpcr+rr +!zz print*, "zhang inde sas=b", cinpcr,ens_sasamp_local,ens_random_seed_local,cinpcr + endif +#endif if(tem1 > cinpcr) then cnvflg(i) = .false. endif @@ -712,6 +799,8 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & totflg = totflg .and. (.not. cnvflg(i)) enddo if(totflg) return + +#ifndef HWRF_SCALESAS !! ! ! turbulent entrainment rate assumed to be proportional @@ -774,6 +863,8 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo + +#endif c c assume that updraft entrainment rate above cloud base is c same as that at cloud base @@ -783,19 +874,21 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & !! \epsilon = \epsilon_0F_0 + d_1\left(1-RH\right)F_1 !! \f] !! where \f$\epsilon_0\f$ is the cloud base entrainment rate, \f$d_1\f$ is a tunable constant, and \f$F_0=\left(\frac{q_s}{q_{s,b}}\right)^2\f$ and \f$F_1=\left(\frac{q_s}{q_{s,b}}\right)^3\f$ where \f$q_s\f$ and \f$q_{s,b}\f$ are the saturation specific humidities at a given level and cloud base, respectively. The detrainment rate in the cloud is assumed to be equal to the entrainment rate at cloud base. -! do i=1,im -! if(cnvflg(i)) then -! xlamx(i) = xlamue(i,kbcon(i)) -! endif -! enddo -! do k = 2, km1 -! do i=1,im -! if(cnvflg(i).and. -! & (k > kbcon(i) .and. k < kmax(i))) then -! xlamue(i,k) = xlamx(i) -! endif -! enddo -! enddo +#ifdef HWRF_SCALESAS + do i=1,im + if(cnvflg(i)) then + xlamx(i) = xlamue(i,kbcon(i)) + endif + enddo + do k = 2, km1 + do i=1,im + if(cnvflg(i).and. & + & (k > kbcon(i) .and. k < kmax(i))) then + xlamue(i,k) = xlamx(i) + endif + enddo + enddo +#endif c c specify detrainment rate for the updrafts c @@ -805,9 +898,11 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & do k = 1, km1 do i=1,im if(cnvflg(i) .and. k < kmax(i)) then -! xlamud(i,k) = xlamx(i) -! xlamud(i,k) = crtlamd +#ifdef HWRF_SCALESAS + xlamud(i,k) = xlamx(i) +#else xlamud(i,k) = 0.001 * clamt(i) +#endif endif enddo enddo @@ -837,8 +932,10 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & & (k > kbcon(i) .and. k < kmax(i))) then tem = cxlame * frh(i,k) * fent2(i,k) xlamue(i,k) = xlamue(i,k)*fent1(i,k) + tem +#ifndef HWRF_SCALESAS tem1 = cxlamd * frh(i,k) xlamud(i,k) = xlamud(i,k) + tem1 +#endif endif enddo enddo @@ -900,6 +997,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & pwavo(i) = 0. endif enddo +#ifndef HWRF_SCALESAS ! for tracers do n = 1, ntr do i = 1, im @@ -909,6 +1007,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo +#endif c c cloud property is modified by the entrainment process c @@ -939,6 +1038,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo +#ifndef HWRF_SCALESAS do n = 1, ntr do k = 2, km1 do i = 1, im @@ -954,6 +1054,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & enddo enddo enddo +#endif c c taking account into convection inhibition due to existence of c dry layers below cloud base @@ -1023,6 +1124,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & do i = 1, im if(cnvflg(i)) then ! +#ifndef HWRF_SCALESAS if(islimsk(i) == 1) then w1 = w1l w2 = w2l @@ -1051,6 +1153,9 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & cinacr = cinacrmx - tem * tem1 ! ! cinacr = cinacrmx +#else + cinacr = cinacrmx +#endif if(cina(i) < cinacr) cnvflg(i) = .false. endif enddo @@ -1137,14 +1242,11 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & ! k = kbcon(i) dp = 1000. * del(i,k) +#ifndef HWRF_SCALEASA xmbmax(i) = dp / (2. * grav * dt2) -! -! xmbmax(i) = dp / (grav * dt2) -! -! mbdt(i) = 0.1 * dp / grav -! -! tem = dp / (grav * dt2) -! xmbmax(i) = min(tem, xmbmax(i)) +#else + xmbmax(i) = dp / (grav * dt2) +#endif endif enddo c @@ -1184,8 +1286,13 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & c if(k >= kbcon(i) .and. dq > 0.) then etah = .5 * (eta(i,k) + eta(i,k-1)) +#ifndef HWRF_SCALESAS dp = 1000. * del(i,k) +#endif if(ncloud > 0 .and. k > jmin(i)) then +#ifdef HWRF_SCALESAS + dp = 1000. * del(i,k) +#endif ptem = c0t(i,k) + c1 qlk = dq / (eta(i,k) + etah * ptem * dz) dellal(i,k) = etah * c1 * dz * qlk * grav / dp @@ -1357,8 +1464,13 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & c if(dq > 0.) then etah = .5 * (eta(i,k) + eta(i,k-1)) +#ifndef HWRF_SCALESAS dp = 1000. * del(i,k) +#endif if(ncloud > 0) then +#ifdef HWRF_SCALESAS + dp = 1000. * del(i,k) +#endif ptem = c0t(i,k) + c1 qlk = dq / (eta(i,k) + etah * ptem * dz) dellal(i,k) = etah * c1 * dz * qlk * grav / dp @@ -1379,30 +1491,23 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & ! compute updraft velocity square(wu2) !> - Calculate updraft velocity square(wu2) according to Han et al.'s (2017) \cite han_et_al_2017 equation 7. ! -! bb1 = 2. * (1.+bet1*cd1) -! bb2 = 2. / (f1*(1.+gam1)) -! -! bb1 = 3.9 -! bb2 = 0.67 -! -! bb1 = 2.0 -! bb2 = 4.0 -! - bb1 = 4.0 - bb2 = 0.8 + bb1 = 4.0 + bb2 = 0.8 +#ifdef HWRF_SCALESAS + do i = 1, im + if (cnvflg(i)) then + k = kbcon1(i) + tem = po(i,k) / (rd * to(i,k)) + wucb = -0.01 * dot(i,k) / (tem * g) + if(wucb.gt.0.) then + wu2(i,k) = wucb * wucb + else + wu2(i,k) = 0. + endif + endif + enddo +#endif ! -! do i = 1, im -! if (cnvflg(i)) then -! k = kbcon1(i) -! tem = po(i,k) / (rd * to(i,k)) -! wucb = -0.01 * dot(i,k) / (tem * grav) -! if(wucb > 0.) then -! wu2(i,k) = wucb * wucb -! else -! wu2(i,k) = 0. -! endif -! endif -! enddo do k = 2, km1 do i = 1, im if (cnvflg(i)) then @@ -1554,6 +1659,10 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & enddo enddo do i = 1, im +#ifdef HWRF_SCALESAS + beta = betas + if(islimsk(i) == 1) beta = betal +#else betamn = betas if(islimsk(i) == 1) betamn = betal if(ntk > 0) then @@ -1569,6 +1678,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & else beta = betamn endif +#endif if(cnvflg(i)) then dz = (sumx(i)+zi(i,1))/float(kbcon(i)) tem = 1./float(kbcon(i)) @@ -1610,6 +1720,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo ! for tracers +#ifndef HWRF_SCALESAS do n = 1, ntr do i = 1, im if(cnvflg(i)) then @@ -1618,6 +1729,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo +#endif cj !> - Calculate the cloud properties as a parcel descends, modified by entrainment and detrainment. Discretization follows Appendix B of Grell (1993) \cite grell_1993 . do k = km1, 1, -1 @@ -1647,6 +1759,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo +#ifndef HWRF_SCALESAS do n = 1, ntr do k = km1, 1, -1 do i = 1, im @@ -1660,6 +1773,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & enddo enddo enddo +#endif c !> - Compute the amount of moisture that is necessary to keep the downdraft saturated. do k = km1, 1, -1 @@ -1762,6 +1876,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo +#ifndef HWRF_SCALESAS do n = 1, ntr do k = 1, km do i = 1, im @@ -1771,6 +1886,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & enddo enddo enddo +#endif do i = 1, im if(cnvflg(i)) then dp = 1000. * del(i,1) @@ -1784,6 +1900,8 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & & - vo(i,1)) * grav / dp endif enddo + +#ifndef HWRF_SCALESAS do n = 1, ntr do i = 1, im if(cnvflg(i)) then @@ -1793,6 +1911,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo +#endif c c--- changed due to subsidence and entrainment c @@ -1857,6 +1976,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo +#ifndef HWRF_SCALESAS do n = 1, ntr do k = 2, km1 do i = 1, im @@ -1878,6 +1998,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & enddo enddo enddo +#endif c c------- cloud top c @@ -1902,6 +2023,8 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & & qlko_ktcon(i) * grav / dp endif enddo + +#ifndef HWRF_SCALESAS do n = 1, ntr do i = 1, im if(cnvflg(i)) then @@ -1912,6 +2035,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo +#endif c c------- final changed variable per unit mass flux c @@ -1942,6 +2066,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo + c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c c--- the above changed environment is now used to calulate the @@ -2282,8 +2407,10 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & if(cnvflg(i)) then tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i)) dtconv(i) = tem / wc(i) +#ifndef HWRF_SCALESAS tfac = 1. + gdx(i) / 75000. dtconv(i) = tfac * dtconv(i) +#endif dtconv(i) = max(dtconv(i),dtmin) dtconv(i) = min(dtconv(i),dtmax) endif @@ -2326,6 +2453,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & xmb(i) = tfac*betaw*rho*wc(i) endif enddo + !> - For the cases where the quasi-equilibrium assumption of Arakawa-Schubert is valid, first calculate the large scale destabilization as in equation 5 of Pan and Wu (1995) \cite pan_and_wu_1995 : !! \f[ !! \frac{\partial A}{\partial t}_{LS}=\frac{A^+-cA^0}{\Delta t_{LS}} @@ -2366,7 +2494,6 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & tfac = tauadv(i) / dtconv(i) tfac = min(tfac, 1.) xmb(i) = -tfac * fld(i) / xk(i) -! xmb(i) = min(xmb(i),xmbmax(i)) endif enddo !! @@ -2377,7 +2504,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & enddo if(totflg) return !! -! + !> - For scale-aware parameterization, the updraft fraction (sigmagfm) is first computed as a function of the lateral entrainment rate at cloud base (see Han et al.'s (2017) \cite han_et_al_2017 equation 4 and 5), following the study by Grell and Freitas (2014) \cite grell_and_freitas_2014. do i = 1, im if(cnvflg(i)) then @@ -2396,6 +2523,9 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & if (gdx(i) < dxcrtuf) then scaldfunc(i) = (1.-sigmagfm(i)) * (1.-sigmagfm(i)) scaldfunc(i) = max(min(scaldfunc(i), 1.0), 0.) +#ifdef HWRF_SCALESAS + sigmuout(i)=sigmagfm(i) +#endif else scaldfunc(i) = 1.0 endif @@ -2404,6 +2534,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo +#ifndef HWRF_SCALESAS !> - If stochastic physics using cellular automata is .true. then perturb the mass-flux here: if(do_ca)then @@ -2420,6 +2551,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & & edto, xlamd, xmb, c0t, eta, etad, zi, xlamue, xlamud, delp, & qtr, qaero) +#endif c c restore to,qo,uo,vo to t1,q1,u1,v1 in case convection stops c @@ -2437,6 +2569,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo +#ifndef HWRF_SCALESAS do n = 1, ntr do k = 1, km do i = 1, im @@ -2446,6 +2579,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & enddo enddo enddo +#endif c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c c--- feedback: simply the changes from the cloud with unit mass flux @@ -2464,11 +2598,13 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & delvbar(i) = 0. qcond(i) = 0. enddo +#ifndef HWRF_SCALESAS do n = 1, ntr do i = 1, im delebar(i,n) = 0. enddo enddo +#endif do k = 1, km do i = 1, im if (cnvflg(i) .and. k <= kmax(i)) then @@ -2491,6 +2627,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo +#ifndef HWRF_SCALESAS do n = 1, ntr kk = n+2 do k = 1, km @@ -2505,6 +2642,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & enddo enddo enddo +#endif !> - Recalculate saturation specific humidity using the updated temperature. do k = 1, km do i = 1, im @@ -2689,6 +2827,8 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo + +#ifndef HWRF_SCALESAS do n = 1, ntr kk = n+2 do k = 1, km @@ -2716,6 +2856,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & enddo enddo endif +#endif ! ! hchuang code change ! @@ -2751,6 +2892,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & ! ! include TKE contribution from deep convection ! +#ifndef HWRF_SCALESAS if (ntk > 0) then ! do k = 2, km1 @@ -2798,6 +2940,8 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & enddo enddo endif + +#endif return end subroutine samfdeepcnv_run From bff254722b7706b7fb50b13f69e31fb125f5f4cc Mon Sep 17 00:00:00 2001 From: Man Zhang Date: Fri, 13 Dec 2019 21:17:18 -0700 Subject: [PATCH 02/14] add preprocessor directives for HWRF in samfshalcnv --- physics/samfshalcnv.f | 159 +++++++++++++++++++++++++++++++++++++----- 1 file changed, 142 insertions(+), 17 deletions(-) diff --git a/physics/samfshalcnv.f b/physics/samfshalcnv.f index ed80a2f54..ae212c98e 100644 --- a/physics/samfshalcnv.f +++ b/physics/samfshalcnv.f @@ -167,8 +167,11 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & parameter(dtke=tkemx-tkemn) parameter(dthk=25.) parameter(cinpcrmx=180.,cinpcrmn=120.) -! parameter(cinacrmx=-120.,cinacrmn=-120.) +#ifdef HWRF_SCALESAS + parameter(cinacrmx=-120.,cinacrmn=-120.) +#else parameter(cinacrmx=-120.,cinacrmn=-80.) +#endif parameter(crtlamd=3.e-4) parameter(dtmax=10800.,dtmin=600.) parameter(bet1=1.875,cd1=.506,f1=2.0,gam1=.5) @@ -202,7 +205,44 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & ! real(kind=kind_phys) tf, tcr, tcrf parameter (tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf)) -! + +#if HWRF==1 + real*8 :: gasdev,ran1 !zhang + real :: rr !zhang + logical,save :: pert_sas_local !zhang + integer,save :: ens_random_seed_local,env_pp_local !zhang + integer :: ensda_physics_pert !zhang + real,save :: ens_sasamp_local !zhang + data ens_random_seed_local/0/ + data env_pp_local/0/ + CHARACTER(len=3) :: env_memb,env_pp + if ( ens_random_seed_local .eq. 0 ) then + CALL nl_get_ensda_physics_pert(1,ensda_physics_pert) + ens_random_seed_local=ens_random_seed + env_pp_local=ensda_physics_pert + pert_sas_local=.false. + ens_sasamp_local=0.0 +! env_pp=1: do physics perturbations for ensda members, ens_random_seed must be 99 + if ( env_pp_local .eq. 1 ) then + if ( ens_random_seed .ne. 99 ) then + pert_sas_local=.true. + ens_sasamp_local=ens_sasamp + else +! ens_random_seed=99 do physics perturbation for ensemble forecasts, env_pp must be zero + ens_random_seed_local=ens_random_seed + pert_sas_local=pert_sas + ens_sasamp_local=ens_sasamp + endif + else + ens_random_seed_local=ens_random_seed + pert_sas_local=pert_sas + ens_sasamp_local=ens_sasamp + endif + + print*, "SHSAS ==", ens_random_seed_local,pert_sas_local,ens_sasamp_local,ensda_physics_pert + endif +#endif + c----------------------------------------------------------------------- ! ! Initialize CCPP error handling variables @@ -216,9 +256,11 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & fact2 = hvap/rv-fact1*t0c c----------------------------------------------------------------------- +#ifndef HWRF_SCALESAS !> ## Determine whether to perform aerosol transport do_aerosols = (itc > 0) .and. (ntc > 0) .and. (ntr > 0) if (do_aerosols) do_aerosols = (ntr >= itc + ntc - 3) +#endif ! !************************************************************************ ! convert input Pa terms to Cb terms -- Moorthi @@ -253,6 +295,11 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & cina(i) = 0. vshear(i) = 0. gdx(i) = sqrt(garea(i)) + +#ifdef HWRF_SCALESAS + scaldfunc(i)=-1.0 ! wang initialized + sigmagfm(i)=-1.0 +#endif enddo !! !> - Return to the calling routine if deep convection is present or the surface buoyancy flux is negative. @@ -265,7 +312,11 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & !> - determine aerosol-aware rain conversion parameter over land do i=1,im if(islimsk(i) == 1) then +#ifdef HWRF_SCALESAS + c0(i) = c0l +#else c0(i) = c0s*asolfac +#endif else c0(i) = c0s endif @@ -303,9 +354,13 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & dt2 = delt ! c model tunable parameters are all here -! clam = .3 -! aafac = .1 +#ifdef HWRF_SCALESAS + clam = .3 + aafac = .1 + pgcon = 0.55 +#else aafac = .05 +#endif c evef = 0.07 evfact = 0.3 evfactl = 0.3 @@ -354,8 +409,16 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & do k = 1, km1 do i=1,im zi(i,k) = 0.5*(zo(i,k)+zo(i,k+1)) +#ifdef HWRF_SCALESAS + xlamue(i,k) = clam / zi(i,k) +#endif enddo enddo +#ifdef HWRF_SCALESAS + do i=1,im + xlamue(i,km) = xlamue(i,km1) + enddo +#endif c c pbl height c @@ -410,6 +473,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & ! ! initialize tracer variables ! +#ifndef HWRF_SCALESAS do n = 3, ntr+2 kk = n-2 do k = 1, km @@ -422,6 +486,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & enddo enddo enddo +#endif !> - Calculate saturation specific humidity and enforce minimum moisture values. do k = 1, km do i=1,im @@ -517,6 +582,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo +#ifndef HWRF_SCALESAS do n = 1, ntr do k = 1, km1 do i=1,im @@ -526,6 +592,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & enddo enddo enddo +#endif c c look for the level of free convection as cloud base c @@ -597,6 +664,18 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & ptem1= .5*(cinpcrmx-cinpcrmn) cinpcr = cinpcrmx - ptem * ptem1 tem1 = pfld(i,kb(i)) - pfld(i,kbcon(i)) +#if HWRF==1 +! randomly perturb the convection trigger +!zzz if( pert_sas_local .and. ens_random_seed_local .gt. 0 ) then + if( pert_sas_local ) then +!zz print*, "zhang inde ens_random_seed=", ens_random_seed,ens_random_seed_local + ens_random_seed_local=ran1(-ens_random_seed_local)*1000 + rr=2.0*ens_sasamp_local*ran1(-ens_random_seed_local)-ens_sasamp_local +!zz print*, "zhang inde shsas=a", cinpcr,ens_sasamp_local,ens_random_seed_local,cinpcr + cinpcr=cinpcr+rr +!zz print*, "zhang inde shsas=b", cinpcr,ens_sasamp_local,ens_random_seed_local,cinpcr + endif +#endif if(tem1 > cinpcr) then cnvflg(i) = .false. endif @@ -612,14 +691,27 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & ! turbulent entrainment rate assumed to be proportional ! to subcloud mean TKE ! - if(ntk > 0) then ! + +#ifdef HWRF_SCALESAS +!c +!c specify the detrainment rate for the updrafts +!c + do i = 1, im + if(cnvflg(i)) then + xlamud(i) = xlamue(i,kbcon(i)) +! xlamud(i) = crtlamd + endif + enddo +#else + if(ntk > 0) then do i= 1, im if(cnvflg(i)) then sumx(i) = 0. tkemean(i) = 0. endif enddo + do k = 1, km1 do i = 1, im if(cnvflg(i)) then @@ -687,6 +779,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & xlamud(i) = 0.001 * clamt(i) endif enddo +#endif c c determine updraft mass flux for the subcloud layers c @@ -742,6 +835,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo ! for tracers +#ifndef HWRF_SCALESAS do n = 1, ntr do i = 1, im if(cnvflg(i)) then @@ -750,6 +844,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo +#endif c ! cm is an enhancement factor in entrainment rates for momentum ! @@ -778,6 +873,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo +#ifndef HWRF_SCALESAS do n = 1, ntr do k = 2, km1 do i = 1, im @@ -793,6 +889,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & enddo enddo enddo +#endif c c taking account into convection inhibition due to existence of c dry layers below cloud base @@ -862,6 +959,9 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & do i = 1, im if(cnvflg(i)) then ! +#ifdef HWRF_SCALESAS + cinacr = cinacrmx +#else if(islimsk(i) == 1) then w1 = w1l w2 = w2l @@ -890,6 +990,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & cinacr = cinacrmx - tem * tem1 ! ! cinacr = cinacrmx +#endif if(cina(i) < cinacr) cnvflg(i) = .false. endif enddo @@ -929,7 +1030,11 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & ! k = kbcon(i) dp = 1000. * del(i,k) +#ifdef HWRF_SCALESAS + xmbmax(i) = dp / (g * dt2) +#else xmbmax(i) = dp / (2. * grav * dt2) +#endif ! ! xmbmax(i) = dp / (grav * dt2) ! @@ -1169,18 +1274,20 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & bb1 = 4.0 bb2 = 0.8 ! -! do i = 1, im -! if (cnvflg(i)) then -! k = kbcon1(i) -! tem = po(i,k) / (rd * to(i,k)) -! wucb = -0.01 * dot(i,k) / (tem * grav) -! if(wucb > 0.) then -! wu2(i,k) = wucb * wucb -! else -! wu2(i,k) = 0. -! endif -! endif -! enddo +#ifdef HWRF_SCALESAS + do i = 1, im + if (cnvflg(i)) then + k = kbcon1(i) + tem = po(i,k) / (rd * to(i,k)) + wucb = -0.01 * dot(i,k) / (tem * grav) + if(wucb > 0.) then + wu2(i,k) = wucb * wucb + else + wu2(i,k) = 0. + endif + endif + enddo +#endif do k = 2, km1 do i = 1, im if (cnvflg(i)) then @@ -1314,6 +1421,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo +#ifndef HWRF_SCALESAS do n = 1, ntr do k = 1, km do i = 1, im @@ -1323,6 +1431,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & enddo enddo enddo +#endif c c--- changed due to subsidence and entrainment c @@ -1367,6 +1476,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo +#ifndef HWRF_SCALESAS do n = 1, ntr do k = 2, km1 do i = 1, im @@ -1383,6 +1493,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & enddo enddo enddo +#endif c c------- cloud top c @@ -1407,6 +1518,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & & qlko_ktcon(i) * grav / dp endif enddo +#ifndef HWRF_SCALESAS do n = 1, ntr do i = 1, im if(cnvflg(i)) then @@ -1417,6 +1529,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo +#endif ! ! compute convective turn-over time ! @@ -1425,8 +1538,10 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & if(cnvflg(i)) then tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i)) dtconv(i) = tem / wc(i) +#ifndef HWRF_SCALESAS tfac = 1. + gdx(i) / 75000. dtconv(i) = tfac * dtconv(i) +#endif dtconv(i) = max(dtconv(i),dtmin) dtconv(i) = max(dtconv(i),dt2) dtconv(i) = min(dtconv(i),dtmax) @@ -1501,6 +1616,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo ! +#ifndef HWRF_SCALESAS !> - Transport aerosols if present ! if (do_aerosols) @@ -1510,6 +1626,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & ! & edto, xlamd, xmb, c0t, eta, etad, zi, xlamue, xlamud, delp, & xmb, c0t, eta, zi, xlamue, xlamud, delp, & qtr, qaero) +#endif ! !> ## For the "feedback control", calculate updated values of the state variables by multiplying the cloud base mass flux and the tendencies calculated per unit cloud base mass flux from the static control. !! - Recalculate saturation specific humidity. @@ -1539,11 +1656,13 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & delvbar(i) = 0. qcond(i) = 0. enddo +#ifndef HWRF_SCALESAS do n = 1, ntr do i = 1, im delebar(i,n) = 0. enddo enddo +#endif do k = 1, km do i = 1, im if (cnvflg(i)) then @@ -1566,6 +1685,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo +#ifndef HWRF_SCALESAS do n = 1, ntr kk = n+2 do k = 1, km @@ -1580,6 +1700,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & enddo enddo enddo +#endif ! !> - Recalculate saturation specific humidity using the updated temperature. do k = 1, km @@ -1750,6 +1871,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & ! endif !> - Store aerosol concentrations if present +#ifndef HWRF_SCALESAS if (do_aerosols) then do n = 1, ntc kk = n + itc - 1 @@ -1762,6 +1884,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & enddo enddo endif +#endif ! ! hchuang code change ! @@ -1787,6 +1910,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & ! ! include TKE contribution from shallow convection ! +#ifndef HWRF_SCALESAS if (ntk > 0) then ! do k = 2, km1 @@ -1804,6 +1928,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & enddo ! endif +#endif !! return end subroutine samfshalcnv_run From a4ac85250abdb82297a8fe3a0034d6c17cc84fbe Mon Sep 17 00:00:00 2001 From: "Man.Zhang" Date: Tue, 17 Dec 2019 17:15:39 -0700 Subject: [PATCH 03/14] fix bugs to pass compilation --- physics/samfdeepcnv.f | 282 +++++++++++++++++---------------------- physics/samfdeepcnv.meta | 8 ++ physics/samfshalcnv.f | 213 ++++++++++++++--------------- physics/samfshalcnv.meta | 8 ++ 4 files changed, 239 insertions(+), 272 deletions(-) diff --git a/physics/samfdeepcnv.f b/physics/samfdeepcnv.f index abd1700c9..49dce2ae9 100644 --- a/physics/samfdeepcnv.f +++ b/physics/samfdeepcnv.f @@ -71,7 +71,7 @@ end subroutine samfdeepcnv_finalize subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & & eps,epsm1,fv,grav,hvap,rd,rv, & & t0c,delt,ntk,ntr,delp, & - & prslp,psp,phil,qtr,q1,t1,u1,v1,fscav, & + & prslp,psp,phil,qtr,q1,t1,u1,v1,fscav,hwrf_samfdeep, & & do_ca,ca_deep,cldwrk,rn,kbot,ktop,kcnv,islimsk,garea, & & dot,ncloud,ud_mf,dd_mf,dt_mf,cnvw,cnvc, & & QLCN, QICN, w_upi, cf_upi, CNV_MFD, & @@ -93,7 +93,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & & prslp(ix,km), garea(im), dot(ix,km), phil(ix,km) real(kind=kind_phys), dimension(:), intent(in) :: fscav real(kind=kind_phys), intent(in) :: ca_deep(ix) - logical, intent(in) :: do_ca + logical, intent(in) :: do_ca, hwrf_samfdeep integer, intent(inout) :: kcnv(im) ! DH* TODO - check dimensions of qtr, ntr+2 correct? *DH @@ -115,7 +115,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & ! *GJF integer :: mp_phys, mp_phys_mg - real(kind=kind_phys), intent(in) :: clam, c0s, c1, & + real(kind=kind_phys), intent(in) :: clam, c0s, c1, & & betal, betas, asolfac, & & evfact, evfactl, pgcon character(len=*), intent(out) :: errmsg @@ -128,8 +128,10 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & real(kind=kind_phys) clamd, tkemx, tkemn, dtke, & beta, dbeta, betamx, betamn, & cxlame, cxlamd, + & cxlamu, & xlamde, xlamdd, - & crtlame, crtlamd + & crtlamu, crtlamd, + & crtlame, c0l ! ! real(kind=kind_phys) detad real(kind=kind_phys) adw, aup, aafac, d0, @@ -180,8 +182,8 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & & rntot(im), vshear(im), xaa0(im), & xlamd(im), xk(im), cina(im), & xmb(im), xmbmax(im), xpwav(im), -! & xpwev(im), xlamx(im), delebar(im,ntr), - & xpwev(im), delebar(im,ntr), + & xpwev(im), xlamx(im), delebar(im,ntr), +! & xpwev(im), delebar(im,ntr), & delubar(im), delvbar(im) ! real(kind=kind_phys) c0(im) @@ -192,17 +194,17 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & ! ! parameters for updraft velocity calculation real(kind=kind_phys) bet1, cd1, f1, gam1, - & bb1, bb2 -! & bb1, bb2, wucb +! & bb1, bb2 + & bb1, bb2, wucb ! c physical parameters - parameter(grav=grav,asolfac=0.89) !HWRF +! parameter(asolfac=0.89) !HWRF ! parameter(grav=grav) ! parameter(elocp=hvap/cp,el2orc=hvap*hvap/(rv*cp)) ! parameter(c0s=.002,c1=.002,d0=.01) ! parameter(d0=.01) parameter(d0=.001) - parameter(c0l=c0s*asolfac) +!mz parameter(c0l=c0s*asolfac) ! ! asolfac: aerosol-aware parameter based on Lim (2011) ! asolfac= cx / c0s(=.002) @@ -232,6 +234,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & ! for updraft velocity calculation real(kind=kind_phys) wu2(im,km), buo(im,km), drag(im,km) real(kind=kind_phys) wc(im), scaldfunc(im), sigmagfm(im) + real(kind=kind_phys) sigmuout(im) ! c cloud water ! real(kind=kind_phys) tvo(im,km) @@ -310,21 +313,18 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & endif #endif + if(.not. hwrf_samfdeep) then + elocp = hvap/cp + el2orc = hvap*hvap/(rv*cp) -#endif - -#ifndef HWRF_SCALESAS - elocp = hvap/cp - el2orc = hvap*hvap/(rv*cp) - - fact1 = (cvap-cliq)/rv - fact2 = hvap/rv-fact1*t0c + fact1 = (cvap-cliq)/rv + fact2 = hvap/rv-fact1*t0c ! c----------------------------------------------------------------------- !> ## Determine whether to perform aerosol transport - do_aerosols = (itc > 0) .and. (ntc > 0) .and. (ntr > 0) - if (do_aerosols) do_aerosols = (ntr >= itc + ntc - 3) -#endif + do_aerosols = (itc > 0) .and. (ntc > 0) .and. (ntr > 0) + if (do_aerosols) do_aerosols = (ntr >= itc + ntc - 3) + endif ! c----------------------------------------------------------------------- !> ## Compute preliminary quantities needed for static, dynamic, and feedback control portions of the algorithm. @@ -373,21 +373,22 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & vshear(i) = 0. gdx(i) = sqrt(garea(i)) -#ifdef HWRF_SCALESAS + if( hwrf_samfdeep ) then scaldfunc(i)=-1.0 ! initialized wang sigmagfm(i)=-1.0 sigmuout(i)=-1.0 -#endif + endif enddo ! + c0l=c0s*asolfac !> - determine aerosol-aware rain conversion parameter over land do i=1,im if(islimsk(i) == 1) then -#ifdef HWRF_SCALESAS + if (hwrf_samfdeep) then c0(i) = c0l -#else + else c0(i) = c0s*asolfac -#endif + endif else c0(i) = c0s endif @@ -453,32 +454,15 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & ! model tunable parameters are all here edtmaxl = .3 edtmaxs = .3 -#ifdef HWRF_SCALESAS - clam = .1 - aafac = .1 - betal = .05 - betas = .05 - evfact = 0.3 - evfactl = 0.3 -#else -! clam = .1 - aafac = .05 -! betal = .15 -! betas = .15 -! betal = .05 -! betas = .05 -! evef = 0.07 -! evfact = 0.3 -! evfactl = 0.3 -#endif -! -#ifdef HWRF_SCALESAS - crtlamu = 1.0e-4 - cxlamu = 1.0e-3 -#else - crtlame = 1.0e-4 - cxlame = 1.0e-4 -#endif + if (hwrf_samfdeep) then + aafac = .1 + crtlamu = 1.0e-4 + cxlamu = 1.0e-3 + else + aafac = .05 + crtlame = 1.0e-4 + cxlame = 1.0e-4 + endif crtlamd = 1.0e-4 cxlamd = 1.0e-4 @@ -534,9 +518,9 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & do k = 1, km1 do i=1,im zi(i,k) = 0.5*(zo(i,k)+zo(i,k+1)) -#ifdef HWRF_SCALESAS - xlamue(i,k) = clam / zi(i,k) -#endif + if (hwrf_samfdeep) then + xlamue(i,k) = clam / zi(i,k) + endif enddo enddo c @@ -584,7 +568,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & ! ! initialize tracer variables ! -#ifndef HWRF_SCALESAS + if(.not.hwrf_samfdeep) then do n = 3, ntr+2 kk = n-2 do k = 1, km @@ -598,7 +582,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & enddo enddo enddo -#endif + endif ! !> - Calculate saturation specific humidity and enforce minimum moisture values. do k = 1, km @@ -696,7 +680,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & enddo enddo -#ifndef HWRF_SCALESAS + if (.not.hwrf_samfdeep) then do n = 1, ntr do k = 1, km1 do i=1,im @@ -706,7 +690,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & enddo enddo enddo -#endif + endif c c look for the level of free convection as cloud base c @@ -800,9 +784,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & enddo if(totflg) return -#ifndef HWRF_SCALESAS -!! -! + if (.not. hwrf_samfdeep) then ! turbulent entrainment rate assumed to be proportional ! to subcloud mean TKE ! @@ -864,7 +846,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & enddo enddo -#endif + endif !(.not.hwrf_samfdeep) c c assume that updraft entrainment rate above cloud base is c same as that at cloud base @@ -874,7 +856,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & !! \epsilon = \epsilon_0F_0 + d_1\left(1-RH\right)F_1 !! \f] !! where \f$\epsilon_0\f$ is the cloud base entrainment rate, \f$d_1\f$ is a tunable constant, and \f$F_0=\left(\frac{q_s}{q_{s,b}}\right)^2\f$ and \f$F_1=\left(\frac{q_s}{q_{s,b}}\right)^3\f$ where \f$q_s\f$ and \f$q_{s,b}\f$ are the saturation specific humidities at a given level and cloud base, respectively. The detrainment rate in the cloud is assumed to be equal to the entrainment rate at cloud base. -#ifdef HWRF_SCALESAS + if (hwrf_samfdeep) then do i=1,im if(cnvflg(i)) then xlamx(i) = xlamue(i,kbcon(i)) @@ -887,8 +869,8 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & xlamue(i,k) = xlamx(i) endif enddo - enddo -#endif + enddo + endif c c specify detrainment rate for the updrafts c @@ -898,11 +880,11 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & do k = 1, km1 do i=1,im if(cnvflg(i) .and. k < kmax(i)) then -#ifdef HWRF_SCALESAS - xlamud(i,k) = xlamx(i) -#else - xlamud(i,k) = 0.001 * clamt(i) -#endif + if (hwrf_samfdeep) then + xlamud(i,k) = xlamx(i) + else + xlamud(i,k) = 0.001 * clamt(i) + endif endif enddo enddo @@ -932,10 +914,10 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & & (k > kbcon(i) .and. k < kmax(i))) then tem = cxlame * frh(i,k) * fent2(i,k) xlamue(i,k) = xlamue(i,k)*fent1(i,k) + tem -#ifndef HWRF_SCALESAS - tem1 = cxlamd * frh(i,k) - xlamud(i,k) = xlamud(i,k) + tem1 -#endif + if (.not.hwrf_samfdeep) then + tem1 = cxlamd * frh(i,k) + xlamud(i,k) = xlamud(i,k) + tem1 + endif endif enddo enddo @@ -997,17 +979,17 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & pwavo(i) = 0. endif enddo -#ifndef HWRF_SCALESAS + if (.not.hwrf_samfdeep) then ! for tracers - do n = 1, ntr + do n = 1, ntr do i = 1, im if(cnvflg(i)) then indx = kb(i) ecko(i,indx,n) = ctro(i,indx,n) endif enddo - enddo -#endif + enddo + endif c c cloud property is modified by the entrainment process c @@ -1038,9 +1020,9 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo -#ifndef HWRF_SCALESAS - do n = 1, ntr - do k = 2, km1 + if (.not.hwrf_samfdeep) then + do n = 1, ntr + do k = 2, km1 do i = 1, im if (cnvflg(i)) then if(k > kb(i) .and. k < kmax(i)) then @@ -1052,9 +1034,9 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & endif endif enddo - enddo - enddo -#endif + enddo + enddo + endif c c taking account into convection inhibition due to existence of c dry layers below cloud base @@ -1124,7 +1106,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & do i = 1, im if(cnvflg(i)) then ! -#ifndef HWRF_SCALESAS + if (.not.hwrf_samfdeep) then if(islimsk(i) == 1) then w1 = w1l w2 = w2l @@ -1151,11 +1133,9 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & tem = 1. - tem tem1= .5*(cinacrmx-cinacrmn) cinacr = cinacrmx - tem * tem1 -! -! cinacr = cinacrmx -#else + else cinacr = cinacrmx -#endif + endif if(cina(i) < cinacr) cnvflg(i) = .false. endif enddo @@ -1238,15 +1218,13 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & !> - Calculate the maximum value of the cloud base mass flux using the CFL-criterion-based formula of Han and Pan (2011) \cite han_and_pan_2011, equation 7. do i = 1, im if(cnvflg(i)) then -! xmbmax(i) = .1 -! k = kbcon(i) dp = 1000. * del(i,k) -#ifndef HWRF_SCALEASA - xmbmax(i) = dp / (2. * grav * dt2) -#else - xmbmax(i) = dp / (grav * dt2) -#endif + if (.not.hwrf_samfdeep) then + xmbmax(i) = dp / (2. * grav * dt2) + else + xmbmax(i) = dp / (grav * dt2) + endif endif enddo c @@ -1286,13 +1264,8 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & c if(k >= kbcon(i) .and. dq > 0.) then etah = .5 * (eta(i,k) + eta(i,k-1)) -#ifndef HWRF_SCALESAS dp = 1000. * del(i,k) -#endif if(ncloud > 0 .and. k > jmin(i)) then -#ifdef HWRF_SCALESAS - dp = 1000. * del(i,k) -#endif ptem = c0t(i,k) + c1 qlk = dq / (eta(i,k) + etah * ptem * dz) dellal(i,k) = etah * c1 * dz * qlk * grav / dp @@ -1464,13 +1437,8 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & c if(dq > 0.) then etah = .5 * (eta(i,k) + eta(i,k-1)) -#ifndef HWRF_SCALESAS dp = 1000. * del(i,k) -#endif if(ncloud > 0) then -#ifdef HWRF_SCALESAS - dp = 1000. * del(i,k) -#endif ptem = c0t(i,k) + c1 qlk = dq / (eta(i,k) + etah * ptem * dz) dellal(i,k) = etah * c1 * dz * qlk * grav / dp @@ -1493,12 +1461,12 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & ! bb1 = 4.0 bb2 = 0.8 -#ifdef HWRF_SCALESAS + if (hwrf_samfdeep) then do i = 1, im if (cnvflg(i)) then k = kbcon1(i) tem = po(i,k) / (rd * to(i,k)) - wucb = -0.01 * dot(i,k) / (tem * g) + wucb = -0.01 * dot(i,k) / (tem * grav) if(wucb.gt.0.) then wu2(i,k) = wucb * wucb else @@ -1506,7 +1474,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & endif endif enddo -#endif + endif ! do k = 2, km1 do i = 1, im @@ -1659,10 +1627,10 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & enddo enddo do i = 1, im -#ifdef HWRF_SCALESAS + if (hwrf_samfdeep) then beta = betas if(islimsk(i) == 1) beta = betal -#else + else betamn = betas if(islimsk(i) == 1) betamn = betal if(ntk > 0) then @@ -1678,7 +1646,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & else beta = betamn endif -#endif + endif if(cnvflg(i)) then dz = (sumx(i)+zi(i,1))/float(kbcon(i)) tem = 1./float(kbcon(i)) @@ -1720,7 +1688,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo ! for tracers -#ifndef HWRF_SCALESAS + if (.not.hwrf_samfdeep) then do n = 1, ntr do i = 1, im if(cnvflg(i)) then @@ -1729,7 +1697,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo -#endif + endif cj !> - Calculate the cloud properties as a parcel descends, modified by entrainment and detrainment. Discretization follows Appendix B of Grell (1993) \cite grell_1993 . do k = km1, 1, -1 @@ -1759,7 +1727,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo -#ifndef HWRF_SCALESAS + if(.not.hwrf_samfdeep) then do n = 1, ntr do k = km1, 1, -1 do i = 1, im @@ -1773,7 +1741,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & enddo enddo enddo -#endif + endif c !> - Compute the amount of moisture that is necessary to keep the downdraft saturated. do k = km1, 1, -1 @@ -1876,7 +1844,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo -#ifndef HWRF_SCALESAS + if (.not.hwrf_samfdeep) then do n = 1, ntr do k = 1, km do i = 1, im @@ -1886,7 +1854,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & enddo enddo enddo -#endif + endif do i = 1, im if(cnvflg(i)) then dp = 1000. * del(i,1) @@ -1901,7 +1869,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo -#ifndef HWRF_SCALESAS + if (.not.hwrf_samfdeep) then do n = 1, ntr do i = 1, im if(cnvflg(i)) then @@ -1911,7 +1879,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo -#endif + endif c c--- changed due to subsidence and entrainment c @@ -1976,7 +1944,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo -#ifndef HWRF_SCALESAS + if (.not.hwrf_samfdeep) then do n = 1, ntr do k = 2, km1 do i = 1, im @@ -1998,7 +1966,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & enddo enddo enddo -#endif + endif c c------- cloud top c @@ -2024,7 +1992,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo -#ifndef HWRF_SCALESAS + if (.not.hwrf_samfdeep) then do n = 1, ntr do i = 1, im if(cnvflg(i)) then @@ -2035,7 +2003,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo -#endif + endif c c------- final changed variable per unit mass flux c @@ -2407,10 +2375,10 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & if(cnvflg(i)) then tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i)) dtconv(i) = tem / wc(i) -#ifndef HWRF_SCALESAS - tfac = 1. + gdx(i) / 75000. - dtconv(i) = tfac * dtconv(i) -#endif + if (.not.hwrf_samfdeep) then + tfac = 1. + gdx(i) / 75000. + dtconv(i) = tfac * dtconv(i) + endif dtconv(i) = max(dtconv(i),dtmin) dtconv(i) = min(dtconv(i),dtmax) endif @@ -2523,9 +2491,9 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & if (gdx(i) < dxcrtuf) then scaldfunc(i) = (1.-sigmagfm(i)) * (1.-sigmagfm(i)) scaldfunc(i) = max(min(scaldfunc(i), 1.0), 0.) -#ifdef HWRF_SCALESAS - sigmuout(i)=sigmagfm(i) -#endif + if (hwrf_samfdeep) then + sigmuout(i)=sigmagfm(i) + endif else scaldfunc(i) = 1.0 endif @@ -2534,7 +2502,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo -#ifndef HWRF_SCALESAS + if (.not.hwrf_samfdeep) then !> - If stochastic physics using cellular automata is .true. then perturb the mass-flux here: if(do_ca)then @@ -2551,7 +2519,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & & edto, xlamd, xmb, c0t, eta, etad, zi, xlamue, xlamud, delp, & qtr, qaero) -#endif + endif c c restore to,qo,uo,vo to t1,q1,u1,v1 in case convection stops c @@ -2569,17 +2537,17 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo -#ifndef HWRF_SCALESAS - do n = 1, ntr - do k = 1, km + if (.not.hwrf_samfdeep) then + do n = 1, ntr + do k = 1, km do i = 1, im if (cnvflg(i) .and. k <= kmax(i)) then ctro(i,k,n) = ctr(i,k,n) endif enddo - enddo - enddo -#endif + enddo + enddo + endif c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c c--- feedback: simply the changes from the cloud with unit mass flux @@ -2598,13 +2566,13 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & delvbar(i) = 0. qcond(i) = 0. enddo -#ifndef HWRF_SCALESAS - do n = 1, ntr - do i = 1, im - delebar(i,n) = 0. - enddo - enddo -#endif + if (.not.hwrf_samfdeep) then + do n = 1, ntr + do i = 1, im + delebar(i,n) = 0. + enddo + enddo + endif do k = 1, km do i = 1, im if (cnvflg(i) .and. k <= kmax(i)) then @@ -2627,10 +2595,10 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo -#ifndef HWRF_SCALESAS - do n = 1, ntr + if (.not.hwrf_samfdeep) then + do n = 1, ntr kk = n+2 - do k = 1, km + do k = 1, km do i = 1, im if (cnvflg(i) .and. k <= kmax(i)) then if(k <= ktcon(i)) then @@ -2640,9 +2608,9 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & endif endif enddo - enddo - enddo -#endif + enddo + enddo + endif !> - Recalculate saturation specific humidity using the updated temperature. do k = 1, km do i = 1, im @@ -2828,7 +2796,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & enddo enddo -#ifndef HWRF_SCALESAS + if (.not.hwrf_samfdeep) then do n = 1, ntr kk = n+2 do k = 1, km @@ -2856,7 +2824,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & enddo enddo endif -#endif + endif ! ! hchuang code change ! @@ -2892,7 +2860,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & ! ! include TKE contribution from deep convection ! -#ifndef HWRF_SCALESAS + if (.not.hwrf_samfdeep) then if (ntk > 0) then ! do k = 2, km1 @@ -2941,7 +2909,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & enddo endif -#endif + endif return end subroutine samfdeepcnv_run diff --git a/physics/samfdeepcnv.meta b/physics/samfdeepcnv.meta index 3b54998fc..1fec047a2 100644 --- a/physics/samfdeepcnv.meta +++ b/physics/samfdeepcnv.meta @@ -265,6 +265,14 @@ kind = kind_phys intent = in optional = F +[hwrf_samfdeep] + standard_name = flag_for_hwrf_samfdeepcnv_scheme + long_name = flag for hwrf samfdeepcnv scheme + units = flag + dimensions = () + type = logical + intent = in + optional = F [do_ca] standard_name = flag_for_cellular_automata long_name = cellular automata main switch diff --git a/physics/samfshalcnv.f b/physics/samfshalcnv.f index ae212c98e..65f19919f 100644 --- a/physics/samfshalcnv.f +++ b/physics/samfshalcnv.f @@ -55,7 +55,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & & prslp,psp,phil,qtr,q1,t1,u1,v1,fscav, & & rn,kbot,ktop,kcnv,islimsk,garea, & & dot,ncloud,hpbl,ud_mf,dt_mf,cnvw,cnvc, & - & clam,c0s,c1,pgcon,asolfac,errmsg,errflg) + & clam,c0s,c1,pgcon,asolfac,hwrf_samfshal,errmsg,errflg) ! use machine , only : kind_phys use funcphys , only : fpvs @@ -82,6 +82,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & ! real(kind=kind_phys), intent(in) :: clam, c0s, c1, & & asolfac, pgcon + logical, intent(in) :: hwrf_samfshal character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg ! @@ -140,8 +141,8 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & ! ! parameters for updraft velocity calculation real(kind=kind_phys) bet1, cd1, f1, gam1, - & bb1, bb2 -! & bb1, bb2, wucb +! & bb1, bb2 + & bb1, bb2, wucb cc c physical parameters @@ -167,11 +168,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & parameter(dtke=tkemx-tkemn) parameter(dthk=25.) parameter(cinpcrmx=180.,cinpcrmn=120.) -#ifdef HWRF_SCALESAS - parameter(cinacrmx=-120.,cinacrmn=-120.) -#else - parameter(cinacrmx=-120.,cinacrmn=-80.) -#endif + parameter(cinacrmx=-120.) parameter(crtlamd=3.e-4) parameter(dtmax=10800.,dtmin=600.) parameter(bet1=1.875,cd1=.506,f1=2.0,gam1=.5) @@ -255,12 +252,18 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & fact1 = (cvap-cliq)/rv fact2 = hvap/rv-fact1*t0c + if (hwrf_samfshal) then + cinacrmn=-120. + else + cinacrmn=-80. + endif + c----------------------------------------------------------------------- -#ifndef HWRF_SCALESAS + if (.not.hwrf_samfshal) then !> ## Determine whether to perform aerosol transport - do_aerosols = (itc > 0) .and. (ntc > 0) .and. (ntr > 0) - if (do_aerosols) do_aerosols = (ntr >= itc + ntc - 3) -#endif + do_aerosols = (itc > 0) .and. (ntc > 0) .and. (ntr > 0) + if (do_aerosols) do_aerosols = (ntr >= itc + ntc - 3) + endif ! !************************************************************************ ! convert input Pa terms to Cb terms -- Moorthi @@ -296,10 +299,10 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & vshear(i) = 0. gdx(i) = sqrt(garea(i)) -#ifdef HWRF_SCALESAS + if (hwrf_samfshal) then scaldfunc(i)=-1.0 ! wang initialized sigmagfm(i)=-1.0 -#endif + endif enddo !! !> - Return to the calling routine if deep convection is present or the surface buoyancy flux is negative. @@ -312,11 +315,11 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & !> - determine aerosol-aware rain conversion parameter over land do i=1,im if(islimsk(i) == 1) then -#ifdef HWRF_SCALESAS + if (hwrf_samfshal) then c0(i) = c0l -#else + else c0(i) = c0s*asolfac -#endif + endif else c0(i) = c0s endif @@ -354,19 +357,15 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & dt2 = delt ! c model tunable parameters are all here -#ifdef HWRF_SCALESAS - clam = .3 - aafac = .1 - pgcon = 0.55 -#else - aafac = .05 -#endif + if (hwrf_samfshal) then + aafac = .1 + else + aafac = .05 + endif c evef = 0.07 evfact = 0.3 evfactl = 0.3 ! -! pgcon = 0.7 ! Gregory et al. (1997, QJRMS) -! pgcon = 0.55 ! Zhang & Wu (2003,JAS) w1l = -8.e-3 w2l = -4.e-2 w3l = -5.e-3 @@ -409,16 +408,16 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & do k = 1, km1 do i=1,im zi(i,k) = 0.5*(zo(i,k)+zo(i,k+1)) -#ifdef HWRF_SCALESAS - xlamue(i,k) = clam / zi(i,k) -#endif + if (hwrf_samfshal) then + xlamue(i,k) = clam / zi(i,k) + endif enddo enddo -#ifdef HWRF_SCALESAS - do i=1,im + if (hwrf_samfshal) then + do i=1,im xlamue(i,km) = xlamue(i,km1) - enddo -#endif + enddo + endif c c pbl height c @@ -473,9 +472,9 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & ! ! initialize tracer variables ! -#ifndef HWRF_SCALESAS - do n = 3, ntr+2 - kk = n-2 + if (.not.hwrf_samfshal) then + do n = 3, ntr+2 + kk = n-2 do k = 1, km do i = 1, im if (cnvflg(i) .and. k <= kmax(i)) then @@ -485,8 +484,8 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo - enddo -#endif + enddo + endif !> - Calculate saturation specific humidity and enforce minimum moisture values. do k = 1, km do i=1,im @@ -582,17 +581,17 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo -#ifndef HWRF_SCALESAS - do n = 1, ntr - do k = 1, km1 + if (.not.hwrf_samfshal) then + do n = 1, ntr + do k = 1, km1 do i=1,im if (cnvflg(i) .and. k <= kmax(i)-1) then ctro(i,k,n) = .5 * (ctro(i,k,n) + ctro(i,k+1,n)) endif enddo - enddo - enddo -#endif + enddo + enddo + endif c c look for the level of free convection as cloud base c @@ -693,17 +692,17 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & ! ! -#ifdef HWRF_SCALESAS !c !c specify the detrainment rate for the updrafts !c + if (hwrf_samfshal) then do i = 1, im if(cnvflg(i)) then xlamud(i) = xlamue(i,kbcon(i)) ! xlamud(i) = crtlamd endif enddo -#else + else if(ntk > 0) then do i= 1, im if(cnvflg(i)) then @@ -779,7 +778,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & xlamud(i) = 0.001 * clamt(i) endif enddo -#endif + endif ! hwrf_samfshal c c determine updraft mass flux for the subcloud layers c @@ -835,7 +834,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo ! for tracers -#ifndef HWRF_SCALESAS + if (.not. hwrf_samfshal) then do n = 1, ntr do i = 1, im if(cnvflg(i)) then @@ -844,7 +843,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo -#endif + endif c ! cm is an enhancement factor in entrainment rates for momentum ! @@ -873,9 +872,9 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo -#ifndef HWRF_SCALESAS - do n = 1, ntr - do k = 2, km1 + if (.not.hwrf_samfshal) then + do n = 1, ntr + do k = 2, km1 do i = 1, im if (cnvflg(i)) then if(k > kb(i) .and. k < kmax(i)) then @@ -887,9 +886,9 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & endif endif enddo - enddo - enddo -#endif + enddo + enddo + endif c c taking account into convection inhibition due to existence of c dry layers below cloud base @@ -959,9 +958,9 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & do i = 1, im if(cnvflg(i)) then ! -#ifdef HWRF_SCALESAS + if (hwrf_samfshal) then cinacr = cinacrmx -#else + else if(islimsk(i) == 1) then w1 = w1l w2 = w2l @@ -988,9 +987,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & tem = 1. - tem tem1= .5*(cinacrmx-cinacrmn) cinacr = cinacrmx - tem * tem1 -! -! cinacr = cinacrmx -#endif + endif if(cina(i) < cinacr) cnvflg(i) = .false. endif enddo @@ -1030,16 +1027,11 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & ! k = kbcon(i) dp = 1000. * del(i,k) -#ifdef HWRF_SCALESAS - xmbmax(i) = dp / (g * dt2) -#else - xmbmax(i) = dp / (2. * grav * dt2) -#endif -! -! xmbmax(i) = dp / (grav * dt2) -! -! tem = dp / (grav * dt2) -! xmbmax(i) = min(tem, xmbmax(i)) + if (hwrf_samfshal) then + xmbmax(i) = dp / (grav * dt2) + else + xmbmax(i) = dp / (2. * grav * dt2) + endif endif enddo c @@ -1261,20 +1253,11 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & ! ! compute updraft velocity square(wu2) !> - Calculate updraft velocity square(wu2) according to Han et al.'s (2017) \cite han_et_al_2017 equation 7. -! -! bb1 = 2. * (1.+bet1*cd1) -! bb2 = 2. / (f1*(1.+gam1)) -! -! bb1 = 3.9 -! bb2 = 0.67 -! -! bb1 = 2.0 -! bb2 = 4.0 ! bb1 = 4.0 bb2 = 0.8 ! -#ifdef HWRF_SCALESAS + if (hwrf_samfshal) then do i = 1, im if (cnvflg(i)) then k = kbcon1(i) @@ -1287,7 +1270,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & endif endif enddo -#endif + endif do k = 2, km1 do i = 1, im if (cnvflg(i)) then @@ -1421,17 +1404,17 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo -#ifndef HWRF_SCALESAS - do n = 1, ntr - do k = 1, km + if (.not.hwrf_samfshal) then + do n = 1, ntr + do k = 1, km do i = 1, im if(cnvflg(i) .and. k <= kmax(i)) then dellae(i,k,n) = 0. endif enddo - enddo - enddo -#endif + enddo + enddo + endif c c--- changed due to subsidence and entrainment c @@ -1476,9 +1459,9 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo -#ifndef HWRF_SCALESAS - do n = 1, ntr - do k = 2, km1 + if(.not.hwrf_samfshal) then + do n = 1, ntr + do k = 2, km1 do i = 1, im if (cnvflg(i)) then if(k > kb(i) .and. k < ktcon(i)) then @@ -1491,9 +1474,9 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & endif endif enddo - enddo - enddo -#endif + enddo + enddo + endif c c------- cloud top c @@ -1518,7 +1501,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & & qlko_ktcon(i) * grav / dp endif enddo -#ifndef HWRF_SCALESAS + if (.not.hwrf_samfshal) then do n = 1, ntr do i = 1, im if(cnvflg(i)) then @@ -1529,7 +1512,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo -#endif + endif ! ! compute convective turn-over time ! @@ -1538,10 +1521,10 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & if(cnvflg(i)) then tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i)) dtconv(i) = tem / wc(i) -#ifndef HWRF_SCALESAS - tfac = 1. + gdx(i) / 75000. - dtconv(i) = tfac * dtconv(i) -#endif + if (.not.hwrf_samfshal) then + tfac = 1. + gdx(i) / 75000. + dtconv(i) = tfac * dtconv(i) + endif dtconv(i) = max(dtconv(i),dtmin) dtconv(i) = max(dtconv(i),dt2) dtconv(i) = min(dtconv(i),dtmax) @@ -1616,17 +1599,17 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo ! -#ifndef HWRF_SCALESAS !> - Transport aerosols if present ! - if (do_aerosols) + if (.not.hwrf_samfshal) then + if (do_aerosols) & call samfshalcnv_aerosols(im, ix, km, itc, ntc, ntr, delt, ! & xlamde, xlamdd, cnvflg, jmin, kb, kmax, kbcon, ktcon, fscav, & cnvflg, kb, kmax, kbcon, ktcon, fscav, ! & edto, xlamd, xmb, c0t, eta, etad, zi, xlamue, xlamud, delp, & xmb, c0t, eta, zi, xlamue, xlamud, delp, & qtr, qaero) -#endif + endif ! !> ## For the "feedback control", calculate updated values of the state variables by multiplying the cloud base mass flux and the tendencies calculated per unit cloud base mass flux from the static control. !! - Recalculate saturation specific humidity. @@ -1656,13 +1639,13 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & delvbar(i) = 0. qcond(i) = 0. enddo -#ifndef HWRF_SCALESAS - do n = 1, ntr - do i = 1, im + if (.not. hwrf_samfshal) then + do n = 1, ntr + do i = 1, im delebar(i,n) = 0. - enddo - enddo -#endif + enddo + enddo + endif do k = 1, km do i = 1, im if (cnvflg(i)) then @@ -1685,7 +1668,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo -#ifndef HWRF_SCALESAS + if (.not.hwrf_samfshal) then do n = 1, ntr kk = n+2 do k = 1, km @@ -1700,7 +1683,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & enddo enddo enddo -#endif + endif ! !> - Recalculate saturation specific humidity using the updated temperature. do k = 1, km @@ -1871,8 +1854,8 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & ! endif !> - Store aerosol concentrations if present -#ifndef HWRF_SCALESAS - if (do_aerosols) then + if (.not. hwrf_samfshal) then + if (do_aerosols) then do n = 1, ntc kk = n + itc - 1 do k = 1, km @@ -1884,7 +1867,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & enddo enddo endif -#endif + endif ! ! hchuang code change ! @@ -1910,7 +1893,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & ! ! include TKE contribution from shallow convection ! -#ifndef HWRF_SCALESAS + if (.not.hwrf_samfshal) then if (ntk > 0) then ! do k = 2, km1 @@ -1928,7 +1911,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & enddo ! endif -#endif + endif !! return end subroutine samfshalcnv_run diff --git a/physics/samfshalcnv.meta b/physics/samfshalcnv.meta index 5189afd95..4e7fd3898 100644 --- a/physics/samfshalcnv.meta +++ b/physics/samfshalcnv.meta @@ -422,6 +422,14 @@ kind = kind_phys intent = in optional = F +[hwrf_samfshal] + standard_name = flag_for_hwrf_samfshalcnv_scheme + long_name = flag for hwrf samfshalcnv scheme + units = flag + dimensions = () + type = logical + intent = in + optional = F [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP From beb3a33f128ee12ed567dbc5fe09645f755456e8 Mon Sep 17 00:00:00 2001 From: "Man.Zhang" Date: Wed, 18 Dec 2019 13:17:10 -0700 Subject: [PATCH 04/14] delete HWRF ensemble capability --- physics/samfdeepcnv.f | 52 ++----------------------------------------- physics/samfshalcnv.f | 49 +--------------------------------------- 2 files changed, 3 insertions(+), 98 deletions(-) diff --git a/physics/samfdeepcnv.f b/physics/samfdeepcnv.f index 49dce2ae9..9a38ef453 100644 --- a/physics/samfdeepcnv.f +++ b/physics/samfdeepcnv.f @@ -270,48 +270,11 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & real(kind=kind_phys) tf, tcr, tcrf parameter (tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf)) -#if HWRF==1 - real*8 :: gasdev,ran1 !zhang - real :: rr !zhang - logical,save :: pert_sas_local !zhang - integer,save :: ens_random_seed_local,env_pp_local !zhang - integer :: ensda_physics_pert !zhang - real,save :: ens_sasamp_local !zhang - data ens_random_seed_local/0/ - data env_pp_local/0/ - CHARACTER(len=3) :: env_memb,env_pp -#endif ! Initialize CCPP error handling variables errmsg = '' errflg = 0 -#if HWRF==1 - if ( ens_random_seed_local .eq. 0 ) then - CALL nl_get_ensda_physics_pert(1,ensda_physics_pert) - ens_random_seed_local=ens_random_seed - env_pp_local=ensda_physics_pert - pert_sas_local=.false. - ens_sasamp_local=0.0 -! env_pp=1: do physics perturbations for ensda members, ens_random_seed must be 99 - if ( env_pp_local .eq. 1 ) then - if ( ens_random_seed .ne. 99 ) then - pert_sas_local=.true. - ens_sasamp_local=ens_sasamp - else -! ens_random_seed=99 do physics perturbation for ensemble forecasts, env_pp must be zero - ens_random_seed_local=ens_random_seed - pert_sas_local=pert_sas - ens_sasamp_local=ens_sasamp - endif - else - ens_random_seed_local=ens_random_seed - pert_sas_local=pert_sas - ens_sasamp_local=ens_sasamp - endif - print*, "DESAS ==", ens_random_seed_local,pert_sas_local,ens_sasamp_local,ensda_physics_pert - endif -#endif if(.not. hwrf_samfdeep) then elocp = hvap/cp @@ -374,7 +337,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & gdx(i) = sqrt(garea(i)) if( hwrf_samfdeep ) then - scaldfunc(i)=-1.0 ! initialized wang + scaldfunc(i)=-1.0 sigmagfm(i)=-1.0 sigmuout(i)=-1.0 endif @@ -760,18 +723,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & ptem1= .5*(cinpcrmx-cinpcrmn) cinpcr = cinpcrmx - ptem * ptem1 tem1 = pfld(i,kb(i)) - pfld(i,kbcon(i)) -#if HWRF==1 -! randomly perturb the convection trigger -!zzz if( pert_sas_local .and. ens_random_seed_local .gt. 0 ) then - if( pert_sas_local ) then -!zz print*,"ens_random_seed==",ens_random_seed,ens_random_seed_local - ens_random_seed_local=ran1(-ens_random_seed_local)*1000 - rr=2.0*ens_sasamp_local*ran1(-ens_random_seed_local)-ens_sasamp_local -!zz print*, "zhang inde sas=a", cinpcr,ens_sasamp_local,ens_random_seed_local,cinpcr - cinpcr=cinpcr+rr -!zz print*, "zhang inde sas=b", cinpcr,ens_sasamp_local,ens_random_seed_local,cinpcr - endif -#endif + if(tem1 > cinpcr) then cnvflg(i) = .false. endif diff --git a/physics/samfshalcnv.f b/physics/samfshalcnv.f index 65f19919f..7fa49a856 100644 --- a/physics/samfshalcnv.f +++ b/physics/samfshalcnv.f @@ -203,42 +203,6 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & real(kind=kind_phys) tf, tcr, tcrf parameter (tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf)) -#if HWRF==1 - real*8 :: gasdev,ran1 !zhang - real :: rr !zhang - logical,save :: pert_sas_local !zhang - integer,save :: ens_random_seed_local,env_pp_local !zhang - integer :: ensda_physics_pert !zhang - real,save :: ens_sasamp_local !zhang - data ens_random_seed_local/0/ - data env_pp_local/0/ - CHARACTER(len=3) :: env_memb,env_pp - if ( ens_random_seed_local .eq. 0 ) then - CALL nl_get_ensda_physics_pert(1,ensda_physics_pert) - ens_random_seed_local=ens_random_seed - env_pp_local=ensda_physics_pert - pert_sas_local=.false. - ens_sasamp_local=0.0 -! env_pp=1: do physics perturbations for ensda members, ens_random_seed must be 99 - if ( env_pp_local .eq. 1 ) then - if ( ens_random_seed .ne. 99 ) then - pert_sas_local=.true. - ens_sasamp_local=ens_sasamp - else -! ens_random_seed=99 do physics perturbation for ensemble forecasts, env_pp must be zero - ens_random_seed_local=ens_random_seed - pert_sas_local=pert_sas - ens_sasamp_local=ens_sasamp - endif - else - ens_random_seed_local=ens_random_seed - pert_sas_local=pert_sas - ens_sasamp_local=ens_sasamp - endif - - print*, "SHSAS ==", ens_random_seed_local,pert_sas_local,ens_sasamp_local,ensda_physics_pert - endif -#endif c----------------------------------------------------------------------- ! @@ -663,18 +627,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & ptem1= .5*(cinpcrmx-cinpcrmn) cinpcr = cinpcrmx - ptem * ptem1 tem1 = pfld(i,kb(i)) - pfld(i,kbcon(i)) -#if HWRF==1 -! randomly perturb the convection trigger -!zzz if( pert_sas_local .and. ens_random_seed_local .gt. 0 ) then - if( pert_sas_local ) then -!zz print*, "zhang inde ens_random_seed=", ens_random_seed,ens_random_seed_local - ens_random_seed_local=ran1(-ens_random_seed_local)*1000 - rr=2.0*ens_sasamp_local*ran1(-ens_random_seed_local)-ens_sasamp_local -!zz print*, "zhang inde shsas=a", cinpcr,ens_sasamp_local,ens_random_seed_local,cinpcr - cinpcr=cinpcr+rr -!zz print*, "zhang inde shsas=b", cinpcr,ens_sasamp_local,ens_random_seed_local,cinpcr - endif -#endif + if(tem1 > cinpcr) then cnvflg(i) = .false. endif From c825f5faeb3c658ffcf32b84b0e7527815862812 Mon Sep 17 00:00:00 2001 From: "Man.Zhang" Date: Fri, 20 Dec 2019 14:56:59 -0700 Subject: [PATCH 05/14] remove if outside of loop per Doms suggestion --- physics/samfdeepcnv.f | 193 +++++++++++++++++++++++++++++++----------- physics/samfshalcnv.f | 100 ++++++++++++++-------- 2 files changed, 207 insertions(+), 86 deletions(-) diff --git a/physics/samfdeepcnv.f b/physics/samfdeepcnv.f index 9a38ef453..fcc63c5d1 100644 --- a/physics/samfdeepcnv.f +++ b/physics/samfdeepcnv.f @@ -305,7 +305,41 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & c c initialize arrays c - do i=1,im + if (.not.hwrf_samfdeep) then + do i=1,im + cnvflg(i) = .true. + rn(i)=0. + mbdt(i)=10. + kbot(i)=km+1 + ktop(i)=0 + kbcon(i)=km + ktcon(i)=1 + ktconn(i)=1 + dtconv(i) = 3600. + cldwrk(i) = 0. + pdot(i) = 0. + lmin(i) = 1 + jmin(i) = 1 + qlko_ktcon(i) = 0. + edt(i) = 0. + edto(i) = 0. + edtx(i) = 0. +! acrt(i) = 0. +! acrtfct(i) = 1. + aa1(i) = 0. + aa2(i) = 0. + xaa0(i) = 0. + cina(i) = 0. + pwavo(i)= 0. + pwevo(i)= 0. + xpwav(i)= 0. + xpwev(i)= 0. + vshear(i) = 0. + gdx(i) = sqrt(garea(i)) + enddo + + else + do i=1,im cnvflg(i) = .true. rn(i)=0. mbdt(i)=10. @@ -336,26 +370,22 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & vshear(i) = 0. gdx(i) = sqrt(garea(i)) - if( hwrf_samfdeep ) then - scaldfunc(i)=-1.0 + !mz*HWRF SAS + scaldfunc(i)=-1.0 sigmagfm(i)=-1.0 sigmuout(i)=-1.0 - endif - enddo + enddo + endif ! - c0l=c0s*asolfac !> - determine aerosol-aware rain conversion parameter over land do i=1,im if(islimsk(i) == 1) then - if (hwrf_samfdeep) then - c0(i) = c0l - else c0(i) = c0s*asolfac - endif else c0(i) = c0s endif enddo + !> - determine rain conversion parameter above the freezing level which exponentially decreases with decreasing temperature from Han et al.'s (2017) \cite han_et_al_2017 equation 8. do k = 1, km do i = 1, im @@ -478,14 +508,21 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & enddo enddo !> - Calculate interface height - do k = 1, km1 + if (hwrf_samfdeep) then + do k = 1, km1 do i=1,im zi(i,k) = 0.5*(zo(i,k)+zo(i,k+1)) - if (hwrf_samfdeep) then - xlamue(i,k) = clam / zi(i,k) - endif + xlamue(i,k) = clam / zi(i,k) enddo - enddo + enddo + else + do k = 1, km1 + do i=1,im + zi(i,k) = 0.5*(zo(i,k)+zo(i,k+1)) + enddo + enddo + endif + c c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c convert surface pressure to mb from cb @@ -860,19 +897,29 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & c organized one depending on the environmental relative humidity c (Bechtold et al., 2008; Derbyshire et al., 2011) c - do k = 2, km1 + if (hwrf_samfdeep) then + do k = 2, km1 + do i=1,im + if(cnvflg(i) .and. + & (k > kbcon(i) .and. k < kmax(i))) then + tem = cxlamu * frh(i,k) * fent2(i,k) + xlamue(i,k) = xlamue(i,k)*fent1(i,k) + tem + endif + enddo + enddo + else + do k = 2, km1 do i=1,im if(cnvflg(i) .and. & (k > kbcon(i) .and. k < kmax(i))) then tem = cxlame * frh(i,k) * fent2(i,k) xlamue(i,k) = xlamue(i,k)*fent1(i,k) + tem - if (.not.hwrf_samfdeep) then - tem1 = cxlamd * frh(i,k) - xlamud(i,k) = xlamud(i,k) + tem1 - endif + tem1 = cxlamd * frh(i,k) + xlamud(i,k) = xlamud(i,k) + tem1 endif enddo - enddo + enddo + endif ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c @@ -1055,10 +1102,17 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & enddo enddo !> - Turn off convection if the CIN is less than a critical value (cinacr) which is inversely proportional to the large-scale vertical velocity. - do i = 1, im + + if(hwrf_samfdeep) then + do i = 1, im + if(cnvflg(i)) then + cinacr = cinacrmx + if(cina(i) < cinacr) cnvflg(i) = .false. + endif + enddo + else !gfs_samfdeep + do i = 1, im if(cnvflg(i)) then -! - if (.not.hwrf_samfdeep) then if(islimsk(i) == 1) then w1 = w1l w2 = w2l @@ -1085,12 +1139,10 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & tem = 1. - tem tem1= .5*(cinacrmx-cinacrmn) cinacr = cinacrmx - tem * tem1 - else - cinacr = cinacrmx - endif if(cina(i) < cinacr) cnvflg(i) = .false. endif - enddo + enddo + endif !hwrf_samfdeep !! totflg = .true. do i=1,im @@ -1168,17 +1220,23 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & c specify upper limit of mass flux at cloud base c !> - Calculate the maximum value of the cloud base mass flux using the CFL-criterion-based formula of Han and Pan (2011) \cite han_and_pan_2011, equation 7. - do i = 1, im + if(hwrf_samfdeep) then + do i = 1, im if(cnvflg(i)) then k = kbcon(i) dp = 1000. * del(i,k) - if (.not.hwrf_samfdeep) then - xmbmax(i) = dp / (2. * grav * dt2) - else - xmbmax(i) = dp / (grav * dt2) - endif + xmbmax(i) = dp / (grav * dt2) endif - enddo + enddo + else + do i = 1, im + if(cnvflg(i)) then + k = kbcon(i) + dp = 1000. * del(i,k) + xmbmax(i) = dp / (2. * grav * dt2) + endif + enddo + endif c c compute cloud moisture property and precipitation c @@ -1578,11 +1636,19 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo - do i = 1, im - if (hwrf_samfdeep) then + + if (hwrf_samfdeep) then + do i = 1, im beta = betas if(islimsk(i) == 1) beta = betal - else + if(cnvflg(i)) then + dz = (sumx(i)+zi(i,1))/float(kbcon(i)) + tem = 1./float(kbcon(i)) + xlamd(i) = (1.-beta**tem)/dz + endif + enddo + else + do i = 1, im betamn = betas if(islimsk(i) == 1) betamn = betal if(ntk > 0) then @@ -1598,13 +1664,13 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & else beta = betamn endif - endif if(cnvflg(i)) then dz = (sumx(i)+zi(i,1))/float(kbcon(i)) tem = 1./float(kbcon(i)) xlamd(i) = (1.-beta**tem)/dz endif - enddo + enddo + endif c c determine downdraft mass flux c @@ -2323,18 +2389,29 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & ! compute convective turn-over time ! !> - Following Bechtold et al. (2008) \cite bechtold_et_al_2008, the convective adjustment time (dtconv) is set to be proportional to the convective turnover time, which is computed using the mean updraft velocity (wc) and the cloud depth. It is also proportional to the grid size (gdx). - do i= 1, im + + if(hwrf_samfdeep) then + do i= 1, im if(cnvflg(i)) then tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i)) dtconv(i) = tem / wc(i) - if (.not.hwrf_samfdeep) then - tfac = 1. + gdx(i) / 75000. - dtconv(i) = tfac * dtconv(i) - endif dtconv(i) = max(dtconv(i),dtmin) dtconv(i) = min(dtconv(i),dtmax) endif - enddo + enddo + else + do i= 1, im + if(cnvflg(i)) then + tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i)) + dtconv(i) = tem / wc(i) + tfac = 1. + gdx(i) / 75000. + dtconv(i) = tfac * dtconv(i) + dtconv(i) = max(dtconv(i),dtmin) + dtconv(i) = min(dtconv(i),dtmax) + endif + enddo + endif + ! !> - Calculate advective time scale (tauadv) using a mean cloud layer wind speed. do i= 1, im @@ -2438,21 +2515,35 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & enddo ! !> - Then, calculate the reduction factor (scaldfunc) of the vertical convective eddy transport of mass flux as a function of updraft fraction from the studies by Arakawa and Wu (2013) \cite arakawa_and_wu_2013 (also see Han et al.'s (2017) \cite han_et_al_2017 equation 1 and 2). The final cloud base mass flux with scale-aware parameterization is obtained from the mass flux when sigmagfm << 1, multiplied by the reduction factor (Han et al.'s (2017) \cite han_et_al_2017 equation 2). - do i = 1, im + if(hwrf_samfdeep) then + do i = 1, im if(cnvflg(i)) then if (gdx(i) < dxcrtuf) then scaldfunc(i) = (1.-sigmagfm(i)) * (1.-sigmagfm(i)) scaldfunc(i) = max(min(scaldfunc(i), 1.0), 0.) - if (hwrf_samfdeep) then - sigmuout(i)=sigmagfm(i) - endif + sigmuout(i)=sigmagfm(i) else scaldfunc(i) = 1.0 endif xmb(i) = xmb(i) * scaldfunc(i) xmb(i) = min(xmb(i),xmbmax(i)) endif - enddo + enddo + + else + do i = 1, im + if(cnvflg(i)) then + if (gdx(i) < dxcrtuf) then + scaldfunc(i) = (1.-sigmagfm(i)) * (1.-sigmagfm(i)) + scaldfunc(i) = max(min(scaldfunc(i), 1.0), 0.) + else + scaldfunc(i) = 1.0 + endif + xmb(i) = xmb(i) * scaldfunc(i) + xmb(i) = min(xmb(i),xmbmax(i)) + endif + enddo + endif if (.not.hwrf_samfdeep) then !> - If stochastic physics using cellular automata is .true. then perturb the mass-flux here: diff --git a/physics/samfshalcnv.f b/physics/samfshalcnv.f index 7fa49a856..e21110bd6 100644 --- a/physics/samfshalcnv.f +++ b/physics/samfshalcnv.f @@ -216,9 +216,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & fact1 = (cvap-cliq)/rv fact2 = hvap/rv-fact1*t0c - if (hwrf_samfshal) then - cinacrmn=-120. - else + if (.not.hwrf_samfshal) then cinacrmn=-80. endif @@ -243,7 +241,8 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & c initialize arrays c !> - Initialize column-integrated and other single-value-per-column variable arrays. - do i=1,im + if(hwrf_samfshal) then + do i=1,im cnvflg(i) = .true. if(kcnv(i) == 1) cnvflg(i) = .false. if(cnvflg(i)) then @@ -262,12 +261,32 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & cina(i) = 0. vshear(i) = 0. gdx(i) = sqrt(garea(i)) - - if (hwrf_samfshal) then scaldfunc(i)=-1.0 ! wang initialized sigmagfm(i)=-1.0 + enddo + + else !gfs_samfshal + do i=1,im + cnvflg(i) = .true. + if(kcnv(i) == 1) cnvflg(i) = .false. + if(cnvflg(i)) then + kbot(i)=km+1 + ktop(i)=0 endif - enddo + rn(i)=0. + kbcon(i)=km + ktcon(i)=1 + ktconn(i)=1 + kb(i)=km + pdot(i) = 0. + qlko_ktcon(i) = 0. + edt(i) = 0. + aa1(i) = 0. + cina(i) = 0. + vshear(i) = 0. + gdx(i) = sqrt(garea(i)) + enddo + endif !! !> - Return to the calling routine if deep convection is present or the surface buoyancy flux is negative. totflg = .true. @@ -279,11 +298,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & !> - determine aerosol-aware rain conversion parameter over land do i=1,im if(islimsk(i) == 1) then - if (hwrf_samfshal) then - c0(i) = c0l - else c0(i) = c0s*asolfac - endif else c0(i) = c0s endif @@ -369,18 +384,22 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & enddo enddo !> - Calculate interface height - do k = 1, km1 + if(hwrf_samfshal) then + do k = 1, km1 do i=1,im zi(i,k) = 0.5*(zo(i,k)+zo(i,k+1)) - if (hwrf_samfshal) then - xlamue(i,k) = clam / zi(i,k) - endif + xlamue(i,k) = clam / zi(i,k) enddo - enddo - if (hwrf_samfshal) then + enddo do i=1,im xlamue(i,km) = xlamue(i,km1) enddo + else + do k = 1, km1 + do i=1,im + zi(i,k) = 0.5*(zo(i,k)+zo(i,k+1)) + enddo + enddo endif c c pbl height @@ -545,6 +564,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo + if (.not.hwrf_samfshal) then do n = 1, ntr do k = 1, km1 @@ -649,12 +669,12 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & !c specify the detrainment rate for the updrafts !c if (hwrf_samfshal) then - do i = 1, im + do i = 1, im if(cnvflg(i)) then xlamud(i) = xlamue(i,kbcon(i)) ! xlamud(i) = crtlamd endif - enddo + enddo else if(ntk > 0) then do i= 1, im @@ -825,6 +845,7 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo + if (.not.hwrf_samfshal) then do n = 1, ntr do k = 2, km1 @@ -908,12 +929,17 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & enddo enddo !> - Turn off convection if the CIN is less than a critical value (cinacr) which is inversely proportional to the large-scale vertical velocity. - do i = 1, im + + if (hwrf_samfshal) then + do i = 1, im + if(cnvflg(i)) then + cinacr = cinacrmx + if(cina(i) < cinacr) cnvflg(i) = .false. + endif + enddo + else + do i = 1, im if(cnvflg(i)) then -! - if (hwrf_samfshal) then - cinacr = cinacrmx - else if(islimsk(i) == 1) then w1 = w1l w2 = w2l @@ -942,8 +968,8 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & cinacr = cinacrmx - tem * tem1 endif if(cina(i) < cinacr) cnvflg(i) = .false. - endif - enddo + enddo + endif !! totflg = .true. do i=1,im @@ -974,19 +1000,23 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & c specify upper limit of mass flux at cloud base c !> - Calculate the maximum value of the cloud base mass flux using the CFL-criterion-based formula of Han and Pan (2011) \cite han_and_pan_2011, equation 7. - do i = 1, im + if(hwrf_samfshal) then + do i = 1, im if(cnvflg(i)) then -! xmbmax(i) = .1 -! k = kbcon(i) dp = 1000. * del(i,k) - if (hwrf_samfshal) then - xmbmax(i) = dp / (grav * dt2) - else - xmbmax(i) = dp / (2. * grav * dt2) - endif + xmbmax(i) = dp / (grav * dt2) endif - enddo + enddo + else + do i = 1, im + if(cnvflg(i)) then + k = kbcon(i) + dp = 1000. * del(i,k) + xmbmax(i) = dp / (2. * grav * dt2) + endif + enddo + endif c c compute cloud moisture property and precipitation c From 029f4489d4f06d48e31601912f2cbfe92435c47e Mon Sep 17 00:00:00 2001 From: "Man.Zhang" Date: Thu, 26 Dec 2019 16:15:13 -0700 Subject: [PATCH 06/14] bug fix --- physics/samfdeepcnv.f | 63 ++++++++++++++++++++++--------------------- 1 file changed, 32 insertions(+), 31 deletions(-) diff --git a/physics/samfdeepcnv.f b/physics/samfdeepcnv.f index fcc63c5d1..07b30db51 100644 --- a/physics/samfdeepcnv.f +++ b/physics/samfdeepcnv.f @@ -130,8 +130,8 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & & cxlame, cxlamd, & cxlamu, & xlamde, xlamdd, - & crtlamu, crtlamd, - & crtlame, c0l + & crtlamd, + & crtlame ! ! real(kind=kind_phys) detad real(kind=kind_phys) adw, aup, aafac, d0, @@ -139,7 +139,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & & dh, dhh, dp, & dq, dqsdp, dqsdt, dt, & dt2, dtmax, dtmin, - & dxcrtas, dxcrtuf, dxcrtuf_hwrf, + & dxcrtas, dxcrtuf, & dv1h, dv2h, dv3h, & dv1q, dv2q, dv3q, & dz, dz1, e1, edtmax, @@ -204,7 +204,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & ! parameter(c0s=.002,c1=.002,d0=.01) ! parameter(d0=.01) parameter(d0=.001) -!mz parameter(c0l=c0s*asolfac) +! parameter(c0l=c0s*asolfac) ! ! asolfac: aerosol-aware parameter based on Lim (2011) ! asolfac= cx / c0s(=.002) @@ -223,7 +223,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & ! parameter(cinacrmx=-120.,cinacrmn=-120.) parameter(cinacrmx=-120.,cinacrmn=-80.) parameter(bet1=1.875,cd1=.506,f1=2.0,gam1=.5) - parameter(betaw=.03,dxcrtas=8.e3,dxcrtuf=15.e3,dxcrtuf_hwrf=25.e3) + parameter(betaw=.03,dxcrtas=8.e3,dxcrtuf=15.e3) ! ! local variables and arrays real(kind=kind_phys) pfld(im,km), to(im,km), qo(im,km), @@ -234,7 +234,6 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & ! for updraft velocity calculation real(kind=kind_phys) wu2(im,km), buo(im,km), drag(im,km) real(kind=kind_phys) wc(im), scaldfunc(im), sigmagfm(im) - real(kind=kind_phys) sigmuout(im) ! c cloud water ! real(kind=kind_phys) tvo(im,km) @@ -370,10 +369,10 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & vshear(i) = 0. gdx(i) = sqrt(garea(i)) - !mz*HWRF SAS + !HWRF SAS scaldfunc(i)=-1.0 sigmagfm(i)=-1.0 - sigmuout(i)=-1.0 +! sigmuout(i)=-1.0 enddo endif ! @@ -449,7 +448,6 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & edtmaxs = .3 if (hwrf_samfdeep) then aafac = .1 - crtlamu = 1.0e-4 cxlamu = 1.0e-3 else aafac = .05 @@ -840,7 +838,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & c assume that updraft entrainment rate above cloud base is c same as that at cloud base c -!> - Calculate the entrainment rate according to Han and Pan (2011) \cite han_and_pan_2011 , equation 8, after Bechtold et al. (2008) \cite bechtold_et_al_2008, equation 2 given by: +!> - In HWRF samfdeep, calculate the entrainment rate according to Han and Pan (2011) \cite han_and_pan_2011 , equation 8, after Bechtold et al. (2008) \cite bechtold_et_al_2008, equation 2 given by: !! \f[ !! \epsilon = \epsilon_0F_0 + d_1\left(1-RH\right)F_1 !! \f] @@ -866,17 +864,23 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & !! (The updraft detrainment rate is set constant and equal to the entrainment rate at cloud base.) !! !> - The updraft detrainment rate is vertically constant and proportional to clamt - do k = 1, km1 + if (hwrf_samfdeep) then + do k = 1, km1 do i=1,im if(cnvflg(i) .and. k < kmax(i)) then - if (hwrf_samfdeep) then xlamud(i,k) = xlamx(i) - else + endif + enddo + enddo + else + do k = 1, km1 + do i=1,im + if(cnvflg(i) .and. k < kmax(i)) then xlamud(i,k) = 0.001 * clamt(i) - endif endif enddo - enddo + enddo + endif c c entrainment functions decreasing with height (fent), c mimicking a cloud ensemble @@ -2503,6 +2507,18 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & !! !> - For scale-aware parameterization, the updraft fraction (sigmagfm) is first computed as a function of the lateral entrainment rate at cloud base (see Han et al.'s (2017) \cite han_et_al_2017 equation 4 and 5), following the study by Grell and Freitas (2014) \cite grell_and_freitas_2014. + if(hwrf_samfdeep) then + do i = 1, im + if(cnvflg(i)) then + tem = min(max(xlamx(i), 7.e-5), 3.e-4) + tem = 0.2 / tem + tem1 = 3.14 * tem * tem + sigmagfm(i) = tem1 / garea(i) + sigmagfm(i) = max(sigmagfm(i), 0.001) + sigmagfm(i) = min(sigmagfm(i), 0.999) + endif + enddo + else do i = 1, im if(cnvflg(i)) then tem = min(max(xlamue(i,kbcon(i)), 7.e-5), 3.e-4) @@ -2513,24 +2529,10 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & sigmagfm(i) = min(sigmagfm(i), 0.999) endif enddo + endif ! !> - Then, calculate the reduction factor (scaldfunc) of the vertical convective eddy transport of mass flux as a function of updraft fraction from the studies by Arakawa and Wu (2013) \cite arakawa_and_wu_2013 (also see Han et al.'s (2017) \cite han_et_al_2017 equation 1 and 2). The final cloud base mass flux with scale-aware parameterization is obtained from the mass flux when sigmagfm << 1, multiplied by the reduction factor (Han et al.'s (2017) \cite han_et_al_2017 equation 2). - if(hwrf_samfdeep) then - do i = 1, im - if(cnvflg(i)) then - if (gdx(i) < dxcrtuf) then - scaldfunc(i) = (1.-sigmagfm(i)) * (1.-sigmagfm(i)) - scaldfunc(i) = max(min(scaldfunc(i), 1.0), 0.) - sigmuout(i)=sigmagfm(i) - else - scaldfunc(i) = 1.0 - endif - xmb(i) = xmb(i) * scaldfunc(i) - xmb(i) = min(xmb(i),xmbmax(i)) - endif - enddo - else do i = 1, im if(cnvflg(i)) then if (gdx(i) < dxcrtuf) then @@ -2543,7 +2545,6 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & xmb(i) = min(xmb(i),xmbmax(i)) endif enddo - endif if (.not.hwrf_samfdeep) then !> - If stochastic physics using cellular automata is .true. then perturb the mass-flux here: From 322f5b17c13015b23e075463459b9077eb8943e3 Mon Sep 17 00:00:00 2001 From: Philip Pegion Date: Fri, 20 Mar 2020 14:03:42 +0000 Subject: [PATCH 07/14] make rain/snow tendency consistent with accumulated rain/snow --- physics/GFS_MP_generic.F90 | 15 ++++----------- 1 file changed, 4 insertions(+), 11 deletions(-) diff --git a/physics/GFS_MP_generic.F90 b/physics/GFS_MP_generic.F90 index f72f9405a..bbf88e24b 100644 --- a/physics/GFS_MP_generic.F90 +++ b/physics/GFS_MP_generic.F90 @@ -341,8 +341,10 @@ subroutine GFS_MP_generic_post_run(im, ix, levs, kdt, nrcm, ncld, nncl, ntcw, nt if (cplflx .or. cplchm) then do i = 1, im - rain_cpl(i) = rain_cpl(i) + rain(i) * (one-srflag(i)) - snow_cpl(i) = snow_cpl(i) + rain(i) * srflag(i) + drain_cpl(i) = rain(i) * (one-srflag(i)) + dsnow_cpl(i) = rain(i) * srflag(i) + rain_cpl(i) = rain_cpl(i) + drain_cpl(i) + snow_cpl(i) = snow_cpl(i) + dsnow_cpl(i) enddo endif @@ -376,15 +378,6 @@ subroutine GFS_MP_generic_post_run(im, ix, levs, kdt, nrcm, ncld, nncl, ntcw, nt if (do_sppt) then !--- radiation heating rate dtdtr(1:im,:) = dtdtr(1:im,:) + dtdtc(1:im,:)*dtf - do i = 1, im - if (t850(i) > 273.16) then -!--- change in change in rain precip - drain_cpl(i) = rain(i) - drain_cpl(i) - else -!--- change in change in snow precip - dsnow_cpl(i) = rain(i) - dsnow_cpl(i) - endif - enddo endif end subroutine GFS_MP_generic_post_run From f143b81eefed08706757a11b1e11cd085ab8aa75 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Tue, 24 Mar 2020 09:26:25 -0600 Subject: [PATCH 08/14] Updates of CCPP code to regain bit-for-bit identical results for coupled model runs --- physics/GFS_MP_generic.F90 | 10 +++--- physics/GFS_PBL_generic.F90 | 21 ++++++----- physics/GFS_debug.F90 | 24 +++++++++++-- physics/GFS_suite_interstitial.F90 | 14 ++++---- physics/GFS_surface_composites.F90 | 56 ++++++++++++++++++------------ physics/GFS_surface_generic.F90 | 13 ++++--- physics/GFS_surface_generic.meta | 18 ++++++++++ physics/gcycle.F90 | 2 +- physics/sfc_nst.f | 12 +++---- physics/sfc_nst.meta | 18 ++++++++++ 10 files changed, 131 insertions(+), 57 deletions(-) diff --git a/physics/GFS_MP_generic.F90 b/physics/GFS_MP_generic.F90 index f72f9405a..e28f535de 100644 --- a/physics/GFS_MP_generic.F90 +++ b/physics/GFS_MP_generic.F90 @@ -191,11 +191,11 @@ subroutine GFS_MP_generic_post_run(im, ix, levs, kdt, nrcm, ncld, nncl, ntcw, nt end if if (lsm==lsm_ruc .or. lsm==lsm_noahmp) then - raincprv(:) = rainc(:) - rainncprv(:) = frain * rain1(:) - iceprv(:) = ice(:) - snowprv(:) = snow(:) - graupelprv(:) = graupel(:) + raincprv(:) = rainc(:) + rainncprv(:) = frain * rain1(:) + iceprv(:) = ice(:) + snowprv(:) = snow(:) + graupelprv(:) = graupel(:) !for NoahMP, calculate precipitation rates from liquid water equivalent thickness for use in next time step !Note (GJF): Precipitation LWE thicknesses are multiplied by the frain factor, and are thus on the dynamics time step, but the conversion as written ! (with dtp in the denominator) assumes the rate is calculated on the physics time step. This only works as expected when dtf=dtp (i.e. when frain=1). diff --git a/physics/GFS_PBL_generic.F90 b/physics/GFS_PBL_generic.F90 index a440836e1..ff59aa465 100644 --- a/physics/GFS_PBL_generic.F90 +++ b/physics/GFS_PBL_generic.F90 @@ -331,7 +331,10 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac, character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg - real(kind=kind_phys), parameter :: huge=1.0d30, epsln = 1.0d-10 + real(kind=kind_phys), parameter :: zero = 0.0d0 + real(kind=kind_phys), parameter :: one = 1.0d0 + real(kind=kind_phys), parameter :: huge = 9.9692099683868690E36 ! NetCDF float FillValue, same as in GFS_typedefs.F90 + real(kind=kind_phys), parameter :: epsln = 1.0d-10 ! same as in GFS_physics_driver.F90 integer :: i, k, kk, k1, n real(kind=kind_phys) :: tem, tem1, rho @@ -486,7 +489,7 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac, if (cplchm) then do i = 1, im tem1 = max(q1(i), 1.e-8) - tem = prsl(i,1) / (rd*t1(i)*(1.0+fvirt*tem1)) + tem = prsl(i,1) / (rd*t1(i)*(one+fvirt*tem1)) ushfsfci(i) = -cp * tem * hflx(i) ! upward sensible heat flux enddo ! dkt_cpl has dimensions (1:im,1:levs), but dkt has (1:im,1:levs-1) @@ -498,22 +501,22 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac, if (cplflx) then do i=1,im - if (oceanfrac(i) > 0.0) then ! Ocean only, NO LAKES - if (fice(i) > 1.-epsln) then ! no open water, use results from CICE + if (oceanfrac(i) > zero) then ! Ocean only, NO LAKES + if (fice(i) > one - epsln) then ! no open water, use results from CICE dusfci_cpl(i) = dusfc_cice(i) dvsfci_cpl(i) = dvsfc_cice(i) dtsfci_cpl(i) = dtsfc_cice(i) dqsfci_cpl(i) = dqsfc_cice(i) - elseif (dry(i) .or. icy(i)) then ! use stress_ocean from sfc_diff for opw component at mixed point + elseif (icy(i) .or. dry(i)) then ! use stress_ocean from sfc_diff for opw component at mixed point tem1 = max(q1(i), 1.e-8) - rho = prsl(i,1) / (rd*t1(i)*(1.0+fvirt*tem1)) - if (wind(i) > 0.0) then + rho = prsl(i,1) / (rd*t1(i)*(one+fvirt*tem1)) + if (wind(i) > zero) then tem = - rho * stress_ocn(i) / wind(i) dusfci_cpl(i) = tem * ugrs1(i) ! U-momentum flux dvsfci_cpl(i) = tem * vgrs1(i) ! V-momentum flux else - dusfci_cpl(i) = 0.0 - dvsfci_cpl(i) = 0.0 + dusfci_cpl(i) = zero + dvsfci_cpl(i) = zero endif dtsfci_cpl(i) = cp * rho * hflx_ocn(i) ! sensible heat flux over open ocean dqsfci_cpl(i) = hvap * rho * evap_ocn(i) ! latent heat flux over open ocean diff --git a/physics/GFS_debug.F90 b/physics/GFS_debug.F90 index df56cc069..6bf39d491 100644 --- a/physics/GFS_debug.F90 +++ b/physics/GFS_debug.F90 @@ -402,7 +402,12 @@ subroutine GFS_diagtoscreen_run (Model, Statein, Stateout, Sfcprop, Coupling, call print_var(mpirank,omprank, blkno, 'Coupling%rain_cpl', Coupling%rain_cpl) call print_var(mpirank,omprank, blkno, 'Coupling%snow_cpl', Coupling%snow_cpl) end if + if (Model%cplwav2atm) then + call print_var(mpirank,omprank, blkno, 'Coupling%zorlwav_cpl' , Coupling%zorlwav_cpl ) + end if if (Model%cplflx) then + call print_var(mpirank,omprank, blkno, 'Coupling%oro_cpl' , Coupling%oro_cpl ) + call print_var(mpirank,omprank, blkno, 'Coupling%slmsk_cpl' , Coupling%slmsk_cpl ) call print_var(mpirank,omprank, blkno, 'Coupling%slimskin_cpl', Coupling%slimskin_cpl ) call print_var(mpirank,omprank, blkno, 'Coupling%dusfcin_cpl ', Coupling%dusfcin_cpl ) call print_var(mpirank,omprank, blkno, 'Coupling%dvsfcin_cpl ', Coupling%dvsfcin_cpl ) @@ -466,11 +471,24 @@ subroutine GFS_diagtoscreen_run (Model, Statein, Stateout, Sfcprop, Coupling, call print_var(mpirank,omprank, blkno, 'Coupling%shum_wts', Coupling%shum_wts) end if if (Model%do_skeb) then - call print_var(mpirank,omprank, blkno, 'Coupling%skebu_wts', Coupling%skebu_wts) - call print_var(mpirank,omprank, blkno, 'Coupling%skebv_wts', Coupling%skebv_wts) + call print_var(mpirank,omprank, blkno, 'Coupling%skebu_wts', Coupling%skebu_wts ) + call print_var(mpirank,omprank, blkno, 'Coupling%skebv_wts', Coupling%skebv_wts ) end if if (Model%do_sfcperts) then - call print_var(mpirank,omprank, blkno, 'Coupling%sfc_wts', Coupling%sfc_wts) + call print_var(mpirank,omprank, blkno, 'Coupling%sfc_wts' , Coupling%sfc_wts ) + end if + if (Model%do_ca) then + call print_var(mpirank,omprank, blkno, 'Coupling%tconvtend', Coupling%tconvtend ) + call print_var(mpirank,omprank, blkno, 'Coupling%qconvtend', Coupling%qconvtend ) + call print_var(mpirank,omprank, blkno, 'Coupling%uconvtend', Coupling%uconvtend ) + call print_var(mpirank,omprank, blkno, 'Coupling%vconvtend', Coupling%vconvtend ) + call print_var(mpirank,omprank, blkno, 'Coupling%ca_out ', Coupling%ca_out ) + call print_var(mpirank,omprank, blkno, 'Coupling%ca_deep ', Coupling%ca_deep ) + call print_var(mpirank,omprank, blkno, 'Coupling%ca_turb ', Coupling%ca_turb ) + call print_var(mpirank,omprank, blkno, 'Coupling%ca_shal ', Coupling%ca_shal ) + call print_var(mpirank,omprank, blkno, 'Coupling%ca_rad ', Coupling%ca_rad ) + call print_var(mpirank,omprank, blkno, 'Coupling%ca_micro ', Coupling%ca_micro ) + call print_var(mpirank,omprank, blkno, 'Coupling%cape ', Coupling%cape ) end if if(Model%imp_physics == Model%imp_physics_thompson .and. Model%ltaerosol) then call print_var(mpirank,omprank, blkno, 'Coupling%nwfa2d', Coupling%nwfa2d) diff --git a/physics/GFS_suite_interstitial.F90 b/physics/GFS_suite_interstitial.F90 index 8abaf24b7..935dd9430 100644 --- a/physics/GFS_suite_interstitial.F90 +++ b/physics/GFS_suite_interstitial.F90 @@ -228,15 +228,15 @@ subroutine GFS_suite_interstitial_2_run (im, levs, lssav, ldiag3d, lsidea, cplfl if (frac_grid) then do i=1,im - tem = one - cice(i) - frland(i) + tem = (one - frland(i)) * cice(i) ! tem = ice fraction wrt whole cell if (flag_cice(i)) then - adjsfculw(i) = adjsfculw_lnd(i) * frland(i) & - + ulwsfc_cice(i) * cice(i) & - + adjsfculw_ocn(i) * tem + adjsfculw(i) = adjsfculw_lnd(i) * frland(i) & + + ulwsfc_cice(i) * tem & + + adjsfculw_ocn(i) * (one - frland(i) - tem) else - adjsfculw(i) = adjsfculw_lnd(i) * frland(i) & - + adjsfculw_ice(i) * cice(i) & - + adjsfculw_ocn(i) * tem + adjsfculw(i) = adjsfculw_lnd(i) * frland(i) & + + adjsfculw_ice(i) * tem & + + adjsfculw_ocn(i) * (one - frland(i) - tem) endif enddo else diff --git a/physics/GFS_surface_composites.F90 b/physics/GFS_surface_composites.F90 index 6cca60ccf..b6d833796 100644 --- a/physics/GFS_surface_composites.F90 +++ b/physics/GFS_surface_composites.F90 @@ -89,7 +89,7 @@ subroutine GFS_surface_composites_pre_run (im, frac_grid, flag_cice, cplflx, cpl endif endif if (cice(i) < one ) then - wet(i)=.true. !there is some open ocean/lake water! + wet(i)=.true. ! some open ocean/lake water exists if (.not. cplflx) tsfco(i) = max(tsfco(i), tisfc(i), tgice) end if else @@ -414,7 +414,7 @@ subroutine GFS_surface_composites_post_run ( fm10(i) = fm10_lnd(i) fh2(i) = fh2_lnd(i) !tsurf(i) = tsurf_lnd(i) - tsfcl(i) = tsfc_lnd(i) + tsfcl(i) = tsfc_lnd(i) ! over land cmm(i) = cmm_lnd(i) chh(i) = chh_lnd(i) gflx(i) = gflx_lnd(i) @@ -426,9 +426,9 @@ subroutine GFS_surface_composites_post_run ( hflx(i) = hflx_lnd(i) qss(i) = qss_lnd(i) tsfc(i) = tsfc_lnd(i) - hice(i) = zero - cice(i) = zero - tisfc(i) = tsfc(i) + !hice(i) = zero + !cice(i) = zero + !tisfc(i) = tsfc(i) elseif (islmsk(i) == 0) then zorl(i) = zorl_ocn(i) cd(i) = cd_ocn(i) @@ -441,7 +441,7 @@ subroutine GFS_surface_composites_post_run ( fm10(i) = fm10_ocn(i) fh2(i) = fh2_ocn(i) !tsurf(i) = tsurf_ocn(i) - tsfco(i) = tsfc_ocn(i) + tsfco(i) = tsfc_ocn(i) ! over lake (and ocean when uncoupled) cmm(i) = cmm_ocn(i) chh(i) = chh_ocn(i) gflx(i) = gflx_ocn(i) @@ -453,10 +453,10 @@ subroutine GFS_surface_composites_post_run ( hflx(i) = hflx_ocn(i) qss(i) = qss_ocn(i) tsfc(i) = tsfc_ocn(i) - hice(i) = zero - cice(i) = zero - tisfc(i) = tsfc(i) - else + !hice(i) = zero + !cice(i) = zero + !tisfc(i) = tsfc(i) + else ! islmsk(i) == 2 zorl(i) = zorl_ice(i) cd(i) = cd_ice(i) cdq(i) = cdq_ice(i) @@ -468,30 +468,42 @@ subroutine GFS_surface_composites_post_run ( fm10(i) = fm10_ice(i) fh2(i) = fh2_ice(i) !tsurf(i) = tsurf_ice(i) + if (.not. flag_cice(i)) then + tisfc(i) = tice(i) ! over lake ice (and sea ice when uncoupled) + endif cmm(i) = cmm_ice(i) chh(i) = chh_ice(i) gflx(i) = gflx_ice(i) ep1d(i) = ep1d_ice(i) weasd(i) = weasd_ice(i) snowd(i) = snowd_ice(i) + !tprcp(i) = cice(i)*tprcp_ice(i) + (one-cice(i))*tprcp_ocn(i) qss(i) = qss_ice(i) - if (flag_cice(i)) then ! this was already done for lake ice in sfc_sice - txi = cice(i) - txo = one - txi - evap(i) = txi * evap_ice(i) + txo * evap_ocn(i) - hflx(i) = txi * hflx_ice(i) + txo * hflx_ocn(i) - tsfc(i) = txi * tsfc_ice(i) + txo * tsfc_ocn(i) - else - evap(i) = evap_ice(i) - hflx(i) = hflx_ice(i) - tsfc(i) = tsfc_ice(i) - tisfc(i) = tice(i) - endif + evap(i) = evap_ice(i) + hflx(i) = hflx_ice(i) + qss(i) = qss_ice(i) + tsfc(i) = tsfc_ice(i) endif zorll(i) = zorl_lnd(i) zorlo(i) = zorl_ocn(i) + if (flag_cice(i) .and. wet(i)) then ! this was already done for lake ice in sfc_sice + txi = cice(i) + txo = one - txi + evap(i) = txi * evap_ice(i) + txo * evap_ocn(i) + hflx(i) = txi * hflx_ice(i) + txo * hflx_ocn(i) + tsfc(i) = txi * tsfc_ice(i) + txo * tsfc_ocn(i) + else + if (islmsk(i) == 2) then + tisfc(i) = tice(i) + else ! over open ocean or land (no ice fraction) + hice(i) = zero + cice(i) = zero + tisfc(i) = tsfc(i) + endif + endif + enddo endif ! if (frac_grid) diff --git a/physics/GFS_surface_generic.F90 b/physics/GFS_surface_generic.F90 index 108d3bee7..98653a052 100644 --- a/physics/GFS_surface_generic.F90 +++ b/physics/GFS_surface_generic.F90 @@ -33,7 +33,7 @@ subroutine GFS_surface_generic_pre_run (im, levs, vfrac, islmsk, isot, ivegsrc, cplflx, flag_cice, islmsk_cice,slimskin_cpl, dusfcin_cpl, dvsfcin_cpl, & dtsfcin_cpl, dqsfcin_cpl, ulwsfcin_cpl, ulwsfc_cice, dusfc_cice, dvsfc_cice, & dtsfc_cice, dqsfc_cice, tisfc, tsfco, fice, hice, & - wind, u1, v1, cnvwind, errmsg, errflg) + wind, u1, v1, cnvwind, smcwlt2, smcref2, errmsg, errflg) use surface_perturbation, only: cdfnor @@ -87,6 +87,8 @@ subroutine GFS_surface_generic_pre_run (im, levs, vfrac, islmsk, isot, ivegsrc, real(kind=kind_phys), dimension(im), intent(in ) :: u1, v1 ! surface wind enhancement due to convection real(kind=kind_phys), dimension(im), intent(inout ) :: cnvwind + ! + real(kind=kind_phys), dimension(im), intent(out) :: smcwlt2, smcref2 ! CCPP error handling character(len=*), intent(out) :: errmsg @@ -173,7 +175,10 @@ subroutine GFS_surface_generic_pre_run (im, levs, vfrac, islmsk, isot, ivegsrc, work3(i) = prsik_1(i) / prslk_1(i) !tsurf(i) = tsfc(i) - zlvl(i) = phil(i,1) * onebg + zlvl(i) = phil(i,1) * onebg + smcwlt2(i) = zero + smcref2(i) = zero + wind(i) = max(sqrt(u1(i)*u1(i) + v1(i)*v1(i)) & + max(zero, min(cnvwind(i), 30.0)), one) !wind(i) = max(sqrt(Statein%ugrs(i,1)*Statein%ugrs(i,1) + & @@ -303,8 +308,8 @@ subroutine GFS_surface_generic_post_run (im, cplflx, cplwav, lssav, icy, wet, dt nlwsfc_cpl (i) = nlwsfc_cpl(i) + nlwsfci_cpl(i)*dtf t2mi_cpl (i) = t2m(i) q2mi_cpl (i) = q2m(i) -! tsfci_cpl (i) = tsfc(i) - tsfci_cpl (i) = tsfc_ocn(i) + tsfci_cpl (i) = tsfc(i) +! tsfci_cpl (i) = tsfc_ocn(i) psurfi_cpl (i) = pgr(i) enddo diff --git a/physics/GFS_surface_generic.meta b/physics/GFS_surface_generic.meta index 6bd18a3b8..d82928d9c 100644 --- a/physics/GFS_surface_generic.meta +++ b/physics/GFS_surface_generic.meta @@ -545,6 +545,24 @@ kind = kind_phys intent = inout optional = F +[smcwlt2] + standard_name = volume_fraction_of_condensed_water_in_soil_at_wilting_point + long_name = wilting point (volumetric) + units = frac + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = out + optional = F +[smcref2] + standard_name = threshold_volume_fraction_of_condensed_water_in_soil + long_name = soil moisture threshold (volumetric) + units = frac + dimensions = (horizontal_dimension) + 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/gcycle.F90 b/physics/gcycle.F90 index bb1730fc2..8c5dd041a 100644 --- a/physics/gcycle.F90 +++ b/physics/gcycle.F90 @@ -200,7 +200,7 @@ SUBROUTINE GCYCLE (nblks, Model, Grid, Sfcprop, Cldprop) Sfcprop(nb)%tref(ix) = TSFFCS (len) ! if ( Model%nstf_name(2) == 0 ) then ! dt_warm = (Sfcprop(nb)%xt(ix) + Sfcprop(nb)%xt(ix) ) & -! / Sfcprop(nb)%xz(ix) +! / Sfcprop(nb)%xz(ix) ! Sfcprop(nb)%tsfco(ix) = Sfcprop(nb)%tref(ix) & ! + dt_warm - Sfcprop(nb)%dt_cool(ix) ! endif diff --git a/physics/sfc_nst.f b/physics/sfc_nst.f index ed6387afb..3d0507ad9 100644 --- a/physics/sfc_nst.f +++ b/physics/sfc_nst.f @@ -676,7 +676,7 @@ end subroutine sfc_nst_pre_finalize !! @{ subroutine sfc_nst_pre_run & (im, wet, tsfc_ocn, tsurf_ocn, tseal, xt, xz, dt_cool, - & z_c, tref, cplflx, errmsg, errflg) + & z_c, tref, cplflx, oceanfrac, errmsg, errflg) use machine , only : kind_phys @@ -686,7 +686,7 @@ subroutine sfc_nst_pre_run integer, intent(in) :: im logical, dimension(im), intent(in) :: wet real (kind=kind_phys), dimension(im), intent(in) :: - & tsfc_ocn, xt, xz, dt_cool, z_c + & tsfc_ocn, xt, xz, dt_cool, z_c, oceanfrac logical, intent(in) :: cplflx ! --- input/outputs: @@ -724,7 +724,7 @@ subroutine sfc_nst_pre_run if (cplflx) then tem1 = half / omz1 do i=1,im - if (wet(i)) then + if (wet(i) .and. oceanfrac(i) > zero) then tem2 = one / xz(i) dt_warm = (xt(i)+xt(i)) * tem2 if ( xz(i) > omz1) then @@ -777,7 +777,7 @@ end subroutine sfc_nst_post_finalize ! \section NSST_detailed_post_algorithm Detailed Algorithm ! @{ subroutine sfc_nst_post_run & - & ( im, rlapse, wet, icy, oro, oro_uf, nstf_name1, & + & ( im, rlapse, tgice, wet, icy, oro, oro_uf, nstf_name1, & & nstf_name4, nstf_name5, xt, xz, dt_cool, z_c, tref, xlon, & & tsurf_ocn, tsfc_ocn, dtzm, errmsg, errflg & & ) @@ -790,7 +790,7 @@ subroutine sfc_nst_post_run & ! --- inputs: integer, intent(in) :: im logical, dimension(im), intent(in) :: wet, icy - real (kind=kind_phys), intent(in) :: rlapse + real (kind=kind_phys), intent(in) :: rlapse, tgice real (kind=kind_phys), dimension(im), intent(in) :: oro, oro_uf integer, intent(in) :: nstf_name1, nstf_name4, nstf_name5 real (kind=kind_phys), dimension(im), intent(in) :: xt, xz, & @@ -838,7 +838,7 @@ subroutine sfc_nst_post_run & ! if (wet(i) .and. .not.icy(i)) then ! if (wet(i) .and. (Model%frac_grid .or. .not. icy(i))) then if (wet(i)) then - tsfc_ocn(i) = max(271.2, tref(i) + dtzm(i)) + tsfc_ocn(i) = max(tgice, tref(i) + dtzm(i)) ! tsfc_ocn(i) = max(271.2, tref(i) + dtzm(i)) - & ! (oro(i)-oro_uf(i))*rlapse endif diff --git a/physics/sfc_nst.meta b/physics/sfc_nst.meta index d74f68c0e..ac75aa05d 100644 --- a/physics/sfc_nst.meta +++ b/physics/sfc_nst.meta @@ -759,6 +759,15 @@ type = logical intent = in optional = F +[oceanfrac] + standard_name = sea_area_fraction + long_name = fraction of horizontal grid area occupied by ocean + units = frac + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP @@ -808,6 +817,15 @@ kind = kind_phys intent = in optional = F +[tgice] + standard_name = freezing_point_temperature_of_seawater + long_name = freezing point temperature of seawater + units = K + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F [wet] standard_name = flag_nonzero_wet_surface_fraction long_name = flag indicating presence of some ocean or lake surface area fraction From 091d47524aedd7e62d5e26afc4da85e7cdd45a1c Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Tue, 24 Mar 2020 16:50:33 -0600 Subject: [PATCH 09/14] physics/GFS_stochastics.F90: update comment --- physics/GFS_stochastics.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/physics/GFS_stochastics.F90 b/physics/GFS_stochastics.F90 index 2a6552f18..99f84e3b1 100644 --- a/physics/GFS_stochastics.F90 +++ b/physics/GFS_stochastics.F90 @@ -79,10 +79,10 @@ subroutine GFS_stochastics_run (im, km, do_sppt, use_zmtnblck, do_shum, do_skeb, real(kind_phys), dimension(1:im), intent(inout) :: totprcpb real(kind_phys), dimension(1:im), intent(inout) :: cnvprcpb logical, intent(in) :: cplflx - ! rain_cpl, snow_cpl only allocated if cplflx == .true. or do_sppt == .true. + ! rain_cpl, snow_cpl only allocated if cplflx == .true. or cplchm == .true. real(kind_phys), dimension(:), intent(inout) :: rain_cpl real(kind_phys), dimension(:), intent(inout) :: snow_cpl - ! drain_cpl, dsnow_cpl only allocated if do_sppt == .true. + ! drain_cpl, dsnow_cpl only allocated if cplflx == .true. or cplchm == .true. real(kind_phys), dimension(:), intent(in) :: drain_cpl real(kind_phys), dimension(:), intent(in) :: dsnow_cpl ! tconvtend ... vconvtend only allocated if isppt_deep == .true. From bde224d6e0c3b092a91e3518ea3a1b238d22a4c7 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Wed, 25 Mar 2020 13:17:08 -0600 Subject: [PATCH 10/14] Remove interstitial variables for seaice coupling --- physics/GFS_PBL_generic.meta | 16 +++--- physics/GFS_surface_generic.F90 | 20 ++----- physics/GFS_surface_generic.meta | 90 -------------------------------- physics/sfc_cice.meta | 16 +++--- 4 files changed, 20 insertions(+), 122 deletions(-) diff --git a/physics/GFS_PBL_generic.meta b/physics/GFS_PBL_generic.meta index 51764e04d..2319f0044 100644 --- a/physics/GFS_PBL_generic.meta +++ b/physics/GFS_PBL_generic.meta @@ -1089,8 +1089,8 @@ intent = in optional = F [dusfc_cice] - standard_name = surface_x_momentum_flux_for_coupling_interstitial - long_name = sfc x momentum flux for coupling interstitial + standard_name = surface_x_momentum_flux_for_coupling + long_name = sfc x momentum flux for coupling units = Pa dimensions = (horizontal_dimension) type = real @@ -1098,8 +1098,8 @@ intent = in optional = F [dvsfc_cice] - standard_name = surface_y_momentum_flux_for_coupling_interstitial - long_name = sfc y momentum flux for coupling interstitial + standard_name = surface_y_momentum_flux_for_coupling + long_name = sfc y momentum flux for coupling units = Pa dimensions = (horizontal_dimension) type = real @@ -1107,8 +1107,8 @@ intent = in optional = F [dtsfc_cice] - standard_name = surface_upward_sensible_heat_flux_for_coupling_interstitial - long_name = sfc sensible heat flux for coupling interstitial + standard_name = surface_upward_sensible_heat_flux_for_coupling + long_name = sfc sensible heat flux for coupling units = W m-2 dimensions = (horizontal_dimension) type = real @@ -1116,8 +1116,8 @@ intent = in optional = F [dqsfc_cice] - standard_name = surface_upward_latent_heat_flux_for_coupling_interstitial - long_name = sfc latent heat flux for coupling interstitial + standard_name = surface_upward_latent_heat_flux_for_coupling + long_name = sfc latent heat flux for coupling units = W m-2 dimensions = (horizontal_dimension) type = real diff --git a/physics/GFS_surface_generic.F90 b/physics/GFS_surface_generic.F90 index 98653a052..3b52677a9 100644 --- a/physics/GFS_surface_generic.F90 +++ b/physics/GFS_surface_generic.F90 @@ -30,9 +30,7 @@ subroutine GFS_surface_generic_pre_run (im, levs, vfrac, islmsk, isot, ivegsrc, sigmaf, soiltyp, vegtype, slopetyp, work3, tsurf, zlvl, do_sppt, dtdtr, & drain_cpl, dsnow_cpl, rain_cpl, snow_cpl, do_sfcperts, nsfcpert, sfc_wts, & pertz0, pertzt, pertshc, pertlai, pertvegf, z01d, zt1d, bexp1d, xlai1d, vegf1d, & - cplflx, flag_cice, islmsk_cice,slimskin_cpl, dusfcin_cpl, dvsfcin_cpl, & - dtsfcin_cpl, dqsfcin_cpl, ulwsfcin_cpl, ulwsfc_cice, dusfc_cice, dvsfc_cice, & - dtsfc_cice, dqsfc_cice, tisfc, tsfco, fice, hice, & + cplflx, flag_cice, islmsk_cice, slimskin_cpl, tisfc, tsfco, fice, hice, & wind, u1, v1, cnvwind, smcwlt2, smcref2, errmsg, errflg) use surface_perturbation, only: cdfnor @@ -76,12 +74,9 @@ subroutine GFS_surface_generic_pre_run (im, levs, vfrac, islmsk, isot, ivegsrc, logical, intent(in) :: cplflx real(kind=kind_phys), dimension(im), intent(in) :: slimskin_cpl logical, dimension(im), intent(inout) :: flag_cice - integer, dimension(im), intent(out) :: islmsk_cice - real(kind=kind_phys), dimension(im), intent(in) ::ulwsfcin_cpl, & - dusfcin_cpl, dvsfcin_cpl, dtsfcin_cpl, dqsfcin_cpl, & + integer, dimension(im), intent(out) :: islmsk_cice + real(kind=kind_phys), dimension(im), intent(in) :: & tisfc, tsfco, fice, hice - real(kind=kind_phys), dimension(im), intent(out) ::ulwsfc_cice, & - dusfc_cice, dvsfc_cice, dtsfc_cice, dqsfc_cice real(kind=kind_phys), dimension(im), intent(out) :: wind real(kind=kind_phys), dimension(im), intent(in ) :: u1, v1 @@ -191,14 +186,7 @@ subroutine GFS_surface_generic_pre_run (im, levs, vfrac, islmsk, isot, ivegsrc, if (cplflx) then do i=1,im islmsk_cice(i) = nint(slimskin_cpl(i)) - if(islmsk_cice(i) == 4)then - flag_cice(i) = .true. - ulwsfc_cice(i) = ulwsfcin_cpl(i) - dusfc_cice(i) = dusfcin_cpl(i) - dvsfc_cice(i) = dvsfcin_cpl(i) - dtsfc_cice(i) = dtsfcin_cpl(i) - dqsfc_cice(i) = dqsfcin_cpl(i) - endif + flag_cice(i) = (islmsk_cice(i) == 4) enddo endif diff --git a/physics/GFS_surface_generic.meta b/physics/GFS_surface_generic.meta index d82928d9c..250f7a2bd 100644 --- a/physics/GFS_surface_generic.meta +++ b/physics/GFS_surface_generic.meta @@ -383,96 +383,6 @@ kind = kind_phys intent = in optional = F -[dusfcin_cpl] - standard_name = surface_x_momentum_flux_for_coupling - long_name = sfc x momentum flux for coupling - units = Pa - dimensions = (horizontal_dimension) - type = real - kind = kind_phys - intent = in - optional = F -[dvsfcin_cpl] - standard_name = surface_y_momentum_flux_for_coupling - long_name = sfc y momentum flux for coupling - units = Pa - dimensions = (horizontal_dimension) - type = real - kind = kind_phys - intent = in - optional = F -[dtsfcin_cpl] - standard_name = surface_upward_sensible_heat_flux_for_coupling - long_name = sfc sensible heat flux input - units = W m-2 - dimensions = (horizontal_dimension) - type = real - kind = kind_phys - intent = in - optional = F -[dqsfcin_cpl] - standard_name = surface_upward_latent_heat_flux_for_coupling - long_name = sfc latent heat flux input for coupling - units = W m-2 - dimensions = (horizontal_dimension) - type = real - kind = kind_phys - intent = in - optional = F -[ulwsfcin_cpl] - standard_name = surface_upwelling_longwave_flux_for_coupling - long_name = surface upwelling LW flux for coupling - units = W m-2 - dimensions = (horizontal_dimension) - type = real - kind = kind_phys - intent = in - optional = F -[ulwsfc_cice] - standard_name = surface_upwelling_longwave_flux_for_coupling_interstitial - long_name = surface upwelling longwave flux for coupling interstitial - units = W m-2 - dimensions = (horizontal_dimension) - type = real - kind = kind_phys - intent = out - optional = F -[dusfc_cice] - standard_name = surface_x_momentum_flux_for_coupling_interstitial - long_name = sfc x momentum flux for coupling interstitial - units = Pa - dimensions = (horizontal_dimension) - type = real - kind = kind_phys - intent = out - optional = F -[dvsfc_cice] - standard_name = surface_y_momentum_flux_for_coupling_interstitial - long_name = sfc y momentum flux for coupling interstitial - units = Pa - dimensions = (horizontal_dimension) - type = real - kind = kind_phys - intent = out - optional = F -[dtsfc_cice] - standard_name = surface_upward_sensible_heat_flux_for_coupling_interstitial - long_name = sfc sensible heat flux for coupling interstitial - units = W m-2 - dimensions = (horizontal_dimension) - type = real - kind = kind_phys - intent = out - optional = F -[dqsfc_cice] - standard_name = surface_upward_latent_heat_flux_for_coupling_interstitial - long_name = sfc latent heat flux for coupling interstitial - units = W m-2 - dimensions = (horizontal_dimension) - type = real - kind = kind_phys - intent = out - optional = F [tisfc] standard_name = sea_ice_temperature long_name = sea-ice surface temperature diff --git a/physics/sfc_cice.meta b/physics/sfc_cice.meta index 543e4d78b..a1c57d4d9 100644 --- a/physics/sfc_cice.meta +++ b/physics/sfc_cice.meta @@ -124,8 +124,8 @@ intent = in optional = F [dqsfc] - standard_name = surface_upward_latent_heat_flux_for_coupling_interstitial - long_name = sfc latent heat flux for coupling interstitial + standard_name = surface_upward_latent_heat_flux_for_coupling + long_name = sfc latent heat flux for coupling units = W m-2 dimensions = (horizontal_dimension) type = real @@ -133,8 +133,8 @@ intent = in optional = F [dtsfc] - standard_name = surface_upward_sensible_heat_flux_for_coupling_interstitial - long_name = sfc sensible heat flux for coupling interstitial + standard_name = surface_upward_sensible_heat_flux_for_coupling + long_name = sfc sensible heat flux for coupling units = W m-2 dimensions = (horizontal_dimension) type = real @@ -142,8 +142,8 @@ intent = in optional = F [dusfc] - standard_name = surface_x_momentum_flux_for_coupling_interstitial - long_name = sfc x momentum flux for coupling interstitial + standard_name = surface_x_momentum_flux_for_coupling + long_name = sfc x momentum flux for coupling units = Pa dimensions = (horizontal_dimension) type = real @@ -151,8 +151,8 @@ intent = in optional = F [dvsfc] - standard_name = surface_y_momentum_flux_for_coupling_interstitial - long_name = sfc y momentum flux for coupling interstitial + standard_name = surface_y_momentum_flux_for_coupling + long_name = sfc y momentum flux for coupling units = Pa dimensions = (horizontal_dimension) type = real From 5c134c17d88f9ed008e0e4c0bbab392b2c1f4d13 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Wed, 25 Mar 2020 14:10:44 -0600 Subject: [PATCH 11/14] physics/GFS_surface_generic.F90: remove old code that no longer exists in IPD --- physics/GFS_surface_generic.F90 | 4 ---- 1 file changed, 4 deletions(-) diff --git a/physics/GFS_surface_generic.F90 b/physics/GFS_surface_generic.F90 index 3b52677a9..ac366ae54 100644 --- a/physics/GFS_surface_generic.F90 +++ b/physics/GFS_surface_generic.F90 @@ -104,10 +104,6 @@ subroutine GFS_surface_generic_pre_run (im, levs, vfrac, islmsk, isot, ivegsrc, ! Set initial quantities for stochastic physics deltas if (do_sppt) then dtdtr = 0.0 - do i=1,im - drain_cpl(i) = rain_cpl (i) - dsnow_cpl(i) = snow_cpl (i) - enddo endif ! Scale random patterns for surface perturbations with perturbation size From bccf301e158c346b051b3d6a9bd392a71d07df25 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Wed, 1 Apr 2020 14:12:13 -0600 Subject: [PATCH 12/14] physics/samfdeepcnv.f: bugfix, ca_deep only allocated when do_ca is .true. --- physics/samfdeepcnv.f | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/physics/samfdeepcnv.f b/physics/samfdeepcnv.f index b0813d98b..f64a0b332 100644 --- a/physics/samfdeepcnv.f +++ b/physics/samfdeepcnv.f @@ -92,8 +92,9 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & real(kind=kind_phys), intent(in) :: psp(im), delp(ix,km), & & prslp(ix,km), garea(im), dot(ix,km), phil(ix,km) real(kind=kind_phys), dimension(:), intent(in) :: fscav - real(kind=kind_phys), intent(in) :: ca_deep(ix) logical, intent(in) :: do_ca, hwrf_samfdeep + ! ca_deep only allocatedd when do_ca is true + real(kind=kind_phys), intent(in) :: ca_deep(:) integer, intent(inout) :: kcnv(im) ! DH* TODO - check dimensions of qtr, ntr+2 correct? *DH @@ -1640,6 +1641,7 @@ subroutine samfdeepcnv_run (im,ix,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo + if (hwrf_samfdeep) then do i = 1, im beta = betas From 71eace19a5f42b60f64816575c425d34624c6018 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Fri, 3 Apr 2020 15:41:13 -0600 Subject: [PATCH 13/14] physics/samfshalcnv.f: bugfix, move assignment inside if block as in previous version --- physics/samfshalcnv.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/physics/samfshalcnv.f b/physics/samfshalcnv.f index e21110bd6..7a6db70f0 100644 --- a/physics/samfshalcnv.f +++ b/physics/samfshalcnv.f @@ -966,8 +966,8 @@ subroutine samfshalcnv_run(im,ix,km,itc,ntc,cliq,cp,cvap, & tem = 1. - tem tem1= .5*(cinacrmx-cinacrmn) cinacr = cinacrmx - tem * tem1 - endif if(cina(i) < cinacr) cnvflg(i) = .false. + endif enddo endif !! From b61ea19fb484acc6bba891b8a12ef430cba8cdb2 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Mon, 6 Apr 2020 17:06:47 -0600 Subject: [PATCH 14/14] physics/GFS_debug.F90: add capability to debug 1-d logical arrays --- physics/GFS_debug.F90 | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/physics/GFS_debug.F90 b/physics/GFS_debug.F90 index 6bf39d491..3670f4ddd 100644 --- a/physics/GFS_debug.F90 +++ b/physics/GFS_debug.F90 @@ -18,6 +18,7 @@ module GFS_diagtoscreen interface print_var module procedure print_logic_0d + module procedure print_logic_1d module procedure print_int_0d module procedure print_int_1d module procedure print_real_0d @@ -116,6 +117,7 @@ subroutine GFS_diagtoscreen_run (Model, Statein, Stateout, Sfcprop, Coupling, do impi=0,mpisize-1 do iomp=0,ompsize-1 if (mpirank==impi .and. omprank==iomp) then + call print_var(mpirank,omprank, blkno, 'Model%kdt' , Model%kdt) ! Sfcprop call print_var(mpirank,omprank, blkno, 'Sfcprop%slmsk' , Sfcprop%slmsk) call print_var(mpirank,omprank, blkno, 'Sfcprop%oceanfrac', Sfcprop%oceanfrac) @@ -557,6 +559,30 @@ subroutine print_int_0d(mpirank,omprank,blkno,name,var) end subroutine print_int_0d + subroutine print_logic_1d(mpirank,omprank,blkno,name,var) + + use machine, only: kind_phys + + implicit none + + integer, intent(in) :: mpirank, omprank, blkno + character(len=*), intent(in) :: name + logical, intent(in) :: var(:) + + integer :: i + +#ifdef PRINT_SUM + write(0,'(2a,3i6,2i8)') 'XXX: ', trim(name), mpirank, omprank, blkno, size(var), count(var) +#elif defined(PRINT_CHKSUM) + write(0,'(2a,3i6,2i8)') 'XXX: ', trim(name), mpirank, omprank, blkno, size(var), count(var) +#else + do i=ISTART,min(IEND,size(var(:))) + write(0,'(2a,3i6,i6,1x,l)') 'XXX: ', trim(name), mpirank, omprank, blkno, i, var(i) + end do +#endif + + end subroutine print_logic_1d + subroutine print_int_1d(mpirank,omprank,blkno,name,var) use machine, only: kind_phys