From 326a7b7ff1b7700971f2cb14911eeda25ca10777 Mon Sep 17 00:00:00 2001 From: "Man.Zhang" Date: Wed, 20 Jun 2018 08:41:59 -0600 Subject: [PATCH 1/6] merge gwdc_run with NEMSfv3gfs, and get B4B results. --- physics/gscond.f | 60 +++++++------- physics/gwdc.f | 142 +++++++++++++++++---------------- physics/gwdps.f | 199 ++++++++++++++++++++++++++++++++++++++++++----- 3 files changed, 287 insertions(+), 114 deletions(-) diff --git a/physics/gscond.f b/physics/gscond.f index 1556db615..fc2525b0f 100644 --- a/physics/gscond.f +++ b/physics/gscond.f @@ -121,7 +121,7 @@ subroutine zhaocarr_gscond_run (im,ix,km,dt,dtf,prsl,ps,q,clw1 & real (kind=kind_phys) qi(im), qint(im), ccrik, e0 &, cond, rdt, us, cclimit, climit &, tmt0, tmt15, qik, cwmik - &, ai, qw, u00ik, tik, pres, pp0, fi + &, ai, qw, u00ik, tik, pres, pp0, fi &, at, aq, ap, fiw, elv, qc, rqik &, rqikk, tx1, tx2, tx3, es, qs &, tsq, delq, condi, cone0, us00, ccrik1 @@ -173,7 +173,7 @@ subroutine zhaocarr_gscond_run (im,ix,km,dt,dtf,prsl,ps,q,clw1 & ! endif ! !************************************************************* -!> -# Begining of grid-scale condensation/evaporation loop (start of +!> -# Begining of grid-scale condensation/evaporation loop (start of !! k-loop, i-loop) !************************************************************* ! @@ -182,7 +182,7 @@ subroutine zhaocarr_gscond_run (im,ix,km,dt,dtf,prsl,ps,q,clw1 & ! vprs(:) = 0.001 * fpvs(t(:,k)) ! fpvs in pa !----------------------------------------------------------------------- !------------------qw, qi and qint-------------------------------------- - do i = 1, im + do i = 1, im tmt0 = t(i,k)-273.16 tmt15 = min(tmt0,cons_m15) qik = max(q(i,k),epsq) @@ -214,23 +214,23 @@ subroutine zhaocarr_gscond_run (im,ix,km,dt,dtf,prsl,ps,q,clw1 & !! cloud identification number IW, which is zero for cloud water and !! unity for cloud ice (Table 2 in !! \cite zhao_and_carr_1997): -!! - All clouds are defined to consist of liquid water below the -!! freezing level (\f$T\geq 0^oC\f$) and of ice particles above the +!! - All clouds are defined to consist of liquid water below the +!! freezing level (\f$T\geq 0^oC\f$) and of ice particles above the !! \f$T=-15^oC\f$ level. -!! - In the temperature region between \f$-15^oC\f$ and \f$0^oC\f$, -!! clouds may be composed of liquid water or ice. If there are cloud +!! - In the temperature region between \f$-15^oC\f$ and \f$0^oC\f$, +!! clouds may be composed of liquid water or ice. If there are cloud !! ice particles above this point at the previous or current time step, -!! or if the cloud at this point at the previous time step consists of -!! ice particles, then the cloud substance at this point is considered +!! or if the cloud at this point at the previous time step consists of +!! ice particles, then the cloud substance at this point is considered !! to be ice particles because of the cloud seeding effect and the -!! memory of its content. Otherwise, all clouds in this region are +!! memory of its content. Otherwise, all clouds in this region are !! considered to contain supercooled cloud water. !-------------------ice-water id number iw------------------------------ if(tmt0.lt.-15.0) then u00ik = u(i,k) - fi = qik - u00ik*qi(i) - if(fi > d00.or.cwmik > climit) then + fi = qik - u00ik*qi(i) + if(fi > d00.or.cwmik > climit) then iw(i,k) = 1 else iw(i,k) = 0 @@ -251,8 +251,8 @@ subroutine zhaocarr_gscond_run (im,ix,km,dt,dtf,prsl,ps,q,clw1 & !> -# Condensation and evaporation of cloud !--------------condensation and evaporation of cloud-------------------- do i = 1, im -!> - Compute the changes in t, q and p (\f$A_{t}\f$,\f$A_{q}\f$ and -!! \f$A_{p}\f$) caused by all the processes except grid-scale +!> - Compute the changes in t, q and p (\f$A_{t}\f$,\f$A_{q}\f$ and +!! \f$A_{p}\f$) caused by all the processes except grid-scale !! condensation and evaporation. !!\f[ !! A_{t}=(t-tp)/dt @@ -274,7 +274,7 @@ subroutine zhaocarr_gscond_run (im,ix,km,dt,dtf,prsl,ps,q,clw1 & at = (tik-tp(i,k)) * rdt aq = (qik-qp(i,k)) * rdt ap = (pres-pp0) * rdt -!> - Calculate the saturation specific humidity \f$q_{s}\f$ and the +!> - Calculate the saturation specific humidity \f$q_{s}\f$ and the !! relative humidity \f$f\f$ using IW. !----------------the satuation specific humidity------------------------ fiw = float(iwik) @@ -283,7 +283,7 @@ subroutine zhaocarr_gscond_run (im,ix,km,dt,dtf,prsl,ps,q,clw1 & ! if (lprnt) print *,' qc=',qc,' qint=',qint(i),' qi=',qi(i) !----------------the relative humidity---------------------------------- if(qc.le.1.0e-10) then - rqik=d00 + rqik=d00 else rqik = qik/qc endif @@ -295,30 +295,30 @@ subroutine zhaocarr_gscond_run (im,ix,km,dt,dtf,prsl,ps,q,clw1 & !! b=1-\left ( \frac{f_{s}-f}{f_{s}-u} \right )^{1/2} !!\f] !! for \f$f>u\f$; and \f$b=0\f$ for \f$f climit) then + if (ccrik <= cclimit.and. cwmik > climit) then ! ! first iteration - increment halved ! @@ -438,7 +438,7 @@ subroutine zhaocarr_gscond_run (im,ix,km,dt,dtf,prsl,ps,q,clw1 & cond = d00 ! if (ccrik .gt. 0.20 .and. qc .gt. epsq) then if (ccrik .gt. cclimit .and. qc .gt. epsq) then - us00 = us - u00ik + us00 = us - u00ik ccrik1 = 1.0 - ccrik aa = eps*elv*pres*qik ab = ccrik*ccrik1*qc*us00 diff --git a/physics/gwdc.f b/physics/gwdc.f index dc67987ad..33d28b9c4 100644 --- a/physics/gwdc.f +++ b/physics/gwdc.f @@ -196,7 +196,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & ! 2013 S. Moorthi - Updated and optimized code for T1534 GFS implementation ! ??? ?? 2015 J. Alpert - reducing the magnitude of tauctmax to fix blow up in L64 GFS ! S. Kar & M. Young -! aug 15 2016 - S. Moorthi - Fix for exessive dissipation which led to blow up in +! aug 15 2016 - S. Moorthi - Fix for exessive dissipation which led to blow up in ! 128 level runs with NEMS/GSM !*********************************************************************** @@ -216,7 +216,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & ! dpmid : midpoint delta p ( pi(k)-pi(k-1) ) ! lat : latitude index ! qmax : deep convective heating -! kcldtop : Vertical level index for cloud top ( mid level ) +! kcldtop : Vertical level index for cloud top ( mid level ) ! kcldbot : Vertical level index for cloud bottom ( mid level ) ! kcnv : (0,1) dependent on whether convection occur or not ! @@ -257,19 +257,19 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & ! parallel to the wind vector at the cloud top ( mid level ) ! tauctx : Wave stress at the cloud top projected in the east ! taucty : Wave stress at the cloud top projected in the north -! qmax : Maximum deep convective heating rate ( K s-1 ) in a +! qmax : Maximum deep convective heating rate ( K s-1 ) in a ! horizontal grid point calculated from cumulus para- ! meterization. ( mid level ) ! wtgwc : Wind tendency in direction to the wind vector at the cloud top level ! due to convectively generated gravity waves ( mid level ) -! utgwcl : Zonal wind tendency due to convectively generated +! utgwcl : Zonal wind tendency due to convectively generated ! gravity waves ( mid level ) ! vtgwcl : Meridional wind tendency due to convectively generated ! gravity waves ( mid level ) ! taugwci : Profile of wave stress calculated using basic-wind -! parallel to the wind vector at the cloud top +! parallel to the wind vector at the cloud top ! taugwcxi : Profile of zonal component of gravity wave stress -! taugwcyi : Profile of meridional component of gravity wave stress +! taugwcyi : Profile of meridional component of gravity wave stress ! ! taugwci, taugwcxi, and taugwcyi are defined at the interface level ! @@ -279,7 +279,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & ! rhom : Air density ( mid level ) ! ti : Temperature ( interface level ) ! basicum : Basic-wind profile. Basic-wind is parallel to the wind -! vector at the cloud top level. (mid level) +! vector at the cloud top level. (mid level) ! basicui : Basic-wind profile. Basic-wind is parallel to the wind ! vector at the cloud top level. ( interface level ) ! riloc : Local Richardson number ( interface level ) @@ -289,9 +289,9 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & ! break : Horizontal location where wave breaking is occurred. ! critic : Horizontal location where critical level filtering is ! occurred. -! dogwdc : Logical flag whether the GWDC parameterization is +! dogwdc : Logical flag whether the GWDC parameterization is ! calculated at a grid point or not. -! +! ! dogwdc is used in order to lessen CPU time for GWDC calculation. ! !----------------------------------------------------------------------- @@ -327,7 +327,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & & ti(:,:), riloc(:,:), & rimin(:,:), pint(:,:) ! real(kind=kind_phys), allocatable :: ugwdc(:,:), vgwdc(:,:), - real(kind=kind_phys), allocatable :: + real(kind=kind_phys), allocatable :: ! & plnmid(:,:), wtgwc(:,:), & plnmid(:,:), taugw(:,:), & utgwcl(:,:), vtgwcl(:,:), @@ -342,7 +342,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & ! ucltop : Zonal wind at the cloud top ( mid level ) ! vcltop : Meridional wind at the cloud top ( mid level ) ! windcltop : Wind speed at the cloud top ( mid level ) -! shear : Vertical shear of basic wind +! shear : Vertical shear of basic wind ! cosphi : Cosine of angle of wind vector at the cloud top ! sinphi : Sine of angle of wind vector at the cloud top ! c1 : Tunable parameter @@ -392,7 +392,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & vtgwc(i,k) = 0.0 ! brunm(i,k) = 0.0 ! rhom(i,k) = 0.0 - enddo + enddo enddo do i=1,im tauctx(i) = 0.0 @@ -468,7 +468,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & & rimin(npt,km+1), pint(npt,km+1)) ! allocate (ugwdc(npt,km), vgwdc(npt,km), - allocate + allocate ! & (plnmid(npt,km), wtgwc(npt,km), & (plnmid(npt,km), velco(npt,km), & utgwcl(npt,km), vtgwcl(npt,km), @@ -580,7 +580,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & ! 4 ======== pint(4) dpint(4) ! 4 -------- pmid(4) dpmid(4) ! ........ -! 17 ======== pint(17) dpint(17) +! 17 ======== pint(17) dpint(17) ! 17 -------- pmid(17) dpmid(17) ! 18 ======== pint(18) dpint(18) ! 18 -------- pmid(18) dpmid(18) @@ -645,11 +645,11 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & qtem = spfh(i,k-1) * tem1 + spfh(i,k) * tem2 rhoi(i,k) = pint(i,k) / ( rd * ti(i,k)*(1.0+fv*qtem) ) dtdp = (t(i,k)-t(i,k-1)) / (pmid(i,k)-pmid(i,k-1)) - n2 = gsqr / ti(i,k) * ( 1./cp - rhoi(i,k)*dtdp ) + n2 = gsqr / ti(i,k) * ( 1./cp - rhoi(i,k)*dtdp ) bruni(i,k) = sqrt (max (n2min, n2)) enddo enddo - + deallocate (spfh) !----------------------------------------------------------------------- ! @@ -660,7 +660,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & do k = 1, km do i = 1, npt dtdp = (ti(i,k+1)-ti(i,k)) / (pint(i,k+1)-pint(i,k)) - n2 = gsqr / t(i,k) * ( 1./cp - rhom(i,k)*dtdp ) + n2 = gsqr / t(i,k) * ( 1./cp - rhom(i,k)*dtdp ) brunm(i,k) = sqrt (max (n2min, n2)) enddo enddo @@ -804,7 +804,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & ! 18 -------- U(18) ! 19 ======== UI(19) rhoi(19) bruni(19) riloc(19) ! -!----------------------------------------------------------------------- +!----------------------------------------------------------------------- do k=2,km do i=1,npt @@ -816,10 +816,10 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & tem = bruni(i,k) / shear riloc(i,k) = tem * tem if (riloc(i,k) >= rimax ) riloc(i,k) = rilarge - end if + end if enddo enddo - + do i=1,npt riloc(i,1) = riloc(i,2) riloc(i,km+1) = riloc(i,km) @@ -860,33 +860,33 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & ! it can be thought that there is deep convective cloud. However, ! deep convective heating between kcldbot and kcldtop is sometimes ! zero in spite of kcldtop less than kcldbot. In this case, -! maximum deep convective heating is assumed to be 1.e-30. +! maximum deep convective heating is assumed to be 1.e-30. ! ! B : kk is the vertical index for interface level cloud top ! ! C : Total convective fractional cover (cldf) is used as the -! convective cloud cover for GWDC calculation instead of +! convective cloud cover for GWDC calculation instead of ! convective cloud cover in each layer (concld). ! a1 = cldf*dlength ! You can see the difference between cldf(i) and concld(i) -! in (4.a.2) in Description of the NCAR Community Climate +! in (4.a.2) in Description of the NCAR Community Climate ! Model (CCM3). ! In NCAR CCM3, cloud fractional cover in each layer in a deep ! cumulus convection is determined assuming total convective -! cloud cover is randomly overlapped in each layer in the +! cloud cover is randomly overlapped in each layer in the ! cumulus convection. ! ! D : Wave stress at cloud top is calculated when the atmosphere ! is dynamically stable at the cloud top ! -! E : Cloud top wave stress and nonlinear parameter are calculated +! E : Cloud top wave stress and nonlinear parameter are calculated ! using density, temperature, and wind that are defined at mid ! level just below the interface level in which cloud top wave ! stress is defined. ! Nonlinct is defined at the interface level. -! +! ! F : If the atmosphere is dynamically unstable at the cloud top, -! GWDC calculation in current horizontal grid is skipped. +! GWDC calculation in current horizontal grid is skipped. ! ! G : If mean wind at the cloud top is less than zero, GWDC ! calculation in current horizontal grid is skipped. @@ -894,12 +894,9 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & !----------------------------------------------------------------------- !D !> - Wave stress at cloud top is calculated when the atmosphere -!! is dynamically stable at the cloud top. - do i=1,npt - kk = kcldtop(i) - if ( abs(basicui(i,kk)) > zero .and. riloc(i,kk) > ricrit) then -!E -!> - The cloud top wave stress and nonlinear parameter are calculated +!! is dynamically stable at the cloud top +!! +!> - The cloud top wave stress and nonlinear parameter are calculated !! using density, temperature, and wind that are defined at mid !! level just below the interface level in which cloud top wave !! stress is defined. @@ -909,11 +906,11 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & !! \mu=\frac{gQ_{0}a_{1}}{c_{p}T_{0}NU^{2}} !! \f] !! where \f$Q_{0}\f$ is the maximum deep convective heating rate in a -!! horizontal grid point calculated from cumulus parameterization. +!! horizontal grid point calculated from cumulus parameterization. !! \f$a_{1}\f$ is the half-width of -!! the forcing function.\f$g\f$ is gravity. \f$c_{p}\f$ is specific -!! heat at constant pressure. \f$T_{0}\f$ is the layer mean -!! temperature (T1). As eqs.(18) and (19) \cite chun_and_baik_1998, +!! the forcing function.\f$g\f$ is gravity. \f$c_{p}\f$ is specific +!! heat at constant pressure. \f$T_{0}\f$ is the layer mean +!! temperature (T1). As eqs.(18) and (19) \cite chun_and_baik_1998, !! the zonal momentum flux is given by !! \f[ !! \tau_{x}=-[\rho U^{3}/(N\triangle x)]G(\mu) @@ -923,13 +920,26 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & !! G(\mu)=c_{1}c_2^2 \mu^{2} !! \f] !! wher \f$\rho\f$ is the local density. -!! The tunable parameter \f$c_1\f$ is related to the horizontal -!! structure of thermal forcing. The tunable parameter \f$c_2\f$ is -!! related to the basic-state wind and stability and the bottom and +!! The tunable parameter \f$c_1\f$ is related to the horizontal +!! structure of thermal forcing. The tunable parameter \f$c_2\f$ is +!! related to the basic-state wind and stability and the bottom and !! top heights of thermal forcing. If the atmosphere is dynamically -!! unstable at the cloud top, the convective GWD calculation is +!! unstable at the cloud top, the convective GWD calculation is !! skipped at that grid point. +!! +! - If mean wind at the cloud top is less than zero, GWDC +! calculation in current horizontal grid is skipped. +! +!> - The stress is capped at tauctmax = - 5\f$n/m^2\f$ +!! in order to prevent numerical instability. +! +!----------------------------------------------------------------------- +!D + do i=1,npt + kk = kcldtop(i) + if ( abs(basicui(i,kk)) > zero .and. riloc(i,kk) > ricrit) then +!E tem = basicum(i,kk) tem1 = tem * tem nonlinct = gqmcldlen(i) / (bruni(i,kk)*t(i,kk)*tem1) ! Mu @@ -956,7 +966,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & tauctxl(i) = zero tauctyl(i) = zero do_gwc(i) = .false. - end if + end if !H enddo @@ -983,11 +993,11 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & ! ! Minimum RI is calculated for the following two cases ! -! (1) RIloc < 1.e+20 +! (1) RIloc < 1.e+20 ! (2) Riloc = 1.e+20 ----> Vertically uniform basic-state wind ! ! RIloc cannot be smaller than zero because N^2 becomes 1.E-32 in the -! case of N^2 < 0.. Thus the sign of RINUM is determined by +! case of N^2 < 0.. Thus the sign of RINUM is determined by ! 1 - nonlin*|c2|. ! !----------------------------------------------------------------------- @@ -1058,20 +1068,20 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & ! 2 ======== - 0.001 -1.e20 ! 2 -------- 0.000 ! 3 ======== - 0.001 -1.e20 -! 3 -------- -.xxx +! 3 -------- -.xxx ! 4 ======== - 0.001 2.600 2.000 ! 4 -------- 0.000 ! 5 ======== - 0.001 2.500 2.000 ! 5 -------- 0.000 ! 6 ======== - 0.001 1.500 0.110 -! 6 -------- +.xxx +! 6 -------- +.xxx ! 7 ======== - 0.005 2.000 3.000 ! 7 -------- 0.000 ! 8 ======== - 0.005 1.000 0.222 ! 8 -------- +.xxx ! 9 ======== - 0.010 1.000 2.000 ! 9 -------- 0.000 -! kcldtopi 10 ======== $$$ - 0.010 +! kcldtopi 10 ======== $$$ - 0.010 ! kcldtop 10 -------- $$$ yyyyy ! 11 ======== $$$ 0 ! 11 -------- $$$ @@ -1086,16 +1096,16 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & ! 16 ======== $$$ 0 ! kcldbot 16 -------- $$$ ! 17 ======== 0 -! 17 -------- +! 17 -------- ! 18 ======== 0 -! 18 -------- +! 18 -------- ! 19 ======== 0 ! !----------------------------------------------------------------------- ! ! Even though the cloud top level obtained in deep convective para- ! meterization is defined in mid-level, the cloud top level for -! the GWDC calculation is assumed to be the interface level just +! the GWDC calculation is assumed to be the interface level just ! above the mid-level cloud top vertical level index. ! !----------------------------------------------------------------------- @@ -1118,22 +1128,22 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & taugwci(i,k) = taugwci(i,k+1) elseif (rimin(i,k) > riminp) then tem = 2.0 + 1.0 / sqrt(riloc(i,k)) - nonlins = (1.0/abs(c2)) * (2.*sqrt(tem) - tem) + nonlins = (1.0/abs(c2)) * (2.*sqrt(tem) - tem) tem1 = basicui(i,k) tem2 = c2*nonlins*tem1 taugwci(i,k) = - rhoi(i,k) * c1 * tem1 * tem2 * tem2 & / (bruni(i,k)*dlen(i)) elseif (rimin(i,k) > riminm) then - taugwci(i,k) = zero -! taugwci(i,k) = taugwci(i,k+1) + taugwci(i,k) = zero +! taugwci(i,k) = taugwci(i,k+1) end if ! RImin else -!> - If the minimum \f$R_{i}\f$ at interface cloud top is less than +!> - If the minimum \f$R_{i}\f$ at interface cloud top is less than !! or equal to 1/4, the convective GWD calculation is skipped at that !! grid point. - taugwci(i,k) = zero + taugwci(i,k) = zero end if ! RIloc else taugwci(i,k) = zero @@ -1179,7 +1189,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & taugw(i,k) = (taugwci(i,k+1) - taugwci(i,k)) / dpmid(i,k) if (taugw(i,k) /= 0.0) then tem = deltim * taugw(i,k) - dtfac(i) = min(dtfac(i), abs(velco(i,k)/tem)) + dtfac(i) = min(dtfac(i), abs(velco(i,k)/tem)) endif else taugw(i,k) = 0.0 @@ -1284,8 +1294,8 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & !----------------------------------------------------------------------- ! -! The GWDC should accelerate the zonal and meridional wind in the -! opposite direction of the previous zonal and meridional wind, +! The GWDC should accelerate the zonal and meridional wind in the +! opposite direction of the previous zonal and meridional wind, ! respectively ! !----------------------------------------------------------------------- @@ -1296,10 +1306,10 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & !-------------------- x-component------------------- -! write(6,'(a)') +! write(6,'(a)') ! + '(GWDC) WARNING: The GWDC should accelerate the zonal wind ' -! write(6,'(a,a,i3,a,i3)') -! + 'in the opposite direction of the previous zonal wind', +! write(6,'(a,a,i3,a,i3)') +! + 'in the opposite direction of the previous zonal wind', ! + ' at I = ',i,' and J = ',lat ! write(6,'(4(1x,e17.10))') u(i,kk),v(i,kk),u(i,k),v(i,k) ! write(6,'(a,1x,e17.10))') 'Vcld . V =', @@ -1309,7 +1319,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & ! do k1=1,km ! write(6,'(i2,36x,2(1x,e17.10))') ! + k1,taugwcxi(i,k1),taugwci(i,k1) -! write(6,'(i2,2(1x,e17.10))') k1,utgwcl(i,k1),u(i,k1) +! write(6,'(i2,2(1x,e17.10))') k1,utgwcl(i,k1),u(i,k1) ! end do ! write(6,'(i2,36x,1x,e17.10)') (km+1),taugwcxi(i,km+1) ! end if @@ -1319,7 +1329,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & ! do k1=1,km ! write(6,'(i2,36x,2(1x,e17.10))') ! + k1,taugwci(i,k1) -! write(6,'(i2,2(1x,e17.10))') k1,wtgwc(i,k1),basicum(i,k1) +! write(6,'(i2,2(1x,e17.10))') k1,wtgwc(i,k1),basicum(i,k1) ! end do ! write(6,'(i2,36x,1x,e17.10)') (km+1),taugwci(i,km+1) @@ -1339,7 +1349,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & ! do k1=1,km ! write(6,'(i2,36x,2(1x,e17.10))') ! + k1,taugwcyi(i,k1),taugwci(i,k1) -! write(6,'(i2,2(1x,e17.10))') k1,vtgwc(i,k1),v(i,k1) +! write(6,'(i2,2(1x,e17.10))') k1,vtgwc(i,k1),v(i,k1) ! end do ! write(6,'(i2,36x,1x,e17.10)') (km+1),taugwcyi(i,km+1) ! end if @@ -1369,7 +1379,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & !----------------------------------------------------------------------- ! -! For GWDC performance analysis +! For GWDC performance analysis ! !----------------------------------------------------------------------- @@ -1385,7 +1395,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & ! if ( abs(taugwci(i,k)-taugwci(i,kk)) > taumin ) then ! break(i) = 1.0 ! go to 2000 -! endif +! endif ! enddo !2000 continue diff --git a/physics/gwdps.f b/physics/gwdps.f index 62ef1b4ce..337111633 100644 --- a/physics/gwdps.f +++ b/physics/gwdps.f @@ -178,12 +178,169 @@ end subroutine gwdps_init !! | me | mpi_rank | rank of the current MPI task | index | 0 | integer | | in | F | !! | lprnt | flag_print | flag for debugging printouts | flag | 0 | logical | | in | F | !! | ipr | horizontal_index_of_printed_column | horizontal index of column used in debugging printouts | index | 0 | integer | | in | F | +!! | rdxzb | level_of_dividing_streamline | level of the dividing streamline | none | 1 | real | kind_phys | out | F | !! | errmsg | error_message | error message for error handling in CCPP | none | 0 | character | len=* | out | F | !! | errflg | error_flag | error flag for error handling in CCPP | flag | 0 | integer | | out | F | !! !> \section gen_gwdps GFS Orographic GWD Scheme General Algorithm !! -# Calculate subgrid mountain blocking !! -# Calculate orographic wave drag +!! +!! The NWP model gravity wave drag (GWD) scheme in the GFS has two +!! main components: how the surface stress is computed, and then how +!! that stress is distributed over a vertical column where it may +!! interact with the models momentum. Each of these depends on the +!! large scale environmental atmospheric state and assumptions about +!! the sub-grid scale processes. In Alpert GWD (1987) based on linear, +!! two-dimensional non-rotating, stably stratified flow over a mountain ridge, +!! sub-grid scale gravity wave motions are assumed which propagate away +!! from the mountain. Described in Alpert (1987), the flux measured over +!! a "low level" vertically averaged layer, in the atmosphere defines a base +!! level flux. "Low level" was taken to be the first 1/3 of the troposphere +!! in the 1987 implementation. This choice was meant to encompass a thick +!! low layer for vertical averages of the environmental (large scale) flow +!! quantities. The vertical momentum flux or gravity wave stress in a +!! grid box due to a single mountain is given as in Pierrehumbert, (1987) (PH): +!! +!! \f$ \tau = \frac {\rho \: U^{3}\: G(F_{r})} {\Delta X \; N } \f$ +!! +!! emetic \f$ \Delta X \f$ is a grid increment, N is the Brunt Viasala frequency +!! +!! +!! \f$ N(\sigma) = \frac{-g \: \sigma \: +!! \frac{\partial\Theta}{\partial\sigma}}{\Theta \:R \:T} \f$ +!! +!! The environmental variables are calculated from a mass weighted vertical +!! average over a base layer. G(Fr) is a monotonically increasing +!! function of Froude number, +!! +!! \f$ F_{r} = \frac{N h^{'}}{U} \f$ +!! +!! where U is the wind speed calculated as a mass weighted vertical average in +!! the base layer, and h', is the vertical displacement caused by the orography +!! variance. An effective mountain length for the gravity wave processes, +!! +!! \f$ l^{*} = \frac{\Delta X}{m} \f$ +!! +!! where m is the number of mountains in a grid box, can then +!! be defined to obtain the form of the base level stress +!! +!! +!! \f$ \tau = \frac {\rho \: U^{3} \: G(F_{r})} {N \;l^{*}} \f$ +!! +!! giving the stress induced from the surface in a model grid box. +!! PH gives the form for the function G(Fr) as +!! +!! +!! \f$ G(F_{r}) = \bar{G}\frac{F^{2}_{r}}{F^{2}_{r}\: + \:a^{2}} \f$ +!! +!! Where \f$ \bar{G} \f$ is an order unity non-dimensional saturation +!! flux set to 1 and 'a' is a function of the mountain aspect ratio also +!!set to 1 in the 1987 implementation of the GFS GWD. Typical values of +!! U=10m/s, N=0.01 1/s, l*=100km, and a=1, gives a flux of 1 Pascal and +!! if this flux is made to go to zero linearly with height then the +!! decelerations would be about 10/m/s/day which is consistent with +!! observations in PH. +!! +!! +!! In Kim, Moorthi, Alpert's (1998, 2001) GWD currently in GFS operations, +!! the GWD scheme has the same physical basis as in Alpert (1987) with the addition +!! of enhancement factors for the amplitude, G, and mountain shape details +!! in G(Fr) to account for effects from the mountain blocking. A factor, +!! E m’, is an enhancement factor on the stress in the Alpert '87 scheme. +!! The E ranges from no enhancement to an upper limit of 3, E=E(OA)[1-3], +!! and is a function of OA, the Orographic Asymmetry defined in KA (1995) as +!! +!! Orographic Asymmetry (OA) = \f$ \frac{ \bar{x} \; - \; +!! \sum\limits_{j=1}^{N_{b}} x_{j} \; n_{j} }{\sigma_{x}} \f$ +!! +!! where Nb is the total number of bottom blocks in the mountain barrier, +!! \f$ \sigma_{x} \f$ is the standard deviation of the horizontal distance defined by +!! +!! \f$ \sigma_{x} = \sqrt{ \frac{\sum\limits_{j=1}^{N_{b}} +!! \; (x_{j} \; - \; \bar{x} )^2}{N_{x}} } \f$ +!! +!! +!! where Nx is the number of grid intervals for the large scale domain being +!! considered. So the term, E(OA)m’/ \f$ \Delta X \f$ in Kim's scheme represents +!! a multiplier on G shown in Alpert's eq (1), where m’ is the number of mountains +!! in a sub-grid scale box. Kim increased the complexity of m’ making it a +!! function of the fractional area of the sub-grid mountain and the asymmetry +!! and convexity statistics which are found from running a gravity wave +!! model for a large number of cases: +!! +!! \f$ m^{'} = C_{m} \Delta X \left[ \frac{1 \; + \; +!! \sum\limits_{x} L_{h} }{\Delta X} \right]^{OA+1} \f$ +!! +!! Where, according to Kim, \f$ \sum \frac{L_{h}}{\Delta X} \f$ is +!! the fractional area covered by the subgrid-scale orography higher than +!! a critical height \f$ h_{c} = Fr_{c} U_{0}/N_{0} \f$ , over the +!! "low level" vertically averaged layer, for a grid box with the interval +!! \f$ \Delta X \f$. Each \f$ L_{n}\f$ is the width of a segment of +!! orography intersection at the critical height: +!! +!! \f$ Fr_{0} = \frac{N_{0} \; h^{'}}{U_{0}} \f$ +!! +!! \f$ G^{'}(OC,Fr_{0}) = \frac{Fr_{0}^{2}}{Fr_{0}^{2} \; + \; a^{2}} \f$ +!! +!! \f$ a^{2} = \frac{C_{G}}{OC} \f$ +!! +!! \f$ E(OA, Fr_{0}) = (OA \; + \; 2)^{\delta} \f$ and \f$ \delta +!! \; = \; \frac{C_{E} \; Fr_{0}}{Fr_{c}} \f$ where \f$ Fr_{c} \f$ +!! is as in Alpert. +!! +!! +!! This represents a closed scheme, somewhat empirical adjustments +!! to the original scheme to calculate the surface stress. +!! +!! Momentum is deposited by the sub-grid scale gravity waves break due +!! to the presence of convective mixing assumed to occur when the +!! minimum Richardson number: +!! +!! Orographic Convexity (OC) = \f$ \frac{ \sum\limits_{j=1}^{N_{x}} +!! \; (h_{j} \; - \; \bar{h})^4 }{N_{x} \;\sigma_{h}^4} \f$ , +!! and where \f$ \sigma_{h} = \sqrt{ \frac{\sum\limits_{j=1}^{N_{x}} +!! \; (h_{j} \; - \; \bar{h} )^2}{N_{x}} } \f$ +!! +!! This represents a closed scheme, somewhat empirical adjustments +!! to the original scheme to calculate the surface stress. +!! +!! Momentum is deposited by the sub-grid scale gravity waves break due +!! to the presence of convective mixing assumed to occur when +!! the minimum Richardson number: +!! +!! \f$ Ri_{m} = \frac{Ri(1 \; - \; Fr)}{(1 \; + \; \sqrt{Ri}Fr)^2} \f$ +!! +!! Is less than 1/4 Or if critical layers are encountered in a layer +!! the the momentum flux will vanish. The critical layer is defined +!! when the base layer wind becomes perpendicular to the environmental +!! wind. Otherwise, wave breaking occurs at a level where the amplification +!! of the wave causes the local Froude number or similarly a truncated +!! (first term of the) Scorer parameter, to be reduced below a critical +!! value by the saturation hypothesis (Lindzen,). This is done through +!! eq 1 which can be written as +!! +!! \f$ \tau = \rho U N k h^{'2} \f$ +!! +!! For small Froude number this is discretized in the vertical so at each +!! level the stress is reduced by ratio of the Froude or truncated Scorer +!! parameter, \f$ \frac{U^{2}}{N^{2}} = \frac{N \tau_{l-1}}{\rho U^{3} k} \f$ , +!! where the stress is from the layer below beginning with that found near +!! the surface. The respective change in momentum is applied in +!! that layer building up from below. +!! +!! An amplitude factor is part of the calibration of this scheme which is +!! a function of the model resolution and the vertical diffusion. This +!! is because the vertical diffusion and the GWD account encompass +!! similar physical processes. Thus, one needs to run the model over +!! and over for various amplitude factors for GWD and vertical diffusion. +!! +!! In addition, there is also mountain blocking from lift and frictional +!! forces. Improved integration between how the GWD is calculated and +!! the mountain blocking of wind flow around sub-grid scale orography +!! is underway at NCEP. The GFS already has convectively forced GWD +!! an independent process. The next step is to test +!! !> \section det_gwdps GFS Orographic GWD Scheme Detailed Algorithm !> @{ subroutine gwdps_run( & @@ -191,7 +348,7 @@ subroutine gwdps_run( & & PRSI,DEL,PRSL,PRSLK,PHII, PHIL,DELTIM,KDT, & & HPRIME,OC,OA4,CLX4,THETA,SIGMA,GAMMA,ELVMAX, & & DUSFC,DVSFC,G, CP, RD, RV, IMX, & - & nmtvr, cdmbgwd, me, lprnt, ipr, errmsg, errflg) + & nmtvr, cdmbgwd, me, lprnt, ipr, rdxzb, errmsg, errflg) ! ! ******************************************************************** ! -----> I M P L E M E N T A T I O N V E R S I O N <---------- @@ -308,7 +465,8 @@ subroutine gwdps_run( & real(kind=kind_phys), intent(inout) :: ELVMAX(IM) real(kind=kind_phys), intent(in) :: & & THETA(IM), SIGMA(IM), GAMMA(IM) - real(kind=kind_phys), intent(out) :: DUSFC(IM), DVSFC(IM) + real(kind=kind_phys), intent(out) :: DUSFC(IM), DVSFC(IM), & + & RDXZB(IX) integer, intent(in) :: nmtvr logical, intent(in) :: lprnt character(len=*), intent(out) :: errmsg @@ -345,19 +503,19 @@ subroutine gwdps_run( & ! real(kind=kind_phys) dpmin,hminmt,hncrit,minwnd,sigfac ! --- for lm mtn blocking -! parameter (cdmb = 1.0) ! non-dim sub grid mtn drag Amp (*j*) - parameter (hncrit=8000.) ! Max value in meters for ELVMAX (*j*) +! parameter (cdmb = 1.0) !< non-dim sub grid mtn drag Amp (*j*) + parameter (hncrit=8000.) !< Max value in meters for ELVMAX (*j*) ! hncrit set to 8000m and sigfac added to enhance elvmax mtn hgt - parameter (sigfac=4.0) ! MB3a expt test for ELVMAX factor (*j*) - parameter (hminmt=50.) ! min mtn height (*j*) - parameter (minwnd=0.1) ! min wind component (*j*) + parameter (sigfac=4.0) !< MB3a expt test for ELVMAX factor (*j*) + parameter (hminmt=50.) !< min mtn height (*j*) + parameter (minwnd=0.1) !< min wind component (*j*) -! parameter (dpmin=00.0) ! Minimum thickness of the reference layer -!! parameter (dpmin=05.0) ! Minimum thickness of the reference layer -! parameter (dpmin=20.0) ! Minimum thickness of the reference layer - ! in centibars - parameter (dpmin=5000.0) ! Minimum thickness of the reference layer - ! in Pa +! parameter (dpmin=00.0) !< Minimum thickness of the reference layer +!! parameter (dpmin=05.0) !< Minimum thickness of the reference layer +! parameter (dpmin=20.0) !< Minimum thickness of the reference layer + !< in centibars + parameter (dpmin=5000.0) !< Minimum thickness of the reference layer + !< in Pa ! real(kind=kind_phys) FDIR integer mdir @@ -372,8 +530,8 @@ subroutine gwdps_run( & ! real(kind=kind_phys) TAUB(IM), XN(IM), YN(IM), UBAR(IM) & &, VBAR(IM), ULOW(IM), OA(IM), CLX(IM) & - &, ROLL(IM), ULOI(IM), DTFAC(IM), XLINV(IM) & - &, DELKS(IM), DELKS1(IM) + &, ROLL(IM), ULOI(IM) & + &, DTFAC(IM), XLINV(IM), DELKS(IM), DELKS1(IM) ! real(kind=kind_phys) BNV2(IM,KM), TAUP(IM,KM+1), ri_n(IM,KM) & &, TAUD(IM,KM), RO(IM,KM), VTK(IM,KM) & @@ -436,6 +594,7 @@ subroutine gwdps_run( & ! IF ( NMTVR .eq. 14) then ! ---- for lm and gwd calculation points + RDXZB(:) = 0. ipt = 0 npt = 0 DO I = 1,IM @@ -613,7 +772,10 @@ subroutine gwdps_run( & EK(I) = 0.5 * UP(I) * UP(I) ! --- Dividing Stream lime is found when PE =exceeds EK. - IF ( PE(I) .ge. EK(I) ) IDXZB(I) = K + IF ( PE(I) .ge. EK(I) ) THEN + IDXZB(I) = K + RDXZB(J) = real(K,kind=kind_phys) + ENDIF ! --- Then mtn blocked flow is between Zb=k(IDXZB(I)) and surface ! !> - The dividing streamline height (idxzb), of a subgrid scale @@ -737,6 +899,7 @@ subroutine gwdps_run( & ! do i=1,npt IDXZB(i) = 0 + RDXZB(i) = 0. enddo ENDIF ! @@ -1194,8 +1357,8 @@ subroutine gwdps_run( & ! if ( ABS(DBIM * U1(J,K)) .gt. .01 ) ! & print *,' in gwdps_lmi.f KDT=',KDT,I,K,DB(I,K), ! & dbim,idxzb(I),U1(J,K),V1(J,K),me - DUSFC(J) = DUSFC(J) - DBIM * V1(J,K) * DEL(J,K) - DVSFC(J) = DVSFC(J) - DBIM * U1(J,K) * DEL(J,K) + DUSFC(J) = DUSFC(J) - DBIM * U1(J,K) * DEL(J,K) + DVSFC(J) = DVSFC(J) - DBIM * V1(J,K) * DEL(J,K) else ! A(J,K) = DTAUY + A(J,K) From 1867fc3e546739144279490d47e23fcb36fd551f Mon Sep 17 00:00:00 2001 From: "Man.Zhang" Date: Sun, 24 Jun 2018 21:54:00 -0600 Subject: [PATCH 2/6] gwdps_run pass ccpp rt tests with option B. --- physics/gwdps.f | 156 ++++++++++++++++++++++++------------------------ 1 file changed, 78 insertions(+), 78 deletions(-) diff --git a/physics/gwdps.f b/physics/gwdps.f index 337111633..d5f9be09d 100644 --- a/physics/gwdps.f +++ b/physics/gwdps.f @@ -131,7 +131,7 @@ end subroutine gwdps_init !! \brief This subroutine includes orographic gravity wave drag and mountain !! blocking. !! -!! The time tendencies of zonal and meridional wind are altered to +!> The time tendencies of zonal and meridional wind are altered to !! include the effect of mountain induced gravity wave drag from !! subgrid scale orography including convective breaking, shear !! breaking and the presence of critical levels. @@ -361,7 +361,7 @@ subroutine gwdps_run( & !----- ALSO INCLUDED IS RI SMOOTH OVER A THICK LOWER LAYER !----- ALSO INCLUDED IS DECREASE IN DE-ACC AT TOP BY 1/2 !----- THE NMC GWD INCORPORATING BOTH GLAS(P&S) AND GFDL(MIGWD) -!----- MOUNTAIN INDUCED GRAVITY WAVE DRAG +!----- MOUNTAIN INDUCED GRAVITY WAVE DRAG !----- CODE FROM .FR30(V3MONNX) FOR MONIN3 !----- THIS VERSION (06 MAR 1987) !----- THIS VERSION (26 APR 1987) 3.G @@ -437,8 +437,8 @@ subroutine gwdps_run( & ! OTHER INPUT VARIABLES UNMODIFIED. ! revision log: ! May 2013 J. Wang change cleff back to opn setting -! Jan 2014 J. Wang merge Henry and Fangin's dissipation heat in gfs to nems -! +! Jan 2014 J. Wang merge Henry and Fangin's dissipation heat in gfs to nems +! ! ! ******************************************************************** USE MACHINE , ONLY : kind_phys @@ -497,8 +497,8 @@ subroutine gwdps_run( & parameter (FRC=1.0, CE=0.8, CEOFRC=CE/FRC, frmax=100., CG=0.5) parameter (GMAX=1.0, VELEPS=1.0, FACTOP=0.5) ! parameter (GMAX=1.0, CRITAC=5.0E-4, VELEPS=1.0, FACTOP=0.5) - parameter (RLOLEV=50000.0) -! parameter (RLOLEV=500.0) + parameter (RLOLEV=50000.0) +! parameter (RLOLEV=500.0) ! parameter (RLOLEV=0.5) ! real(kind=kind_phys) dpmin,hminmt,hncrit,minwnd,sigfac @@ -592,13 +592,13 @@ subroutine gwdps_run( & LCAPP1 = LCAP + 1 ! ! - IF ( NMTVR .eq. 14) then + IF ( NMTVR .eq. 14) then ! ---- for lm and gwd calculation points - RDXZB(:) = 0. + RDXZB(:) = 0. ! CCPP: it is a bug ipt = 0 npt = 0 DO I = 1,IM - IF ( (elvmax(i) .GT. HMINMT) + IF ( (elvmax(i) .GT. HMINMT) & .and. (hprime(i) .GT. hpmin) ) then npt = npt + 1 ipt(npt) = i @@ -616,7 +616,7 @@ subroutine gwdps_run( & ! do i=1,npt iwklm(i) = 2 - IDXZB(i) = 0 + IDXZB(i) = 0 kreflm(i) = 0 enddo ! if (lprnt) @@ -632,7 +632,7 @@ subroutine gwdps_run( & ! then do not need hncrit -- test with large hncrit first. ! KMLL = km / 2 ! maximum mtnlm height : # of vertical levels / 2 KMLL = kmm1 -! --- No mtn should be as high as KMLL (so we do not have to start at +! --- No mtn should be as high as KMLL (so we do not have to start at ! --- the top of the model but could do calc for all levels). ! DO I = 1, npt @@ -648,12 +648,12 @@ subroutine gwdps_run( & pkp1log = phil(j,k+1) / G pklog = phil(j,k) / G !!!------- ELVMAX(J) = min (ELVMAX(J) + sigfac * hprime(j), hncrit) - if ( ( ELVMAX(j) .le. pkp1log ) .and. + if ( ( ELVMAX(j) .le. pkp1log ) .and. & ( ELVMAX(j) .ge. pklog ) ) THEN ! print *,' in gwdps_lm.f 1 =',k,ELVMAX(j),pklog,pkp1log,me -! --- wk for diags but can be saved and reused. +! --- wk for diags but can be saved and reused. wk(i) = G * ELVMAX(j) / ( phil(j,k+1) - phil(j,k) ) - iwklm(I) = MAX(iwklm(I), k+1 ) + iwklm(I) = MAX(iwklm(I), k+1 ) ! print *,' in gwdps_lm.f 2 npt=',npt,i,j,wk(i),iwklm(i),me endif ! @@ -671,16 +671,16 @@ subroutine gwdps_run( & ! jhit = 0 ! do i = 1, npt ! j=ipt(i) -! if ( iwklm(i) .gt. ihit ) then +! if ( iwklm(i) .gt. ihit ) then ! ihit = iwklm(i) ! jhit = j ! endif ! enddo ! print *, ' mb: kdt,max(iwklm),jhit,phil,me=', ! & kdt,ihit,jhit,phil(jhit,ihit),me - + klevm1 = KMLL - 1 - DO K = 1, klevm1 + DO K = 1, klevm1 DO I = 1, npt j = ipt(i) RDZ = g / ( phil(j,k+1) - phil(j,k) ) @@ -705,7 +705,7 @@ subroutine gwdps_run( & BNV2bar(I) = (PRSL(J,1)-PRSL(J,2)) * DELKS1(I) * BNV2LM(I,1) ENDDO -! --- find the dividing stream line height +! --- find the dividing stream line height ! --- starting from the level above the max mtn downward ! --- iwklm(i) is the k-index of mtn elvmax elevation !> - Find the dividing streamline height starting from the level above @@ -723,14 +723,14 @@ subroutine gwdps_run( & ! --- make averages, guess dividing stream (DS) line layer. ! --- This is not used in the first cut except for testing and ! --- is the vert ave of quantities from the surface to mtn top. -! +! DO I = 1, npt DO K = 1, Kreflm(I) J = ipt(i) RDELKS = DEL(J,K) * DELKS(I) - UBAR(I) = UBAR(I) + RDELKS * U1(J,K) ! trial Mean U below - VBAR(I) = VBAR(I) + RDELKS * V1(J,K) ! trial Mean V below - ROLL(I) = ROLL(I) + RDELKS * RO(I,K) ! trial Mean RO below + UBAR(I) = UBAR(I) + RDELKS * U1(J,K) ! trial Mean U below + VBAR(I) = VBAR(I) + RDELKS * V1(J,K) ! trial Mean V below + ROLL(I) = ROLL(I) + RDELKS * RO(I,K) ! trial Mean RO below RDELKS = (PRSL(J,K)-PRSL(J,K+1)) * DELKS1(I) BNV2bar(I) = BNV2bar(I) + BNV2lm(I,K) * RDELKS ! --- these vert ave are for diags, testing and GWD to follow (*j*). @@ -739,7 +739,7 @@ subroutine gwdps_run( & ! print *,' in gwdps_lm.f 5 =',i,kreflm(npt),BNV2bar(npt),me ! ! --- integrate to get PE in the trial layer. -! --- Need the first layer where PE>EK - as soon as +! --- Need the first layer where PE>EK - as soon as ! --- IDXZB is not 0 we have a hit and Zb is found. ! DO I = 1, npt @@ -755,21 +755,21 @@ subroutine gwdps_run( & !!\f[ !! UDS=\max(\sqrt{U1^2+V1^2},minwnd) !!\f] -!! where \f$ minwnd=0.1 \f$, \f$U1\f$ and \f$V1\f$ are zonal and +!! where \f$ minwnd=0.1 \f$, \f$U1\f$ and \f$V1\f$ are zonal and !! meridional wind components of model layer wind. - UDS(I,K) = + UDS(I,K) = & MAX(SQRT(U1(J,K)*U1(J,K) + V1(J,K)*V1(J,K)), minwnd) ! --- Test to see if we found Zb previously IF (IDXZB(I) .eq. 0 ) then - PE(I) = PE(I) + BNV2lm(I,K) * - & ( G * ELVMAX(J) - phil(J,K) ) * + PE(I) = PE(I) + BNV2lm(I,K) * + & ( G * ELVMAX(J) - phil(J,K) ) * & ( PHII(J,K+1) - PHII(J,K) ) / (G*G) ! --- KE ! --- Wind projected on the line perpendicular to mtn range, U(Zb(K)). ! --- kenetic energy is at the layer Zb ! --- THETA ranges from -+90deg |_ to the mtn "largest topo variations" UP(I) = UDS(I,K) * cos(ANG(I,K)) - EK(I) = 0.5 * UP(I) * UP(I) + EK(I) = 0.5 * UP(I) * UP(I) ! --- Dividing Stream lime is found when PE =exceeds EK. IF ( PE(I) .ge. EK(I) ) THEN @@ -778,18 +778,18 @@ subroutine gwdps_run( & ENDIF ! --- Then mtn blocked flow is between Zb=k(IDXZB(I)) and surface ! -!> - The dividing streamline height (idxzb), of a subgrid scale -!! obstable, is found by comparing the potential (PE) and kinetic +!> - The dividing streamline height (idxzb), of a subgrid scale +!! obstable, is found by comparing the potential (PE) and kinetic !! energies (EK) of the upstream large scale wind and subgrid scale air !! parcel movements. the dividing streamline is found when -!! \f$PE\geq EK\f$. Mountain-blocked flow is defined to exist between -!! the surface and the dividing streamline height (\f$h_d\f$), which -!! can be found by solving an integral equation for \f$h_d\f$: +!! \f$PE\geq EK\f$. Mountain-blocked flow is defined to exist between +!! the surface and the dividing streamline height (\f$h_d\f$), which +!! can be found by solving an integral equation for \f$h_d\f$: !!\f[ !! \frac{U^{2}(h_{d})}{2}=\int_{h_{d}}^{H} N^{2}(z)(H-z)dz !!\f] !! where \f$H\f$ is the maximum subgrid scale elevation within the grid -!! box of actual orography, \f$h\f$, obtained from the GTOPO30 dataset +!! box of actual orography, \f$h\f$, obtained from the GTOPO30 dataset !! from the U.S. Geological Survey. ENDIF ENDDO @@ -818,12 +818,12 @@ subroutine gwdps_run( & ! endif ! ! --- The drag for mtn blocked flow -! +! DO I = 1, npt J = ipt(i) ZLEN = 0. ! print *,' in gwdps_lm.f 9 =',i,j,IDXZB(i),me - IF ( IDXZB(I) .gt. 0 ) then + IF ( IDXZB(I) .gt. 0 ) then DO K = IDXZB(I), 1, -1 IF ( PHIL(J,IDXZB(I)) .gt. PHIL(J,K) ) then @@ -832,9 +832,9 @@ subroutine gwdps_run( & !!\f[ !! ZLEN=\sqrt{[\frac{h_{d}-z}{z+h'}]} !!\f] -!! where \f$z\f$ is the height, \f$h'\f$ is the orographic standard +!! where \f$z\f$ is the height, \f$h'\f$ is the orographic standard !! deviation (HPRIME). - ZLEN = SQRT( ( PHIL(J,IDXZB(I)) - PHIL(J,K) ) / + ZLEN = SQRT( ( PHIL(J,IDXZB(I)) - PHIL(J,K) ) / & ( PHIL(J,K ) + G * hprime(J) ) ) ! --- lm eq 14: !> - Calculate the drag coefficient to vary with the aspect ratio of @@ -843,29 +843,29 @@ subroutine gwdps_run( & !!\f[ !! R=\frac{\cos^{2}\psi+\gamma\sin^{2}\psi}{\gamma\cos^{2}\psi+\sin^{2}\psi} !!\f] -!! where \f$\psi\f$, which is derived from THETA, is the angle between -!! the incident flow direction and the normal ridge direcion. +!! where \f$\psi\f$, which is derived from THETA, is the angle between +!! the incident flow direction and the normal ridge direcion. !! \f$\gamma\f$ is the orographic anisotropy (GAMMA). - R = (cos(ANG(I,K))**2 + GAMMA(J) * sin(ANG(I,K))**2) / + R = (cos(ANG(I,K))**2 + GAMMA(J) * sin(ANG(I,K))**2) / & (gamma(J) * cos(ANG(I,K))**2 + sin(ANG(I,K))**2) ! --- (negitive of DB -- see sign at tendency) -!> - In each model layer below the dividing streamlines, a drag from +!> - In each model layer below the dividing streamlines, a drag from !! the blocked flow is exerted by the obstacle on the large scale flow. !! The drag per unit area and per unit height is written (eq.15 in !! \cite lott_and_miller_1997): !!\f[ !! D_{b}(z)=-C_{d}\max(2-\frac{1}{R},0)\rho\frac{\sigma}{2h'}ZLEN\max(\cos\psi,\gamma\sin\psi)\frac{UDS}{2} !!\f] -!! where \f$C_{d}\f$ is a specified constant, \f$\sigma\f$ is the -!! orographic slope. +!! where \f$C_{d}\f$ is a specified constant, \f$\sigma\f$ is the +!! orographic slope. DBTMP = 0.25 * CDmb * & MAX( 2. - 1. / R, 0. ) * sigma(J) * & MAX(cos(ANG(I,K)), gamma(J)*sin(ANG(I,K))) * - & ZLEN / hprime(J) - DB(I,K) = DBTMP * UDS(I,K) + & ZLEN / hprime(J) + DB(I,K) = DBTMP * UDS(I,K) ! -! if(lprnt .and. i .eq. npr) then +! if(lprnt .and. i .eq. npr) then ! print *,' in gwdps_lmi.f 10 npt=',npt,i,j,idxzb(i) ! &, DBTMP,R' ang=',ang(i,k),' gamma=',gamma(j),' K=',K ! print *,' in gwdps_lmi.f 11 K=',k,ZLEN,cos(ANG(I,K)) @@ -876,13 +876,13 @@ subroutine gwdps_run( & ! if(lprnt) print *,' @K=1,ZLEN,DBTMP=',K,ZLEN,DBTMP endif ENDDO -! +! !............................. !............................. ! end mtn blocking section ! - ELSEIF ( NMTVR .ne. 14) then -! ---- for mb not present and gwd (nmtvr .ne .14) + ELSEIF ( NMTVR .ne. 14) then +! ---- for mb not present and gwd (nmtvr .ne .14) ipt = 0 npt = 0 DO I = 1,IM @@ -982,13 +982,13 @@ subroutine gwdps_run( & enddo enddo ! -!> - Calculate the reference level index: kref=max(2,KPBL+1). where +!> - Calculate the reference level index: kref=max(2,KPBL+1). where !! KPBL is the index for the PBL top layer. KBPS = 1 KMPS = KM DO I=1,npt J = ipt(i) - kref(I) = MAX(IWK(I), KPBL(J)+1 ) ! reference level + kref(I) = MAX(IWK(I), KPBL(J)+1 ) ! reference level DELKS(I) = 1.0 / (PRSI(J,1) - PRSI(J,kref(I))) DELKS1(I) = 1.0 / (PRSL(J,1) - PRSL(J,kref(I))) UBAR (I) = 0.0 @@ -1023,8 +1023,8 @@ subroutine gwdps_run( & ! NWD 1 2 3 4 5 6 7 8 ! WD W S SW NW E N NE SE ! -!> - Calculate low-level horizontal wind direction, the derived -!! orographic asymmetry parameter (OA), and the derived Lx (CLX). +!> - Calculate low-level horizontal wind direction, the derived +!! orographic asymmetry parameter (OA), and the derived Lx (CLX). DO I = 1,npt J = ipt(i) wdir = atan2(UBAR(I),VBAR(I)) + pi @@ -1053,7 +1053,7 @@ subroutine gwdps_run( & ULOW (I) = 0.0 DTFAC(I) = 1.0 ICRILV(I) = .FALSE. ! INITIALIZE CRITICAL LEVEL CONTROL VECTOR - + ! !----COMPUTE THE "LOW LEVEL" WIND MAGNITUDE (M/S) ! @@ -1072,7 +1072,7 @@ subroutine gwdps_run( & ! ENDIF ENDDO ENDDO -! +! ! ! find the interface level of the projected wind where ! low levels & upper levels meet above pbl @@ -1110,7 +1110,7 @@ subroutine gwdps_run( & ! RATIO CONST. USE SIMPLIFIED RELATIONSHIP BETWEEN STANDARD ! DEVIATION & CRITICAL HGT ! -!> - Calculate enhancement factor (E),number of mountans (m') and +!> - Calculate enhancement factor (E),number of mountans (m') and !! aspect ratio constant. !!\n As in eq.(4.9),(4.10),(4.11) in !! \cite kim_and_arakawa_1995, we define m' and E in such a way that they @@ -1135,7 +1135,7 @@ subroutine gwdps_run( & !!\f[ !! a^{2}=C_{G}OC^{-1} !!\f] -!! where \f$F_{r_{c}}(=1)\f$ is the critical Froude number, +!! where \f$F_{r_{c}}(=1)\f$ is the critical Froude number, !! \f$F_{r_{0}}\f$ is the Froude number. \f$C_{E}\f$,\f$C_{m}\f$, !! \f$C_{G}\f$ are constants. @@ -1144,8 +1144,8 @@ subroutine gwdps_run( & !!\f[ !! \tau_0=E\frac{m'}{\triangle x}\frac{\rho_{0}U_0^3}{N_{0}}G' !!\f] -!! where \f$E\f$,\f$m'\f$, and \f$G'\f$ are the enhancement factor, -!! "the number of mountains", and the flux function defined above, +!! where \f$E\f$,\f$m'\f$, and \f$G'\f$ are the enhancement factor, +!! "the number of mountains", and the flux function defined above, !! respectively. EFACT = (OA(I) + 2.) ** (CEOFRC*FR) @@ -1169,7 +1169,7 @@ subroutine gwdps_run( & SCOR(I) = BNV2(I,K) / TEM ! Scorer parameter below ref level ENDDO ! if(lprnt) print *,' taub=',taub -! +! !----SET UP BOTTOM VALUES OF STRESS ! DO K = 1, KBPS @@ -1215,7 +1215,7 @@ subroutine gwdps_run( & SCORK = BNV2(I,K) * TEMV * TEMV RSCOR = MIN(1.0, SCORK / SCOR(I)) SCOR(I) = SCORK - ELSE + ELSE RSCOR = 1. ENDIF ! @@ -1224,11 +1224,11 @@ subroutine gwdps_run( & !! \tau=\frac{m'}{\triangle x}\rho NUh_d^2 !!\f] !! where \f$h_{d}\f$ is the displacement wave amplitude. In the absence -!! of wave breaking, the displacement amplitude for the \f$i^{th}\f$ -!! layer can be expressed using the drag for the layer immediately +!! of wave breaking, the displacement amplitude for the \f$i^{th}\f$ +!! layer can be expressed using the drag for the layer immediately !! below. Thus, assuming \f$\tau_i=\tau_{i+1}\f$, we can get: !!\f[ -!! h_{d_i}^2=\frac{\triangle x}{m'}\frac{\tau_{i+1}}{\rho_{i}N_{i}U_{i}} +!! h_{d_i}^2=\frac{\triangle x}{m'}\frac{\tau_{i+1}}{\rho_{i}N_{i}U_{i}} !!\f] BRVF = SQRT(BNV2(I,K)) ! Brunt-Vaisala Frequency @@ -1240,9 +1240,9 @@ subroutine gwdps_run( & ! ! RIM is the MINIMUM-RICHARDSON NUMBER BY SHUTTS (1985) ! -!> - The minimum Richardson number (\f$Ri_{m}\f$) or local +!> - The minimum Richardson number (\f$Ri_{m}\f$) or local !! wave-modified Richardson number, which determines the onset of wave -!! breaking, is expressed in terms of \f$R_{i}\f$ and +!! breaking, is expressed in terms of \f$R_{i}\f$ and !! \f$F_{r_{d}}=Nh_{d}/U\f$: !!\f[ !! Ri_{m}=\frac{Ri(1-Fr_{d})}{(1+\sqrt{Ri}\cdot Fr_{d})^{2}} @@ -1259,17 +1259,17 @@ subroutine gwdps_run( & !> - Check stability to employ the 'saturation hypothesis' of !! \cite lindzen_1981 except at tropospheric downstream regions. !! \n Wave breaking occurs when \f$Ri_{m} - Calculate outputs: A, B, DUSFC, DVSFC (see parameter description). -!! - Below the dividing streamline height (k < idxzb), mountain -!! blocking(\f$D_{b}\f$) is applied. +!! - Below the dividing streamline height (k < idxzb), mountain +!! blocking(\f$D_{b}\f$) is applied. !! - Otherwise (k>= idxzb), orographic GWD (\f$\tau\f$) is applied. DO K = 1,KM DO I = 1,npt @@ -1354,7 +1354,7 @@ subroutine gwdps_run( & A(J,K) = - DBIM * V1(J,K) + A(J,K) B(J,K) = - DBIM * U1(J,K) + B(J,K) ENG1 = ENG0*(1.0-DBIM*DELTIM)*(1.0-DBIM*DELTIM) -! if ( ABS(DBIM * U1(J,K)) .gt. .01 ) +! if ( ABS(DBIM * U1(J,K)) .gt. .01 ) ! & print *,' in gwdps_lmi.f KDT=',KDT,I,K,DB(I,K), ! & dbim,idxzb(I),U1(J,K),V1(J,K),me DUSFC(J) = DUSFC(J) - DBIM * U1(J,K) * DEL(J,K) @@ -1384,7 +1384,7 @@ subroutine gwdps_run( & DUSFC(J) = TEM * DUSFC(J) DVSFC(J) = TEM * DVSFC(J) ENDDO -! +! ! MONITOR FOR EXCESSIVE GRAVITY WAVE DRAG TENDENCIES IF NCNT>0 ! ! IF(NCNT.GT.0) THEN From 1d277e3623440d69521c08455d35ff97e910b5c1 Mon Sep 17 00:00:00 2001 From: "Man.Zhang" Date: Mon, 25 Jun 2018 16:13:48 -0600 Subject: [PATCH 3/6] omitting white spaces in fortran codes. --- physics/gwdc.f | 106 +++++++++++++++++++++----------------------- physics/gwdps.f | 114 +++++++++++++++++++++++------------------------- 2 files changed, 104 insertions(+), 116 deletions(-) diff --git a/physics/gwdc.f b/physics/gwdc.f index 33d28b9c4..ab257e5f0 100644 --- a/physics/gwdc.f +++ b/physics/gwdc.f @@ -1,6 +1,6 @@ !> \file gwdc.f This file is the original code for parameterization of -!! stationary convection forced gravity wave drag based on -!! \cite chun_and_baik_1998 . +!! stationary convection forced gravity wave drag based on +!! \cite chun_and_baik_1998. module gwdc_pre contains @@ -101,8 +101,6 @@ end subroutine gwdc_pre_finalize end module gwdc_pre - - module gwdc contains @@ -196,7 +194,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & ! 2013 S. Moorthi - Updated and optimized code for T1534 GFS implementation ! ??? ?? 2015 J. Alpert - reducing the magnitude of tauctmax to fix blow up in L64 GFS ! S. Kar & M. Young -! aug 15 2016 - S. Moorthi - Fix for exessive dissipation which led to blow up in +! aug 15 2016 - S. Moorthi - Fix for exessive dissipation which led to blow up in ! 128 level runs with NEMS/GSM !*********************************************************************** @@ -216,7 +214,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & ! dpmid : midpoint delta p ( pi(k)-pi(k-1) ) ! lat : latitude index ! qmax : deep convective heating -! kcldtop : Vertical level index for cloud top ( mid level ) +! kcldtop : Vertical level index for cloud top ( mid level ) ! kcldbot : Vertical level index for cloud bottom ( mid level ) ! kcnv : (0,1) dependent on whether convection occur or not ! @@ -257,19 +255,19 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & ! parallel to the wind vector at the cloud top ( mid level ) ! tauctx : Wave stress at the cloud top projected in the east ! taucty : Wave stress at the cloud top projected in the north -! qmax : Maximum deep convective heating rate ( K s-1 ) in a +! qmax : Maximum deep convective heating rate ( K s-1 ) in a ! horizontal grid point calculated from cumulus para- ! meterization. ( mid level ) ! wtgwc : Wind tendency in direction to the wind vector at the cloud top level ! due to convectively generated gravity waves ( mid level ) -! utgwcl : Zonal wind tendency due to convectively generated +! utgwcl : Zonal wind tendency due to convectively generated ! gravity waves ( mid level ) ! vtgwcl : Meridional wind tendency due to convectively generated ! gravity waves ( mid level ) ! taugwci : Profile of wave stress calculated using basic-wind -! parallel to the wind vector at the cloud top +! parallel to the wind vector at the cloud top ! taugwcxi : Profile of zonal component of gravity wave stress -! taugwcyi : Profile of meridional component of gravity wave stress +! taugwcyi : Profile of meridional component of gravity wave stress ! ! taugwci, taugwcxi, and taugwcyi are defined at the interface level ! @@ -279,7 +277,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & ! rhom : Air density ( mid level ) ! ti : Temperature ( interface level ) ! basicum : Basic-wind profile. Basic-wind is parallel to the wind -! vector at the cloud top level. (mid level) +! vector at the cloud top level. (mid level) ! basicui : Basic-wind profile. Basic-wind is parallel to the wind ! vector at the cloud top level. ( interface level ) ! riloc : Local Richardson number ( interface level ) @@ -289,7 +287,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & ! break : Horizontal location where wave breaking is occurred. ! critic : Horizontal location where critical level filtering is ! occurred. -! dogwdc : Logical flag whether the GWDC parameterization is +! dogwdc : Logical flag whether the GWDC parameterization is ! calculated at a grid point or not. ! ! dogwdc is used in order to lessen CPU time for GWDC calculation. @@ -342,7 +340,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & ! ucltop : Zonal wind at the cloud top ( mid level ) ! vcltop : Meridional wind at the cloud top ( mid level ) ! windcltop : Wind speed at the cloud top ( mid level ) -! shear : Vertical shear of basic wind +! shear : Vertical shear of basic wind ! cosphi : Cosine of angle of wind vector at the cloud top ! sinphi : Sine of angle of wind vector at the cloud top ! c1 : Tunable parameter @@ -645,7 +643,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & qtem = spfh(i,k-1) * tem1 + spfh(i,k) * tem2 rhoi(i,k) = pint(i,k) / ( rd * ti(i,k)*(1.0+fv*qtem) ) dtdp = (t(i,k)-t(i,k-1)) / (pmid(i,k)-pmid(i,k-1)) - n2 = gsqr / ti(i,k) * ( 1./cp - rhoi(i,k)*dtdp ) + n2 = gsqr / ti(i,k) * ( 1./cp - rhoi(i,k)*dtdp ) bruni(i,k) = sqrt (max (n2min, n2)) enddo enddo @@ -741,7 +739,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & !----------------------------------------------------------------------- ! !> - Calculate the basic state wind profiles projected in the direction of the -!! cloud top wind at mid level and interface level. +!! cloud top wind at mid level and interface level. ! \n U : Basic-wind speed profile. Basic-wind is parallel to the wind ! vector at the cloud top level. (mid level) ! \n UI: Basic-wind speed profile. Basic-wind is parallel to the wind @@ -804,7 +802,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & ! 18 -------- U(18) ! 19 ======== UI(19) rhoi(19) bruni(19) riloc(19) ! -!----------------------------------------------------------------------- +!----------------------------------------------------------------------- do k=2,km do i=1,npt @@ -860,33 +858,33 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & ! it can be thought that there is deep convective cloud. However, ! deep convective heating between kcldbot and kcldtop is sometimes ! zero in spite of kcldtop less than kcldbot. In this case, -! maximum deep convective heating is assumed to be 1.e-30. +! maximum deep convective heating is assumed to be 1.e-30. ! ! B : kk is the vertical index for interface level cloud top ! ! C : Total convective fractional cover (cldf) is used as the -! convective cloud cover for GWDC calculation instead of +! convective cloud cover for GWDC calculation instead of ! convective cloud cover in each layer (concld). ! a1 = cldf*dlength ! You can see the difference between cldf(i) and concld(i) -! in (4.a.2) in Description of the NCAR Community Climate +! in (4.a.2) in Description of the NCAR Community Climate ! Model (CCM3). ! In NCAR CCM3, cloud fractional cover in each layer in a deep ! cumulus convection is determined assuming total convective -! cloud cover is randomly overlapped in each layer in the +! cloud cover is randomly overlapped in each layer in the ! cumulus convection. ! ! D : Wave stress at cloud top is calculated when the atmosphere ! is dynamically stable at the cloud top ! -! E : Cloud top wave stress and nonlinear parameter are calculated +! E : Cloud top wave stress and nonlinear parameter are calculated ! using density, temperature, and wind that are defined at mid ! level just below the interface level in which cloud top wave ! stress is defined. ! Nonlinct is defined at the interface level. ! ! F : If the atmosphere is dynamically unstable at the cloud top, -! GWDC calculation in current horizontal grid is skipped. +! GWDC calculation in current horizontal grid is skipped. ! ! G : If mean wind at the cloud top is less than zero, GWDC ! calculation in current horizontal grid is skipped. @@ -906,11 +904,11 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & !! \mu=\frac{gQ_{0}a_{1}}{c_{p}T_{0}NU^{2}} !! \f] !! where \f$Q_{0}\f$ is the maximum deep convective heating rate in a -!! horizontal grid point calculated from cumulus parameterization. +!! horizontal grid point calculated from cumulus parameterization. !! \f$a_{1}\f$ is the half-width of -!! the forcing function.\f$g\f$ is gravity. \f$c_{p}\f$ is specific -!! heat at constant pressure. \f$T_{0}\f$ is the layer mean -!! temperature (T1). As eqs.(18) and (19) \cite chun_and_baik_1998, +!! the forcing function.\f$g\f$ is gravity. \f$c_{p}\f$ is specific +!! heat at constant pressure. \f$T_{0}\f$ is the layer mean +!! temperature (T1). As eqs.(18) and (19) \cite chun_and_baik_1998, !! the zonal momentum flux is given by !! \f[ !! \tau_{x}=-[\rho U^{3}/(N\triangle x)]G(\mu) @@ -920,11 +918,11 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & !! G(\mu)=c_{1}c_2^2 \mu^{2} !! \f] !! wher \f$\rho\f$ is the local density. -!! The tunable parameter \f$c_1\f$ is related to the horizontal -!! structure of thermal forcing. The tunable parameter \f$c_2\f$ is -!! related to the basic-state wind and stability and the bottom and +!! The tunable parameter \f$c_1\f$ is related to the horizontal +!! structure of thermal forcing. The tunable parameter \f$c_2\f$ is +!! related to the basic-state wind and stability and the bottom and !! top heights of thermal forcing. If the atmosphere is dynamically -!! unstable at the cloud top, the convective GWD calculation is +!! unstable at the cloud top, the convective GWD calculation is !! skipped at that grid point. !! ! - If mean wind at the cloud top is less than zero, GWDC @@ -993,11 +991,11 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & ! ! Minimum RI is calculated for the following two cases ! -! (1) RIloc < 1.e+20 +! (1) RIloc < 1.e+20 ! (2) Riloc = 1.e+20 ----> Vertically uniform basic-state wind ! ! RIloc cannot be smaller than zero because N^2 becomes 1.E-32 in the -! case of N^2 < 0.. Thus the sign of RINUM is determined by +! case of N^2 < 0.. Thus the sign of RINUM is determined by ! 1 - nonlin*|c2|. ! !----------------------------------------------------------------------- @@ -1068,20 +1066,20 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & ! 2 ======== - 0.001 -1.e20 ! 2 -------- 0.000 ! 3 ======== - 0.001 -1.e20 -! 3 -------- -.xxx +! 3 -------- -.xxx ! 4 ======== - 0.001 2.600 2.000 ! 4 -------- 0.000 ! 5 ======== - 0.001 2.500 2.000 ! 5 -------- 0.000 ! 6 ======== - 0.001 1.500 0.110 -! 6 -------- +.xxx +! 6 -------- +.xxx ! 7 ======== - 0.005 2.000 3.000 ! 7 -------- 0.000 ! 8 ======== - 0.005 1.000 0.222 ! 8 -------- +.xxx ! 9 ======== - 0.010 1.000 2.000 ! 9 -------- 0.000 -! kcldtopi 10 ======== $$$ - 0.010 +! kcldtopi 10 ======== $$$ - 0.010 ! kcldtop 10 -------- $$$ yyyyy ! 11 ======== $$$ 0 ! 11 -------- $$$ @@ -1128,22 +1126,22 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & taugwci(i,k) = taugwci(i,k+1) elseif (rimin(i,k) > riminp) then tem = 2.0 + 1.0 / sqrt(riloc(i,k)) - nonlins = (1.0/abs(c2)) * (2.*sqrt(tem) - tem) + nonlins = (1.0/abs(c2)) * (2.*sqrt(tem) - tem) tem1 = basicui(i,k) tem2 = c2*nonlins*tem1 taugwci(i,k) = - rhoi(i,k) * c1 * tem1 * tem2 * tem2 & / (bruni(i,k)*dlen(i)) elseif (rimin(i,k) > riminm) then - taugwci(i,k) = zero + taugwci(i,k) = zero ! taugwci(i,k) = taugwci(i,k+1) end if ! RImin else -!> - If the minimum \f$R_{i}\f$ at interface cloud top is less than +!> - If the minimum \f$R_{i}\f$ at interface cloud top is less than !! or equal to 1/4, the convective GWD calculation is skipped at that !! grid point. - taugwci(i,k) = zero + taugwci(i,k) = zero end if ! RIloc else taugwci(i,k) = zero @@ -1189,7 +1187,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & taugw(i,k) = (taugwci(i,k+1) - taugwci(i,k)) / dpmid(i,k) if (taugw(i,k) /= 0.0) then tem = deltim * taugw(i,k) - dtfac(i) = min(dtfac(i), abs(velco(i,k)/tem)) + dtfac(i) = min(dtfac(i), abs(velco(i,k)/tem)) endif else taugw(i,k) = 0.0 @@ -1200,8 +1198,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & enddo enddo -!!!!!! Vertical differentiation -!!!!!! +!!!!!! Vertical differentiation!!!!!! !> - Calculate wind tendency in direction to the wind vector,zonal !! wind tendency and meridional wind tendency above the cloud top !! level due to convectively generated gravity waves. @@ -1294,22 +1291,19 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & !----------------------------------------------------------------------- ! -! The GWDC should accelerate the zonal and meridional wind in the -! opposite direction of the previous zonal and meridional wind, +! The GWDC should accelerate the zonal and meridional wind in the +! opposite direction of the previous zonal and meridional wind, ! respectively ! !----------------------------------------------------------------------- ! do k=1,kcldtop(i)-1 - ! if (utgwcl(i,k)*u(i,k) .gt. 0.0) then - !-------------------- x-component------------------- - -! write(6,'(a)') +! write(6,'(a)') ! + '(GWDC) WARNING: The GWDC should accelerate the zonal wind ' -! write(6,'(a,a,i3,a,i3)') -! + 'in the opposite direction of the previous zonal wind', +! write(6,'(a,a,i3,a,i3)') +! + 'in the opposite direction of the previous zonal wind', ! + ' at I = ',i,' and J = ',lat ! write(6,'(4(1x,e17.10))') u(i,kk),v(i,kk),u(i,k),v(i,k) ! write(6,'(a,1x,e17.10))') 'Vcld . V =', @@ -1319,7 +1313,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & ! do k1=1,km ! write(6,'(i2,36x,2(1x,e17.10))') ! + k1,taugwcxi(i,k1),taugwci(i,k1) -! write(6,'(i2,2(1x,e17.10))') k1,utgwcl(i,k1),u(i,k1) +! write(6,'(i2,2(1x,e17.10))') k1,utgwcl(i,k1),u(i,k1) ! end do ! write(6,'(i2,36x,1x,e17.10)') (km+1),taugwcxi(i,km+1) ! end if @@ -1329,7 +1323,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & ! do k1=1,km ! write(6,'(i2,36x,2(1x,e17.10))') ! + k1,taugwci(i,k1) -! write(6,'(i2,2(1x,e17.10))') k1,wtgwc(i,k1),basicum(i,k1) +! write(6,'(i2,2(1x,e17.10))') k1,wtgwc(i,k1),basicum(i,k1) ! end do ! write(6,'(i2,36x,1x,e17.10)') (km+1),taugwci(i,km+1) @@ -1349,7 +1343,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & ! do k1=1,km ! write(6,'(i2,36x,2(1x,e17.10))') ! + k1,taugwcyi(i,k1),taugwci(i,k1) -! write(6,'(i2,2(1x,e17.10))') k1,vtgwc(i,k1),v(i,k1) +! write(6,'(i2,2(1x,e17.10))') k1,vtgwc(i,k1),v(i,k1) ! end do ! write(6,'(i2,36x,1x,e17.10)') (km+1),taugwcyi(i,km+1) ! end if @@ -1379,7 +1373,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & !----------------------------------------------------------------------- ! -! For GWDC performance analysis +! For GWDC performance analysis ! !----------------------------------------------------------------------- @@ -1395,7 +1389,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & ! if ( abs(taugwci(i,k)-taugwci(i,kk)) > taumin ) then ! break(i) = 1.0 ! go to 2000 -! endif +! endif ! enddo !2000 continue @@ -1490,8 +1484,6 @@ end subroutine gwdc_finalize end module gwdc - - module gwdc_post contains diff --git a/physics/gwdps.f b/physics/gwdps.f index d5f9be09d..5f883ed70 100644 --- a/physics/gwdps.f +++ b/physics/gwdps.f @@ -361,7 +361,7 @@ subroutine gwdps_run( & !----- ALSO INCLUDED IS RI SMOOTH OVER A THICK LOWER LAYER !----- ALSO INCLUDED IS DECREASE IN DE-ACC AT TOP BY 1/2 !----- THE NMC GWD INCORPORATING BOTH GLAS(P&S) AND GFDL(MIGWD) -!----- MOUNTAIN INDUCED GRAVITY WAVE DRAG +!----- MOUNTAIN INDUCED GRAVITY WAVE DRAG !----- CODE FROM .FR30(V3MONNX) FOR MONIN3 !----- THIS VERSION (06 MAR 1987) !----- THIS VERSION (26 APR 1987) 3.G @@ -437,8 +437,7 @@ subroutine gwdps_run( & ! OTHER INPUT VARIABLES UNMODIFIED. ! revision log: ! May 2013 J. Wang change cleff back to opn setting -! Jan 2014 J. Wang merge Henry and Fangin's dissipation heat in gfs to nems -! +! Jan 2014 J. Wang merge Henry and Fangin's dissipation heat in gfs to nems ! ! ******************************************************************** USE MACHINE , ONLY : kind_phys @@ -497,8 +496,8 @@ subroutine gwdps_run( & parameter (FRC=1.0, CE=0.8, CEOFRC=CE/FRC, frmax=100., CG=0.5) parameter (GMAX=1.0, VELEPS=1.0, FACTOP=0.5) ! parameter (GMAX=1.0, CRITAC=5.0E-4, VELEPS=1.0, FACTOP=0.5) - parameter (RLOLEV=50000.0) -! parameter (RLOLEV=500.0) + parameter (RLOLEV=50000.0) +! parameter (RLOLEV=500.0) ! parameter (RLOLEV=0.5) ! real(kind=kind_phys) dpmin,hminmt,hncrit,minwnd,sigfac @@ -592,13 +591,13 @@ subroutine gwdps_run( & LCAPP1 = LCAP + 1 ! ! - IF ( NMTVR .eq. 14) then + IF ( NMTVR .eq. 14) then ! ---- for lm and gwd calculation points - RDXZB(:) = 0. ! CCPP: it is a bug + RDXZB(:) = 0. ! CCPP: it was a bug ipt = 0 npt = 0 DO I = 1,IM - IF ( (elvmax(i) .GT. HMINMT) + IF ( (elvmax(i) .GT. HMINMT) & .and. (hprime(i) .GT. hpmin) ) then npt = npt + 1 ipt(npt) = i @@ -616,7 +615,7 @@ subroutine gwdps_run( & ! do i=1,npt iwklm(i) = 2 - IDXZB(i) = 0 + IDXZB(i) = 0 kreflm(i) = 0 enddo ! if (lprnt) @@ -632,7 +631,7 @@ subroutine gwdps_run( & ! then do not need hncrit -- test with large hncrit first. ! KMLL = km / 2 ! maximum mtnlm height : # of vertical levels / 2 KMLL = kmm1 -! --- No mtn should be as high as KMLL (so we do not have to start at +! --- No mtn should be as high as KMLL (so we do not have to start at ! --- the top of the model but could do calc for all levels). ! DO I = 1, npt @@ -648,12 +647,12 @@ subroutine gwdps_run( & pkp1log = phil(j,k+1) / G pklog = phil(j,k) / G !!!------- ELVMAX(J) = min (ELVMAX(J) + sigfac * hprime(j), hncrit) - if ( ( ELVMAX(j) .le. pkp1log ) .and. + if ( ( ELVMAX(j) .le. pkp1log ) .and. & ( ELVMAX(j) .ge. pklog ) ) THEN ! print *,' in gwdps_lm.f 1 =',k,ELVMAX(j),pklog,pkp1log,me -! --- wk for diags but can be saved and reused. +! --- wk for diags but can be saved and reused. wk(i) = G * ELVMAX(j) / ( phil(j,k+1) - phil(j,k) ) - iwklm(I) = MAX(iwklm(I), k+1 ) + iwklm(I) = MAX(iwklm(I), k+1 ) ! print *,' in gwdps_lm.f 2 npt=',npt,i,j,wk(i),iwklm(i),me endif ! @@ -671,16 +670,15 @@ subroutine gwdps_run( & ! jhit = 0 ! do i = 1, npt ! j=ipt(i) -! if ( iwklm(i) .gt. ihit ) then +! if ( iwklm(i) .gt. ihit ) then ! ihit = iwklm(i) ! jhit = j ! endif ! enddo ! print *, ' mb: kdt,max(iwklm),jhit,phil,me=', ! & kdt,ihit,jhit,phil(jhit,ihit),me - klevm1 = KMLL - 1 - DO K = 1, klevm1 + DO K = 1, klevm1 DO I = 1, npt j = ipt(i) RDZ = g / ( phil(j,k+1) - phil(j,k) ) @@ -705,7 +703,7 @@ subroutine gwdps_run( & BNV2bar(I) = (PRSL(J,1)-PRSL(J,2)) * DELKS1(I) * BNV2LM(I,1) ENDDO -! --- find the dividing stream line height +! --- find the dividing stream line height ! --- starting from the level above the max mtn downward ! --- iwklm(i) is the k-index of mtn elvmax elevation !> - Find the dividing streamline height starting from the level above @@ -723,14 +721,13 @@ subroutine gwdps_run( & ! --- make averages, guess dividing stream (DS) line layer. ! --- This is not used in the first cut except for testing and ! --- is the vert ave of quantities from the surface to mtn top. -! DO I = 1, npt DO K = 1, Kreflm(I) J = ipt(i) RDELKS = DEL(J,K) * DELKS(I) - UBAR(I) = UBAR(I) + RDELKS * U1(J,K) ! trial Mean U below - VBAR(I) = VBAR(I) + RDELKS * V1(J,K) ! trial Mean V below - ROLL(I) = ROLL(I) + RDELKS * RO(I,K) ! trial Mean RO below + UBAR(I) = UBAR(I) + RDELKS * U1(J,K) ! trial Mean U below + VBAR(I) = VBAR(I) + RDELKS * V1(J,K) ! trial Mean V below + ROLL(I) = ROLL(I) + RDELKS * RO(I,K) ! trial Mean RO below RDELKS = (PRSL(J,K)-PRSL(J,K+1)) * DELKS1(I) BNV2bar(I) = BNV2bar(I) + BNV2lm(I,K) * RDELKS ! --- these vert ave are for diags, testing and GWD to follow (*j*). @@ -739,7 +736,7 @@ subroutine gwdps_run( & ! print *,' in gwdps_lm.f 5 =',i,kreflm(npt),BNV2bar(npt),me ! ! --- integrate to get PE in the trial layer. -! --- Need the first layer where PE>EK - as soon as +! --- Need the first layer where PE>EK - as soon as ! --- IDXZB is not 0 we have a hit and Zb is found. ! DO I = 1, npt @@ -755,21 +752,21 @@ subroutine gwdps_run( & !!\f[ !! UDS=\max(\sqrt{U1^2+V1^2},minwnd) !!\f] -!! where \f$ minwnd=0.1 \f$, \f$U1\f$ and \f$V1\f$ are zonal and +!! where \f$ minwnd=0.1 \f$, \f$U1\f$ and \f$V1\f$ are zonal and !! meridional wind components of model layer wind. - UDS(I,K) = + UDS(I,K) = & MAX(SQRT(U1(J,K)*U1(J,K) + V1(J,K)*V1(J,K)), minwnd) ! --- Test to see if we found Zb previously IF (IDXZB(I) .eq. 0 ) then - PE(I) = PE(I) + BNV2lm(I,K) * - & ( G * ELVMAX(J) - phil(J,K) ) * + PE(I) = PE(I) + BNV2lm(I,K) * + & ( G * ELVMAX(J) - phil(J,K) ) * & ( PHII(J,K+1) - PHII(J,K) ) / (G*G) ! --- KE ! --- Wind projected on the line perpendicular to mtn range, U(Zb(K)). ! --- kenetic energy is at the layer Zb ! --- THETA ranges from -+90deg |_ to the mtn "largest topo variations" UP(I) = UDS(I,K) * cos(ANG(I,K)) - EK(I) = 0.5 * UP(I) * UP(I) + EK(I) = 0.5 * UP(I) * UP(I) ! --- Dividing Stream lime is found when PE =exceeds EK. IF ( PE(I) .ge. EK(I) ) THEN @@ -777,14 +774,13 @@ subroutine gwdps_run( & RDXZB(J) = real(K,kind=kind_phys) ENDIF ! --- Then mtn blocked flow is between Zb=k(IDXZB(I)) and surface -! -!> - The dividing streamline height (idxzb), of a subgrid scale -!! obstable, is found by comparing the potential (PE) and kinetic +!> - The dividing streamline height (idxzb), of a subgrid scale +!! obstable, is found by comparing the potential (PE) and kinetic !! energies (EK) of the upstream large scale wind and subgrid scale air !! parcel movements. the dividing streamline is found when -!! \f$PE\geq EK\f$. Mountain-blocked flow is defined to exist between -!! the surface and the dividing streamline height (\f$h_d\f$), which -!! can be found by solving an integral equation for \f$h_d\f$: +!! \f$PE\geq EK\f$. Mountain-blocked flow is defined to exist between +!! the surface and the dividing streamline height (\f$h_d\f$), which +!! can be found by solving an integral equation for \f$h_d\f$: !!\f[ !! \frac{U^{2}(h_{d})}{2}=\int_{h_{d}}^{H} N^{2}(z)(H-z)dz !!\f] @@ -823,7 +819,7 @@ subroutine gwdps_run( & J = ipt(i) ZLEN = 0. ! print *,' in gwdps_lm.f 9 =',i,j,IDXZB(i),me - IF ( IDXZB(I) .gt. 0 ) then + IF ( IDXZB(I) .gt. 0 ) then DO K = IDXZB(I), 1, -1 IF ( PHIL(J,IDXZB(I)) .gt. PHIL(J,K) ) then @@ -832,31 +828,31 @@ subroutine gwdps_run( & !!\f[ !! ZLEN=\sqrt{[\frac{h_{d}-z}{z+h'}]} !!\f] -!! where \f$z\f$ is the height, \f$h'\f$ is the orographic standard +!! where \f$z\f$ is the height, \f$h'\f$ is the orographic standard !! deviation (HPRIME). - ZLEN = SQRT( ( PHIL(J,IDXZB(I)) - PHIL(J,K) ) / + ZLEN = SQRT( ( PHIL(J,IDXZB(I)) - PHIL(J,K) ) / & ( PHIL(J,K ) + G * hprime(J) ) ) ! --- lm eq 14: !> - Calculate the drag coefficient to vary with the aspect ratio of -!! the obstable as seen by the incident flow (see eq.14 in +!! the obstable as seen by the incident flow (see eq.14 in !! \cite lott_and_miller_1997) !!\f[ !! R=\frac{\cos^{2}\psi+\gamma\sin^{2}\psi}{\gamma\cos^{2}\psi+\sin^{2}\psi} !!\f] -!! where \f$\psi\f$, which is derived from THETA, is the angle between -!! the incident flow direction and the normal ridge direcion. +!! where \f$\psi\f$, which is derived from THETA, is the angle between +!! the incident flow direction and the normal ridge direcion. !! \f$\gamma\f$ is the orographic anisotropy (GAMMA). - R = (cos(ANG(I,K))**2 + GAMMA(J) * sin(ANG(I,K))**2) / + R = (cos(ANG(I,K))**2 + GAMMA(J) * sin(ANG(I,K))**2) / & (gamma(J) * cos(ANG(I,K))**2 + sin(ANG(I,K))**2) ! --- (negitive of DB -- see sign at tendency) -!> - In each model layer below the dividing streamlines, a drag from +!> - In each model layer below the dividing streamlines, a drag from !! the blocked flow is exerted by the obstacle on the large scale flow. !! The drag per unit area and per unit height is written (eq.15 in !! \cite lott_and_miller_1997): !!\f[ !! D_{b}(z)=-C_{d}\max(2-\frac{1}{R},0)\rho\frac{\sigma}{2h'}ZLEN\max(\cos\psi,\gamma\sin\psi)\frac{UDS}{2} !!\f] -!! where \f$C_{d}\f$ is a specified constant, \f$\sigma\f$ is the +!! where \f$C_{d}\f$ is a specified constant, \f$\sigma\f$ is the !! orographic slope. DBTMP = 0.25 * CDmb * @@ -881,8 +877,8 @@ subroutine gwdps_run( & !............................. ! end mtn blocking section ! - ELSEIF ( NMTVR .ne. 14) then -! ---- for mb not present and gwd (nmtvr .ne .14) + ELSEIF ( NMTVR .ne. 14) then +! ---- for mb not present and gwd (nmtvr .ne .14) ipt = 0 npt = 0 DO I = 1,IM @@ -982,13 +978,13 @@ subroutine gwdps_run( & enddo enddo ! -!> - Calculate the reference level index: kref=max(2,KPBL+1). where +!> - Calculate the reference level index: kref=max(2,KPBL+1). where !! KPBL is the index for the PBL top layer. KBPS = 1 KMPS = KM DO I=1,npt J = ipt(i) - kref(I) = MAX(IWK(I), KPBL(J)+1 ) ! reference level + kref(I) = MAX(IWK(I), KPBL(J)+1 ) ! reference level DELKS(I) = 1.0 / (PRSI(J,1) - PRSI(J,kref(I))) DELKS1(I) = 1.0 / (PRSL(J,1) - PRSL(J,kref(I))) UBAR (I) = 0.0 @@ -1024,7 +1020,7 @@ subroutine gwdps_run( & ! WD W S SW NW E N NE SE ! !> - Calculate low-level horizontal wind direction, the derived -!! orographic asymmetry parameter (OA), and the derived Lx (CLX). +!! orographic asymmetry parameter (OA), and the derived Lx (CLX). DO I = 1,npt J = ipt(i) wdir = atan2(UBAR(I),VBAR(I)) + pi @@ -1135,7 +1131,7 @@ subroutine gwdps_run( & !!\f[ !! a^{2}=C_{G}OC^{-1} !!\f] -!! where \f$F_{r_{c}}(=1)\f$ is the critical Froude number, +!! where \f$F_{r_{c}}(=1)\f$ is the critical Froude number, !! \f$F_{r_{0}}\f$ is the Froude number. \f$C_{E}\f$,\f$C_{m}\f$, !! \f$C_{G}\f$ are constants. @@ -1144,7 +1140,7 @@ subroutine gwdps_run( & !!\f[ !! \tau_0=E\frac{m'}{\triangle x}\frac{\rho_{0}U_0^3}{N_{0}}G' !!\f] -!! where \f$E\f$,\f$m'\f$, and \f$G'\f$ are the enhancement factor, +!! where \f$E\f$,\f$m'\f$, and \f$G'\f$ are the enhancement factor, !! "the number of mountains", and the flux function defined above, !! respectively. @@ -1228,7 +1224,7 @@ subroutine gwdps_run( & !! layer can be expressed using the drag for the layer immediately !! below. Thus, assuming \f$\tau_i=\tau_{i+1}\f$, we can get: !!\f[ -!! h_{d_i}^2=\frac{\triangle x}{m'}\frac{\tau_{i+1}}{\rho_{i}N_{i}U_{i}} +!! h_{d_i}^2=\frac{\triangle x}{m'}\frac{\tau_{i+1}}{\rho_{i}N_{i}U_{i}} !!\f] BRVF = SQRT(BNV2(I,K)) ! Brunt-Vaisala Frequency @@ -1240,9 +1236,9 @@ subroutine gwdps_run( & ! ! RIM is the MINIMUM-RICHARDSON NUMBER BY SHUTTS (1985) ! -!> - The minimum Richardson number (\f$Ri_{m}\f$) or local +!> - The minimum Richardson number (\f$Ri_{m}\f$) or local !! wave-modified Richardson number, which determines the onset of wave -!! breaking, is expressed in terms of \f$R_{i}\f$ and +!! breaking, is expressed in terms of \f$R_{i}\f$ and !! \f$F_{r_{d}}=Nh_{d}/U\f$: !!\f[ !! Ri_{m}=\frac{Ri(1-Fr_{d})}{(1+\sqrt{Ri}\cdot Fr_{d})^{2}} @@ -1259,17 +1255,17 @@ subroutine gwdps_run( & !> - Check stability to employ the 'saturation hypothesis' of !! \cite lindzen_1981 except at tropospheric downstream regions. !! \n Wave breaking occurs when \f$Ri_{m} Date: Tue, 26 Jun 2018 10:39:45 -0600 Subject: [PATCH 4/6] gwdps_post pass CCPP rt tests with option B --- physics/gwdps.f | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/physics/gwdps.f b/physics/gwdps.f index 5f883ed70..11f0ee871 100644 --- a/physics/gwdps.f +++ b/physics/gwdps.f @@ -785,7 +785,7 @@ subroutine gwdps_run( & !! \frac{U^{2}(h_{d})}{2}=\int_{h_{d}}^{H} N^{2}(z)(H-z)dz !!\f] !! where \f$H\f$ is the maximum subgrid scale elevation within the grid -!! box of actual orography, \f$h\f$, obtained from the GTOPO30 dataset +!! box of actual orography, \f$h\f$, obtained from the GTOPO30 dataset !! from the U.S. Geological Survey. ENDIF ENDDO @@ -858,10 +858,10 @@ subroutine gwdps_run( & DBTMP = 0.25 * CDmb * & MAX( 2. - 1. / R, 0. ) * sigma(J) * & MAX(cos(ANG(I,K)), gamma(J)*sin(ANG(I,K))) * - & ZLEN / hprime(J) - DB(I,K) = DBTMP * UDS(I,K) + & ZLEN / hprime(J) + DB(I,K) = DBTMP * UDS(I,K) ! -! if(lprnt .and. i .eq. npr) then +! if(lprnt .and. i .eq. npr) then ! print *,' in gwdps_lmi.f 10 npt=',npt,i,j,idxzb(i) ! &, DBTMP,R' ang=',ang(i,k),' gamma=',gamma(j),' K=',K ! print *,' in gwdps_lmi.f 11 K=',k,ZLEN,cos(ANG(I,K)) @@ -1019,7 +1019,7 @@ subroutine gwdps_run( & ! NWD 1 2 3 4 5 6 7 8 ! WD W S SW NW E N NE SE ! -!> - Calculate low-level horizontal wind direction, the derived +!> - Calculate low-level horizontal wind direction, the derived !! orographic asymmetry parameter (OA), and the derived Lx (CLX). DO I = 1,npt J = ipt(i) @@ -1108,7 +1108,7 @@ subroutine gwdps_run( & ! !> - Calculate enhancement factor (E),number of mountans (m') and !! aspect ratio constant. -!!\n As in eq.(4.9),(4.10),(4.11) in +!!\n As in eq.(4.9),(4.10),(4.11) in !! \cite kim_and_arakawa_1995, we define m' and E in such a way that they !! depend on the geometry and location of the subgrid-scale orography !! through OA and the nonlinearity of flow above the orography through @@ -1141,7 +1141,7 @@ subroutine gwdps_run( & !! \tau_0=E\frac{m'}{\triangle x}\frac{\rho_{0}U_0^3}{N_{0}}G' !!\f] !! where \f$E\f$,\f$m'\f$, and \f$G'\f$ are the enhancement factor, -!! "the number of mountains", and the flux function defined above, +!! "the number of mountains", and the flux function defined above, !! respectively. EFACT = (OA(I) + 2.) ** (CEOFRC*FR) @@ -1220,8 +1220,8 @@ subroutine gwdps_run( & !! \tau=\frac{m'}{\triangle x}\rho NUh_d^2 !!\f] !! where \f$h_{d}\f$ is the displacement wave amplitude. In the absence -!! of wave breaking, the displacement amplitude for the \f$i^{th}\f$ -!! layer can be expressed using the drag for the layer immediately +!! of wave breaking, the displacement amplitude for the \f$i^{th}\f$ +!! layer can be expressed using the drag for the layer immediately !! below. Thus, assuming \f$\tau_i=\tau_{i+1}\f$, we can get: !!\f[ !! h_{d_i}^2=\frac{\triangle x}{m'}\frac{\tau_{i+1}}{\rho_{i}N_{i}U_{i}} @@ -1252,7 +1252,7 @@ subroutine gwdps_run( & ! CHECK STABILITY TO EMPLOY THE 'SATURATION HYPOTHESIS' ! OF LINDZEN (1981) EXCEPT AT TROPOSPHERIC DOWNSTREAM REGIONS ! -!> - Check stability to employ the 'saturation hypothesis' of +!> - Check stability to employ the 'saturation hypothesis' of !! \cite lindzen_1981 except at tropospheric downstream regions. !! \n Wave breaking occurs when \f$Ri_{m} Date: Tue, 26 Jun 2018 16:55:36 -0600 Subject: [PATCH 5/6] gscond.f pass ccpp rt tests. --- physics/gscond.f | 70 +++++++++++++++++++++++------------------------- 1 file changed, 33 insertions(+), 37 deletions(-) diff --git a/physics/gscond.f b/physics/gscond.f index fc2525b0f..e5b39873b 100644 --- a/physics/gscond.f +++ b/physics/gscond.f @@ -1,6 +1,6 @@ !> \file gscond.f !! This file contains the subroutine that calculates grid-scale -!! condensation and evaporation for use in +!! condensation and evaporation for use in !! \cite zhao_and_carr_1997 scheme. module zhaocarr_gscond @@ -57,7 +57,7 @@ end subroutine zhaocarr_gscond_finalize !! cloud water and cloud ice (table2 of \cite zhao_and_carr_1997). !! -# Calculate the changes in \f$t\f$, \f$q\f$ and \f$p\f$ due to all the processes except microphysics. !! -# Calculate cloud evaporation rate (\f$E_c\f$, eq. 19 of \cite zhao_and_carr_1997) -!! -# Calculate cloud condensation rate (\f$C_g\f$, eq.8 of \cite zhao_and_carr_1997) +!! -# Calculate cloud condensation rate (\f$C_g\f$, eq.8 of \cite zhao_and_carr_1997) !! -# update t,q,cwm due to cloud evaporation and condensation process !> \section Zhao-Carr_cond_detailed GFS gscond Scheme Detailed Algorithm !> @{ @@ -121,7 +121,7 @@ subroutine zhaocarr_gscond_run (im,ix,km,dt,dtf,prsl,ps,q,clw1 & real (kind=kind_phys) qi(im), qint(im), ccrik, e0 &, cond, rdt, us, cclimit, climit &, tmt0, tmt15, qik, cwmik - &, ai, qw, u00ik, tik, pres, pp0, fi + &, ai, qw, u00ik, tik, pres, pp0, fi &, at, aq, ap, fiw, elv, qc, rqik &, rqikk, tx1, tx2, tx3, es, qs &, tsq, delq, condi, cone0, us00, ccrik1 @@ -173,7 +173,7 @@ subroutine zhaocarr_gscond_run (im,ix,km,dt,dtf,prsl,ps,q,clw1 & ! endif ! !************************************************************* -!> -# Begining of grid-scale condensation/evaporation loop (start of +!> -# Begining of grid-scale condensation/evaporation loop (start of !! k-loop, i-loop) !************************************************************* ! @@ -182,7 +182,7 @@ subroutine zhaocarr_gscond_run (im,ix,km,dt,dtf,prsl,ps,q,clw1 & ! vprs(:) = 0.001 * fpvs(t(:,k)) ! fpvs in pa !----------------------------------------------------------------------- !------------------qw, qi and qint-------------------------------------- - do i = 1, im + do i = 1, im tmt0 = t(i,k)-273.16 tmt15 = min(tmt0,cons_m15) qik = max(q(i,k),epsq) @@ -212,25 +212,25 @@ subroutine zhaocarr_gscond_run (im,ix,km,dt,dtf,prsl,ps,q,clw1 & !> -# Compute ice-water identification number IW. !!\n The distinction between cloud water and cloud ice is made by the !! cloud identification number IW, which is zero for cloud water and -!! unity for cloud ice (Table 2 in +!! unity for cloud ice (Table 2 in !! \cite zhao_and_carr_1997): -!! - All clouds are defined to consist of liquid water below the -!! freezing level (\f$T\geq 0^oC\f$) and of ice particles above the +!! - All clouds are defined to consist of liquid water below the +!! freezing level (\f$T\geq 0^oC\f$) and of ice particles above the !! \f$T=-15^oC\f$ level. -!! - In the temperature region between \f$-15^oC\f$ and \f$0^oC\f$, -!! clouds may be composed of liquid water or ice. If there are cloud +!! - In the temperature region between \f$-15^oC\f$ and \f$0^oC\f$, +!! clouds may be composed of liquid water or ice. If there are cloud !! ice particles above this point at the previous or current time step, -!! or if the cloud at this point at the previous time step consists of -!! ice particles, then the cloud substance at this point is considered +!! or if the cloud at this point at the previous time step consists of +!! ice particles, then the cloud substance at this point is considered !! to be ice particles because of the cloud seeding effect and the -!! memory of its content. Otherwise, all clouds in this region are +!! memory of its content. Otherwise, all clouds in this region are !! considered to contain supercooled cloud water. !-------------------ice-water id number iw------------------------------ if(tmt0.lt.-15.0) then u00ik = u(i,k) - fi = qik - u00ik*qi(i) - if(fi > d00.or.cwmik > climit) then + fi = qik - u00ik*qi(i) + if(fi > d00.or.cwmik > climit) then iw(i,k) = 1 else iw(i,k) = 0 @@ -251,8 +251,8 @@ subroutine zhaocarr_gscond_run (im,ix,km,dt,dtf,prsl,ps,q,clw1 & !> -# Condensation and evaporation of cloud !--------------condensation and evaporation of cloud-------------------- do i = 1, im -!> - Compute the changes in t, q and p (\f$A_{t}\f$,\f$A_{q}\f$ and -!! \f$A_{p}\f$) caused by all the processes except grid-scale +!> - Compute the changes in t, q and p (\f$A_{t}\f$,\f$A_{q}\f$ and +!! \f$A_{p}\f$) caused by all the processes except grid-scale !! condensation and evaporation. !!\f[ !! A_{t}=(t-tp)/dt @@ -274,7 +274,7 @@ subroutine zhaocarr_gscond_run (im,ix,km,dt,dtf,prsl,ps,q,clw1 & at = (tik-tp(i,k)) * rdt aq = (qik-qp(i,k)) * rdt ap = (pres-pp0) * rdt -!> - Calculate the saturation specific humidity \f$q_{s}\f$ and the +!> - Calculate the saturation specific humidity \f$q_{s}\f$ and the !! relative humidity \f$f\f$ using IW. !----------------the satuation specific humidity------------------------ fiw = float(iwik) @@ -283,7 +283,7 @@ subroutine zhaocarr_gscond_run (im,ix,km,dt,dtf,prsl,ps,q,clw1 & ! if (lprnt) print *,' qc=',qc,' qint=',qint(i),' qi=',qi(i) !----------------the relative humidity---------------------------------- if(qc.le.1.0e-10) then - rqik=d00 + rqik=d00 else rqik = qik/qc endif @@ -295,30 +295,30 @@ subroutine zhaocarr_gscond_run (im,ix,km,dt,dtf,prsl,ps,q,clw1 & !! b=1-\left ( \frac{f_{s}-f}{f_{s}-u} \right )^{1/2} !!\f] !! for \f$f>u\f$; and \f$b=0\f$ for \f$f -# End of the condensation/evaporation loop (end of i-loop,k-loop). !********************************************************************* ! - !> -# Store \f$t\f$, \f$q\f$, \f$ps\f$ for next time step. if (dt > dtf+0.001) then ! three time level @@ -519,8 +518,5 @@ subroutine zhaocarr_gscond_run (im,ix,km,dt,dtf,prsl,ps,q,clw1 & end subroutine zhaocarr_gscond_run !> @} ! @} - - ! @} - end module zhaocarr_gscond From d52a0beb0ad06ac54f100cdc2769ce05657369ad Mon Sep 17 00:00:00 2001 From: "Man.Zhang" Date: Wed, 27 Jun 2018 16:02:17 -0600 Subject: [PATCH 6/6] delete some white spaces. --- physics/gscond.f | 2 +- physics/gwdc.f | 28 ++++++++++++++-------------- physics/gwdps.f | 18 +++++++----------- 3 files changed, 22 insertions(+), 26 deletions(-) diff --git a/physics/gscond.f b/physics/gscond.f index e5b39873b..a9191a2bf 100644 --- a/physics/gscond.f +++ b/physics/gscond.f @@ -389,7 +389,7 @@ subroutine zhaocarr_gscond_run (im,ix,km,dt,dtf,prsl,ps,q,clw1 & ! if no cloud exists then evaporate any existing cloud condensate !----------------evaporation of cloud water----------------------------- e0 = d00 - if (ccrik <= cclimit.and. cwmik > climit) then + if (ccrik <= cclimit.and. cwmik > climit) then ! ! first iteration - increment halved ! diff --git a/physics/gwdc.f b/physics/gwdc.f index ab257e5f0..25036dfb9 100644 --- a/physics/gwdc.f +++ b/physics/gwdc.f @@ -289,7 +289,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & ! occurred. ! dogwdc : Logical flag whether the GWDC parameterization is ! calculated at a grid point or not. -! +! ! dogwdc is used in order to lessen CPU time for GWDC calculation. ! !----------------------------------------------------------------------- @@ -325,7 +325,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & & ti(:,:), riloc(:,:), & rimin(:,:), pint(:,:) ! real(kind=kind_phys), allocatable :: ugwdc(:,:), vgwdc(:,:), - real(kind=kind_phys), allocatable :: + real(kind=kind_phys), allocatable :: ! & plnmid(:,:), wtgwc(:,:), & plnmid(:,:), taugw(:,:), & utgwcl(:,:), vtgwcl(:,:), @@ -390,7 +390,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & vtgwc(i,k) = 0.0 ! brunm(i,k) = 0.0 ! rhom(i,k) = 0.0 - enddo + enddo enddo do i=1,im tauctx(i) = 0.0 @@ -466,7 +466,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & & rimin(npt,km+1), pint(npt,km+1)) ! allocate (ugwdc(npt,km), vgwdc(npt,km), - allocate + allocate ! & (plnmid(npt,km), wtgwc(npt,km), & (plnmid(npt,km), velco(npt,km), & utgwcl(npt,km), vtgwcl(npt,km), @@ -578,7 +578,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & ! 4 ======== pint(4) dpint(4) ! 4 -------- pmid(4) dpmid(4) ! ........ -! 17 ======== pint(17) dpint(17) +! 17 ======== pint(17) dpint(17) ! 17 -------- pmid(17) dpmid(17) ! 18 ======== pint(18) dpint(18) ! 18 -------- pmid(18) dpmid(18) @@ -658,7 +658,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & do k = 1, km do i = 1, npt dtdp = (ti(i,k+1)-ti(i,k)) / (pint(i,k+1)-pint(i,k)) - n2 = gsqr / t(i,k) * ( 1./cp - rhom(i,k)*dtdp ) + n2 = gsqr / t(i,k) * ( 1./cp - rhom(i,k)*dtdp ) brunm(i,k) = sqrt (max (n2min, n2)) enddo enddo @@ -814,10 +814,10 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & tem = bruni(i,k) / shear riloc(i,k) = tem * tem if (riloc(i,k) >= rimax ) riloc(i,k) = rilarge - end if + end if enddo enddo - + do i=1,npt riloc(i,1) = riloc(i,2) riloc(i,km+1) = riloc(i,km) @@ -882,7 +882,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & ! level just below the interface level in which cloud top wave ! stress is defined. ! Nonlinct is defined at the interface level. -! +! ! F : If the atmosphere is dynamically unstable at the cloud top, ! GWDC calculation in current horizontal grid is skipped. ! @@ -964,7 +964,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & tauctxl(i) = zero tauctyl(i) = zero do_gwc(i) = .false. - end if + end if !H enddo @@ -1094,16 +1094,16 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & ! 16 ======== $$$ 0 ! kcldbot 16 -------- $$$ ! 17 ======== 0 -! 17 -------- +! 17 -------- ! 18 ======== 0 -! 18 -------- +! 18 -------- ! 19 ======== 0 ! !----------------------------------------------------------------------- ! ! Even though the cloud top level obtained in deep convective para- ! meterization is defined in mid-level, the cloud top level for -! the GWDC calculation is assumed to be the interface level just +! the GWDC calculation is assumed to be the interface level just ! above the mid-level cloud top vertical level index. ! !----------------------------------------------------------------------- @@ -1133,7 +1133,7 @@ subroutine gwdc_run (im,ix,km,lat,u1,v1,t1,q1,deltim, & & / (bruni(i,k)*dlen(i)) elseif (rimin(i,k) > riminm) then taugwci(i,k) = zero -! taugwci(i,k) = taugwci(i,k+1) +! taugwci(i,k) = taugwci(i,k+1) end if ! RImin else diff --git a/physics/gwdps.f b/physics/gwdps.f index 11f0ee871..7a9c14ce5 100644 --- a/physics/gwdps.f +++ b/physics/gwdps.f @@ -814,7 +814,7 @@ subroutine gwdps_run( & ! endif ! ! --- The drag for mtn blocked flow -! +! DO I = 1, npt J = ipt(i) ZLEN = 0. @@ -853,7 +853,7 @@ subroutine gwdps_run( & !! D_{b}(z)=-C_{d}\max(2-\frac{1}{R},0)\rho\frac{\sigma}{2h'}ZLEN\max(\cos\psi,\gamma\sin\psi)\frac{UDS}{2} !!\f] !! where \f$C_{d}\f$ is a specified constant, \f$\sigma\f$ is the -!! orographic slope. +!! orographic slope. DBTMP = 0.25 * CDmb * & MAX( 2. - 1. / R, 0. ) * sigma(J) * @@ -872,7 +872,6 @@ subroutine gwdps_run( & ! if(lprnt) print *,' @K=1,ZLEN,DBTMP=',K,ZLEN,DBTMP endif ENDDO -! !............................. !............................. ! end mtn blocking section @@ -1049,7 +1048,6 @@ subroutine gwdps_run( & ULOW (I) = 0.0 DTFAC(I) = 1.0 ICRILV(I) = .FALSE. ! INITIALIZE CRITICAL LEVEL CONTROL VECTOR - ! !----COMPUTE THE "LOW LEVEL" WIND MAGNITUDE (M/S) ! @@ -1068,7 +1066,6 @@ subroutine gwdps_run( & ! ENDIF ENDDO ENDDO -! ! ! find the interface level of the projected wind where ! low levels & upper levels meet above pbl @@ -1106,7 +1103,7 @@ subroutine gwdps_run( & ! RATIO CONST. USE SIMPLIFIED RELATIONSHIP BETWEEN STANDARD ! DEVIATION & CRITICAL HGT ! -!> - Calculate enhancement factor (E),number of mountans (m') and +!> - Calculate enhancement factor (E),number of mountans (m') and !! aspect ratio constant. !!\n As in eq.(4.9),(4.10),(4.11) in !! \cite kim_and_arakawa_1995, we define m' and E in such a way that they @@ -1165,7 +1162,6 @@ subroutine gwdps_run( & SCOR(I) = BNV2(I,K) / TEM ! Scorer parameter below ref level ENDDO ! if(lprnt) print *,' taub=',taub -! !----SET UP BOTTOM VALUES OF STRESS ! DO K = 1, KBPS @@ -1211,7 +1207,7 @@ subroutine gwdps_run( & SCORK = BNV2(I,K) * TEMV * TEMV RSCOR = MIN(1.0, SCORK / SCOR(I)) SCOR(I) = SCORK - ELSE + ELSE RSCOR = 1. ENDIF ! @@ -1334,8 +1330,8 @@ subroutine gwdps_run( & ! endif !> - Calculate outputs: A, B, DUSFC, DVSFC (see parameter description). -!! - Below the dividing streamline height (k < idxzb), mountain -!! blocking(\f$D_{b}\f$) is applied. +!! - Below the dividing streamline height (k < idxzb), mountain +!! blocking(\f$D_{b}\f$) is applied. !! - Otherwise (k>= idxzb), orographic GWD (\f$\tau\f$) is applied. DO K = 1,KM DO I = 1,npt @@ -1350,7 +1346,7 @@ subroutine gwdps_run( & A(J,K) = - DBIM * V1(J,K) + A(J,K) B(J,K) = - DBIM * U1(J,K) + B(J,K) ENG1 = ENG0*(1.0-DBIM*DELTIM)*(1.0-DBIM*DELTIM) -! if ( ABS(DBIM * U1(J,K)) .gt. .01 ) +! if ( ABS(DBIM * U1(J,K)) .gt. .01 ) ! & print *,' in gwdps_lmi.f KDT=',KDT,I,K,DB(I,K), ! & dbim,idxzb(I),U1(J,K),V1(J,K),me DUSFC(J) = DUSFC(J) - DBIM * U1(J,K) * DEL(J,K)