| cplchm | gfs_control_type | logical flag for chemistry collection | .false.
diff --git a/physics/drag_suite.F90 b/physics/drag_suite.F90
new file mode 100644
index 000000000..eb371adb1
--- /dev/null
+++ b/physics/drag_suite.F90
@@ -0,0 +1,1473 @@
+!> \File drag_suite.F90
+!! This file is the parameterization of orographic gravity wave
+!! drag, mountain blocking, and form drag.
+
+!> This module contains the CCPP-compliant orographic gravity wave
+!! drag pre interstitial codes.
+ module drag_suite_pre
+
+ contains
+
+!> \section arg_table_drag_suite_pre_init Argument Table
+!!
+ subroutine drag_suite_pre_init()
+ end subroutine drag_suite_pre_init
+
+!> \section arg_table_drag_suite_pre_run Argument Table
+!! \htmlinclude drag_suite_pre_run.html
+!!
+!! \section general General Algorithm
+!! \section detailed Detailed Algorithm
+!! @{
+ subroutine drag_suite_pre_run( &
+ & im, nmtvr, mntvar, &
+ & hprime, oc, oa4, clx, theta, &
+ & sigma, gamma, elvmax, errmsg, errflg)
+
+ use machine, only : kind_phys
+ implicit none
+
+ integer, intent(in) :: im, nmtvr
+ real(kind=kind_phys), intent(in) :: mntvar(im,nmtvr)
+
+ real(kind=kind_phys), intent(out) :: &
+ & hprime(im), oc(im), oa4(im,4), clx(im,4), &
+ & theta(im), sigma(im), gamma(im), elvmax(im)
+
+ character(len=*), intent(out) :: errmsg
+ integer, intent(out) :: errflg
+
+ ! Initialize CCPP error handling variables
+ errmsg = ''
+ errflg = 0
+
+ if (nmtvr == 14) then ! current operational - as of 2014
+ hprime(:) = mntvar(:,1)
+ oc(:) = mntvar(:,2)
+ oa4(:,1) = mntvar(:,3)
+ oa4(:,2) = mntvar(:,4)
+ oa4(:,3) = mntvar(:,5)
+ oa4(:,4) = mntvar(:,6)
+ clx(:,1) = mntvar(:,7)
+ clx(:,2) = mntvar(:,8)
+ clx(:,3) = mntvar(:,9)
+ clx(:,4) = mntvar(:,10)
+ theta(:) = mntvar(:,11)
+ gamma(:) = mntvar(:,12)
+ sigma(:) = mntvar(:,13)
+ elvmax(:) = mntvar(:,14)
+ elseif (nmtvr == 10) then
+ hprime(:) = mntvar(:,1)
+ oc(:) = mntvar(:,2)
+ oa4(:,1) = mntvar(:,3)
+ oa4(:,2) = mntvar(:,4)
+ oa4(:,3) = mntvar(:,5)
+ oa4(:,4) = mntvar(:,6)
+ clx(:,1) = mntvar(:,7)
+ clx(:,2) = mntvar(:,8)
+ clx(:,3) = mntvar(:,9)
+ clx(:,4) = mntvar(:,10)
+ elseif (nmtvr == 6) then
+ hprime(:) = mntvar(:,1)
+ oc(:) = mntvar(:,2)
+ oa4(:,1) = mntvar(:,3)
+ oa4(:,2) = mntvar(:,4)
+ oa4(:,3) = mntvar(:,5)
+ oa4(:,4) = mntvar(:,6)
+ clx(:,1) = 0.0
+ clx(:,2) = 0.0
+ clx(:,3) = 0.0
+ clx(:,4) = 0.0
+ else
+ hprime = 0
+ oc = 0
+ oa4 = 0
+ clx = 0
+ theta = 0
+ gamma = 0
+ sigma = 0
+ elvmax = 0
+ endif ! end if_nmtvr
+
+ end subroutine drag_suite_pre_run
+!> @}
+
+! \ingroup GFS_ogwd
+! \brief Brief description of the subroutine
+!
+!> \section arg_table_drag_suite_pre_finalize Argument Table
+!!
+ subroutine drag_suite_pre_finalize()
+ end subroutine drag_suite_pre_finalize
+
+ end module drag_suite_pre
+
+!> This module contains the CCPP-compliant orographic gravity wave dray scheme.
+ module drag_suite
+
+ contains
+
+!> \section arg_table_drag_suite_init Argument Table
+!!
+ subroutine drag_suite_init()
+ end subroutine drag_suite_init
+
+! \defgroup GFS_ogwd GFS Orographic Gravity Wave Drag
+!> \defgroup gfs_drag_suite GFS drag_suite Main
+!! \brief This subroutine includes orographic gravity wave drag, mountain
+!! blocking, and form drag.
+!!
+!> 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.
+!!
+!> \section arg_table_drag_suite_run Argument Table
+!! \htmlinclude drag_suite_run.html
+!!
+!> \section gen_drag_suite 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_drag_suite GFS Orographic GWD Scheme Detailed Algorithm
+!> @{
+! subroutine drag_suite_run( &
+! & IM,IX,KM,A,B,C,U1,V1,T1,Q1,KPBL, &
+! & 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, rdxzb, errmsg, errflg)
+!
+ subroutine drag_suite_run( &
+ & IM,IX,KM,dvdt,dudt,dtdt,U1,V1,T1,Q1,KPBL, &
+ & PRSI,DEL,PRSL,PRSLK,PHII,PHIL,DELTIM,KDT, &
+ & VAR,oc1,oa4,ol4, &
+! & varss,oc1ss,oa4ss,ol4ss, &
+ & THETA,SIGMA,GAMMA,ELVMAX, &
+ & dtaux2d_ls,dtauy2d_ls,dtaux2d_bl,dtauy2d_bl, &
+ & dtaux2d_ss,dtauy2d_ss,dtaux2d_fd,dtauy2d_fd, &
+ & dusfc,dvsfc, &
+ & dusfc_ls,dvsfc_ls,dusfc_bl,dvsfc_bl, &
+ & dusfc_ss,dvsfc_ss,dusfc_fd,dvsfc_fd, &
+ & slmsk,br1,hpbl, &
+ & g, cp, rd, rv, fv, pi, imx, cdmbgwd, me, master, &
+ & lprnt, ipr, rdxzb, dx, gwd_opt, errmsg, errflg )
+
+! ********************************************************************
+! -----> I M P L E M E N T A T I O N V E R S I O N <----------
+!
+! ----- This code -----
+!begin WRF code
+
+! this code handles the time tendencies of u v due to the effect of mountain
+! induced gravity wave drag from sub-grid scale orography. this routine
+! not only treats the traditional upper-level wave breaking due to mountain
+! variance (alpert 1988), but also the enhanced lower-tropospheric wave
+! breaking due to mountain convexity and asymmetry (kim and arakawa 1995).
+! thus, in addition to the terrain height data in a model grid box,
+! additional 10-2d topographic statistics files are needed, including
+! orographic standard deviation (var), convexity (oc1), asymmetry (oa4)
+! and ol (ol4). these data sets are prepared based on the 30 sec usgs orography
+! hong (1999). the current scheme was implmented as in hong et al.(2008)
+!
+! Originally coded by song-you hong and young-joon kim and implemented by song-you hong
+!
+! program history log:
+! 2014-10-01 Hyun-Joo Choi (from KIAPS) flow-blocking drag of kim and doyle
+! with blocked height by dividing streamline theory
+! 2017-04-06 Joseph Olson (from Gert-Jan Steeneveld) added small-scale
+! orographic grabity wave drag:
+! 2017-09-15 Joseph Olson, with some bug fixes from Michael Toy: added the
+! topographic form drag of Beljaars et al. (2004, QJRMS)
+! Activation of each component is done by specifying the integer-parameters
+! (defined below) to 0: inactive or 1: active
+! gwd_opt_ls = 0 or 1: large-scale
+! gwd_opt_bl = 0 or 1: blocking drag
+! gwd_opt_ss = 0 or 1: small-scale gravity wave drag
+! gwd_opt_fd = 0 or 1: topographic form drag
+! 2017-09-25 Michael Toy (from NCEP GFS model) added dissipation heating
+! gsd_diss_ht_opt = 0: dissipation heating off
+! gsd_diss_ht_opt = 1: dissipation heating on
+!
+! References:
+! Hong et al. (2008), wea. and forecasting
+! Kim and Doyle (2005), Q. J. R. Meteor. Soc.
+! Kim and Arakawa (1995), j. atmos. sci.
+! Alpert et al. (1988), NWP conference.
+! Hong (1999), NCEP office note 424.
+! Steeneveld et al (2008), JAMC
+! Tsiringakis et al. (2017), Q. J. R. Meteor. Soc.
+! Beljaars et al. (2004), Q. J. R. Meteor. Soc.
+!
+! notice : comparible or lower resolution orography files than model resolution
+! are desirable in preprocess (wps) to prevent weakening of the drag
+!-------------------------------------------------------------------------------
+!
+! input
+! dudt (im,km) non-lin tendency for u wind component
+! dvdt (im,km) non-lin tendency for v wind component
+! u1(im,km) zonal wind / sqrt(rcl) m/sec at t0-dt
+! v1(im,km) meridional wind / sqrt(rcl) m/sec at t0-dt
+! t1(im,km) temperature deg k at t0-dt
+! q1(im,km) specific humidity at t0-dt
+! deltim time step secs
+! del(km) positive increment of pressure across layer (pa)
+! KPBL(IM) is the index of the top layer of the PBL
+! ipr & lprnt for diagnostics
+!
+! output
+! dudt, dvdt wind tendency due to gwdo
+! dTdt
+!
+!-------------------------------------------------------------------------------
+
+!end wrf code
+!----------------------------------------------------------------------C
+! USE
+! ROUTINE IS CALLED FROM CCPP (AFTER CALLING PBL SCHEMES)
+!
+! PURPOSE
+! USING THE GWD PARAMETERIZATIONS OF PS-GLAS AND PH-
+! GFDL TECHNIQUE. THE TIME TENDENCIES OF U V
+! ARE ALTERED TO INCLUDE THE EFFECT OF MOUNTAIN INDUCED
+! GRAVITY WAVE DRAG FROM SUB-GRID SCALE OROGRAPHY INCLUDING
+! CONVECTIVE BREAKING, SHEAR BREAKING AND THE PRESENCE OF
+! CRITICAL LEVELS
+!
+!
+! ********************************************************************
+ USE MACHINE , ONLY : kind_phys
+ implicit none
+
+ ! Interface variables
+ integer, intent(in) :: im, ix, km, imx, kdt, ipr, me, master
+ integer, intent(in) :: gwd_opt
+ logical, intent(in) :: lprnt
+ integer, intent(in) :: KPBL(im)
+ real(kind=kind_phys), intent(in) :: deltim, G, CP, RD, RV, cdmbgwd(2)
+
+ integer :: kpblmax
+ integer, parameter :: ims=1, kms=1, its=1, kts=1
+ real(kind=kind_phys), intent(in) :: fv, pi
+ real(kind=kind_phys) :: rcl, cdmb
+ real(kind=kind_phys) :: g_inv
+
+ real(kind=kind_phys), intent(out) :: &
+ & dudt(im,km),dvdt(im,km), &
+ & dtdt(im,km), rdxzb(im)
+ real(kind=kind_phys), intent(in) :: &
+ & u1(im,km),v1(im,km), &
+ & t1(im,km),q1(im,km), &
+ & PHII(im,km+1),prsl(im,km), &
+ & prslk(im,km),PHIL(im,km)
+ real(kind=kind_phys), intent(in) :: prsi(im,km+1), &
+ & del(im,km)
+ real(kind=kind_phys), intent(in) :: var(im),oc1(im), &
+ & oa4(im,4),ol4(im,4), &
+ & dx(im)
+ !real(kind=kind_phys), intent(in) :: varss(im),oc1ss(im), &
+ real(kind=kind_phys) :: varss(im),oc1ss(im), &
+ & oa4ss(im,4),ol4ss(im,4)
+ real(kind=kind_phys), intent(in) :: THETA(im),SIGMA(im), &
+ & GAMMA(im),ELVMAX(im)
+
+! added for small-scale orographic wave drag
+ real(kind=kind_phys), dimension(im,km) :: utendwave,vtendwave,thx,thvx
+ real(kind=kind_phys), intent(in) :: br1(im), &
+ & hpbl(im), &
+ & slmsk(im)
+ real(kind=kind_phys), dimension(im) :: govrth,xland
+ real(kind=kind_phys), dimension(im,km) :: dz2
+ real(kind=kind_phys) :: tauwavex0,tauwavey0, &
+ & XNBV,density,tvcon,hpbl2
+ integer :: kpbl2,kvar
+ real(kind=kind_phys), dimension(im,km+1) :: zq ! = PHII/g
+ real(kind=kind_phys), dimension(im,km) :: zl ! = PHIL/g
+
+!SPP
+ real(kind=kind_phys), dimension(im) :: rstoch
+
+!Output:
+ real(kind=kind_phys), intent(out) :: &
+ & dusfc(im), dvsfc(im)
+!Output (optional):
+ real(kind=kind_phys), intent(out) :: &
+ & dusfc_ls(:),dvsfc_ls(:), &
+ & dusfc_bl(:),dvsfc_bl(:), &
+ & dusfc_ss(:),dvsfc_ss(:), &
+ & dusfc_fd(:),dvsfc_fd(:)
+ real(kind=kind_phys), intent(out) :: &
+ & dtaux2d_ls(:,:),dtauy2d_ls(:,:), &
+ & dtaux2d_bl(:,:),dtauy2d_bl(:,:), &
+ & dtaux2d_ss(:,:),dtauy2d_ss(:,:), &
+ & dtaux2d_fd(:,:),dtauy2d_fd(:,:)
+
+!Misc arrays
+ real(kind=kind_phys), dimension(im,km) :: dtaux2d, dtauy2d
+
+!-------------------------------------------------------------------------
+! Flags to regulate the activation of specific components of drag suite:
+! Each component is tapered off automatically as a function of dx, so best to
+! keep them activated (=1).
+ integer, parameter :: &
+ gwd_opt_ls = 1, & ! large-scale gravity wave drag
+ gwd_opt_bl = 1, & ! blocking drag
+ gwd_opt_ss = 1, & ! small-scale gravity wave drag (Steeneveld et al. 2008)
+ gwd_opt_fd = 1, & ! form drag (Beljaars et al. 2004, QJRMS)
+ gsd_diss_ht_opt = 0
+
+! Parameters for bounding the scale-adaptive variability:
+! Small-scale GWD + turbulent form drag
+ real(kind=kind_phys), parameter :: dxmin_ss = 1000., &
+ & dxmax_ss = 12000. ! min,max range of tapering (m)
+! Large-scale GWD + blocking
+ real(kind=kind_phys), parameter :: dxmin_ls = 3000., &
+ & dxmax_ls = 13000. ! min,max range of tapering (m)
+ real(kind=kind_phys) :: ss_taper, ls_taper ! small- and large-scale tapering factors (-)
+!
+! Variables for limiting topographic standard deviation (var)
+ real(kind=kind_phys), parameter :: varmax_ss = 50., &
+ varmax_fd = 150., &
+ beta_ss = 0.1, &
+ beta_fd = 0.2
+ real(kind=kind_phys) :: var_temp, var_temp2
+
+! added Beljaars orographic form drag
+ real(kind=kind_phys), dimension(im,km) :: utendform,vtendform
+ real(kind=kind_phys) :: a1,a2,wsp
+ real(kind=kind_phys) :: H_efold
+
+! critical richardson number for wave breaking : ! larger drag with larger value
+ real(kind=kind_phys), parameter :: ric = 0.25
+ real(kind=kind_phys), parameter :: dw2min = 1.
+ real(kind=kind_phys), parameter :: rimin = -100.
+ real(kind=kind_phys), parameter :: bnv2min = 1.0e-5
+ real(kind=kind_phys), parameter :: efmin = 0.0
+ real(kind=kind_phys), parameter :: efmax = 10.0
+ real(kind=kind_phys), parameter :: xl = 4.0e4
+ real(kind=kind_phys), parameter :: critac = 1.0e-5
+ real(kind=kind_phys), parameter :: gmax = 1.
+ real(kind=kind_phys), parameter :: veleps = 1.0
+ real(kind=kind_phys), parameter :: factop = 0.5
+ real(kind=kind_phys), parameter :: frc = 1.0
+ real(kind=kind_phys), parameter :: ce = 0.8
+ real(kind=kind_phys), parameter :: cg = 0.5
+ integer,parameter :: kpblmin = 2
+
+!
+! local variables
+!
+ integer :: i,j,k,lcap,lcapp1,nwd,idir, &
+ klcap,kp1,ikount,kk
+!
+ real(kind=kind_phys) :: rcs,rclcs,csg,fdir,cleff,cleff_ss,cs, &
+ rcsks,wdir,ti,rdz,temp,tem2,dw2,shr2, &
+ bvf2,rdelks,wtkbj,tem,gfobnv,hd,fro, &
+ rim,temc,tem1,efact,temv,dtaux,dtauy, &
+ dtauxb,dtauyb,eng0,eng1
+!
+ logical :: ldrag(im),icrilv(im), &
+ flag(im),kloop1(im)
+!
+ real(kind=kind_phys) :: taub(im),taup(im,km+1), &
+ xn(im),yn(im), &
+ ubar(im),vbar(im), &
+ fr(im),ulow(im), &
+ rulow(im),bnv(im), &
+ oa(im),ol(im), &
+ oass(im),olss(im), &
+ roll(im),dtfac(im), &
+ brvf(im),xlinv(im), &
+ delks(im),delks1(im), &
+ bnv2(im,km),usqj(im,km), &
+ taud_ls(im,km),taud_bl(im,km), &
+ ro(im,km), &
+ vtk(im,km),vtj(im,km), &
+ zlowtop(im),velco(im,km-1), &
+ coefm(im),coefm_ss(im)
+!
+ integer :: kbl(im),klowtop(im)
+ logical :: iope
+ integer,parameter :: mdir=8
+ !integer :: nwdir(mdir)
+ !data nwdir/6,7,5,8,2,3,1,4/
+ integer, parameter :: nwdir(8) = (/6,7,5,8,2,3,1,4/)
+!
+! variables for flow-blocking drag
+!
+ real(kind=kind_phys),parameter :: frmax = 10.
+ real(kind=kind_phys),parameter :: olmin = 1.0e-5
+ real(kind=kind_phys),parameter :: odmin = 0.1
+ real(kind=kind_phys),parameter :: odmax = 10.
+ real(kind=kind_phys),parameter :: erad = 6371.315e+3
+ integer :: komax(im)
+ integer :: kblk
+ real(kind=kind_phys) :: cd
+ real(kind=kind_phys) :: zblk,tautem
+ real(kind=kind_phys) :: pe,ke
+ real(kind=kind_phys) :: delx,dely,dxy4(4),dxy4p(4)
+ real(kind=kind_phys) :: dxy(im),dxyp(im)
+ real(kind=kind_phys) :: ol4p(4),olp(im),od(im)
+ real(kind=kind_phys) :: taufb(im,km+1)
+
+ character(len=*), intent(out) :: errmsg
+ integer, intent(out) :: errflg
+
+ ! Calculate inverse of gravitational acceleration
+ g_inv = 1./G
+
+ ! Initialize CCPP error handling variables
+ errmsg = ''
+ errflg = 0
+
+if (me==master) print *,"Running drag suite"
+!--------------------------------------------------------------------
+! SCALE-ADPTIVE PARAMETER FROM GFS GWD SCHEME
+!--------------------------------------------------------------------
+! parameter (cdmb = 1.0) ! non-dim sub grid mtn drag Amp (*j*)
+! non-dim sub grid mtn drag Amp (*j*)
+! cdmb = 1.0/float(IMX/192)
+! cdmb = 192.0/float(IMX)
+ cdmb = 4.0 * 192.0/float(IMX)
+ if (cdmbgwd(1) >= 0.0) cdmb = cdmb * cdmbgwd(1)
+
+!>-# Orographic Gravity Wave Drag Section
+ kpblmax = km / 2 ! maximum pbl height : # of vertical levels / 2
+!
+! Scale cleff between IM=384*2 and 192*2 for T126/T170 and T62
+!
+ if (imx > 0) then
+! cleff = 1.0E-5 * SQRT(FLOAT(IMX)/384.0)
+! cleff = 1.0E-5 * SQRT(FLOAT(IMX)/192.0) ! this is inverse of CLEFF!
+! cleff = 0.5E-5 * SQRT(FLOAT(IMX)/192.0) ! this is inverse of CLEFF!
+! cleff = 1.0E-5 * SQRT(FLOAT(IMX)/192)/float(IMX/192)
+! cleff = 1.0E-5 / SQRT(FLOAT(IMX)/192.0) ! this is inverse of CLEFF!
+ cleff = 0.5E-5 / SQRT(FLOAT(IMX)/192.0) ! this is inverse of CLEFF!
+! hmhj for ndsl
+! jw cleff = 0.1E-5 / SQRT(FLOAT(IMX)/192.0) ! this is inverse of CLEFF!
+! cleff = 2.0E-5 * SQRT(FLOAT(IMX)/192.0) ! this is inverse of CLEFF!
+! cleff = 2.5E-5 * SQRT(FLOAT(IMX)/192.0) ! this is inverse of CLEFF!
+ endif
+ if (cdmbgwd(2) >= 0.0) cleff = cleff * cdmbgwd(2)
+!--------------------------------------------------------------------
+! END SCALE-ADPTIVE PARAMETER SECTION
+!--------------------------------------------------------------------
+!
+!---- constants
+!
+ rcl = 1.
+ rcs = sqrt(rcl)
+ cs = 1. / sqrt(rcl)
+ csg = cs * g
+ lcap = km
+ lcapp1 = lcap + 1
+ fdir = mdir / (2.0*pi)
+
+ do i=1,im
+ if (slmsk(i)==1. .or. slmsk(i)==2.) then !sea/land/ice mask (=0/1/2) in FV3
+ xland(i)=1.0 !but land/water = (1/2) in this module
+ else
+ xland(i)=2.0
+ endif
+ RDXZB(i) = 0.0
+ enddo
+
+!temporary use of large-scale data:
+ do i=1,im
+ varss(i)=var(i)
+ oc1ss(i)=oc1(i)
+ do j=1,4
+ oa4ss(i,j)=oa4(i,j)
+ ol4ss(i,j)=ol4(i,j)
+ enddo
+ enddo
+!
+!--- calculate scale-aware tapering factors
+!NOTE: if dx(1) is not representative of most/all dx, this needs to change...
+if ( dx(1) .ge. dxmax_ls ) then
+ ls_taper = 1.
+else
+ if ( dx(1) .le. dxmin_ls) then
+ ls_taper = 0.
+ else
+ ls_taper = 0.5 * ( SIN(pi*(dx(1)-0.5*(dxmax_ls+dxmin_ls))/ &
+ (dxmax_ls-dxmin_ls)) + 1. )
+ end if
+end if
+if (me==master) print *,"in Drag Suite, dx(1:2):",dx(1),dx(2)
+if ( dx(1) .ge. dxmax_ss ) then
+ ss_taper = 1.
+else
+ if ( dx(1) .le. dxmin_ss) then
+ ss_taper = 0.
+ else
+ ss_taper = dxmax_ss * (1. - dxmin_ss/dx(1))/(dxmax_ss-dxmin_ss)
+ end if
+end if
+if (me==master) print *,"in Drag Suite, ss_taper:",ss_taper
+
+!--- calculate length of grid for flow-blocking drag
+!
+ delx = dx(1)
+ dely = dx(1)
+ dxy4(1) = delx
+ dxy4(2) = dely
+ dxy4(3) = sqrt(delx*delx + dely*dely)
+ dxy4(4) = dxy4(3)
+ dxy4p(1) = dxy4(2)
+ dxy4p(2) = dxy4(1)
+ dxy4p(3) = dxy4(4)
+ dxy4p(4) = dxy4(3)
+!
+!-----initialize arrays
+!
+ dtaux = 0.0
+ dtauy = 0.0
+ do i = its,im
+ klowtop(i) = 0
+ kbl(i) = 0
+ enddo
+!
+ do i = its,im
+ xn(i) = 0.0
+ yn(i) = 0.0
+ ubar (i) = 0.0
+ vbar (i) = 0.0
+ roll (i) = 0.0
+ taub (i) = 0.0
+ oa(i) = 0.0
+ ol(i) = 0.0
+ oass(i) = 0.0
+ olss(i) = 0.0
+ ulow (i) = 0.0
+ dtfac(i) = 1.0
+ ldrag(i) = .false.
+ icrilv(i) = .false.
+ flag(i) = .true.
+ enddo
+
+ do k = kts,km
+ do i = its,im
+ usqj(i,k) = 0.0
+ bnv2(i,k) = 0.0
+ vtj(i,k) = 0.0
+ vtk(i,k) = 0.0
+ taup(i,k) = 0.0
+ taud_ls(i,k) = 0.0
+ taud_bl(i,k) = 0.0
+ dtaux2d(i,k) = 0.0
+ dtauy2d(i,k) = 0.0
+ enddo
+ enddo
+!
+ if (gwd_opt == 33) then
+ do i = its,im
+ dusfc_ls(i) = 0.0
+ dvsfc_ls(i) = 0.0
+ dusfc_bl(i) = 0.0
+ dvsfc_bl(i) = 0.0
+ dusfc_ss(i) = 0.0
+ dvsfc_ss(i) = 0.0
+ dusfc_fd(i) = 0.0
+ dvsfc_fd(i) = 0.0
+ enddo
+ do k = kts,km
+ do i = its,im
+ dtaux2d_ls(i,k)= 0.0
+ dtauy2d_ls(i,k)= 0.0
+ dtaux2d_bl(i,k)= 0.0
+ dtauy2d_bl(i,k)= 0.0
+ dtaux2d_ss(i,k)= 0.0
+ dtauy2d_ss(i,k)= 0.0
+ dtaux2d_fd(i,k)= 0.0
+ dtauy2d_fd(i,k)= 0.0
+ enddo
+ enddo
+ endif
+
+ do i = its,im
+ taup(i,km+1) = 0.0
+ xlinv(i) = 1.0/xl
+ dusfc(i) = 0.0
+ dvsfc(i) = 0.0
+ enddo
+!
+! initialize array for flow-blocking drag
+!
+ taufb(1:im,1:km+1) = 0.0
+ komax(1:im) = 0
+!
+ do k = kts,km
+ do i = its,im
+ vtj(i,k) = t1(i,k) * (1.+fv*q1(i,k))
+ vtk(i,k) = vtj(i,k) / prslk(i,k)
+ ro(i,k) = 1./rd * prsl(i,k) / vtj(i,k) ! density kg/m**3
+ enddo
+ enddo
+!
+! determine reference level: maximum of 2*var and pbl heights
+!
+ do i = its,im
+ zlowtop(i) = 2. * var(i)
+ enddo
+!
+ do i = its,im
+ kloop1(i) = .true.
+ enddo
+!
+ do k = kts+1,km
+ do i = its,im
+ if(kloop1(i).and.zl(i,k)-zl(i,1).ge.zlowtop(i)) then
+ klowtop(i) = k+1
+ kloop1(i) = .false.
+ endif
+ enddo
+ enddo
+!
+ do i = its,im
+ kbl(i) = max(kpbl(i), klowtop(i))
+ kbl(i) = max(min(kbl(i),kpblmax),kpblmin)
+ enddo
+!
+! determine the level of maximum orographic height
+!
+ ! komax(:) = kbl(:)
+ komax(:) = klowtop(:) - 1 ! modification by NOAA/GSD March 2018
+!
+ do i = its,im
+ delks(i) = 1.0 / (prsi(i,1) - prsi(i,kbl(i)))
+ delks1(i) = 1.0 / (prsl(i,1) - prsl(i,kbl(i)))
+ enddo
+!
+! compute low level averages within pbl
+!
+ do k = kts,kpblmax
+ do i = its,im
+ if (k.lt.kbl(i)) then
+ rcsks = rcs * del(i,k) * delks(i)
+ rdelks = del(i,k) * delks(i)
+ ubar(i) = ubar(i) + rcsks * u1(i,k) ! pbl u mean
+ vbar(i) = vbar(i) + rcsks * v1(i,k) ! pbl v mean
+ roll(i) = roll(i) + rdelks * ro(i,k) ! ro mean
+ endif
+ enddo
+ enddo
+!
+! figure out low-level horizontal wind direction
+!
+! nwd 1 2 3 4 5 6 7 8
+! wd w s sw nw e n ne se
+!
+ do i = its,im
+ wdir = atan2(ubar(i),vbar(i)) + pi
+ idir = mod(nint(fdir*wdir),mdir) + 1
+ nwd = nwdir(idir)
+ oa(i) = (1-2*int( (nwd-1)/4 )) * oa4(i,mod(nwd-1,4)+1)
+ ol(i) = ol4(i,mod(nwd-1,4)+1)
+ ! Repeat for small-scale gwd
+ oass(i) = (1-2*int( (nwd-1)/4 )) * oa4ss(i,mod(nwd-1,4)+1)
+ olss(i) = ol4ss(i,mod(nwd-1,4)+1)
+
+!
+!----- compute orographic width along (ol) and perpendicular (olp)
+!----- the direction of wind
+!
+ ol4p(1) = ol4(i,2)
+ ol4p(2) = ol4(i,1)
+ ol4p(3) = ol4(i,4)
+ ol4p(4) = ol4(i,3)
+ olp(i) = ol4p(mod(nwd-1,4)+1)
+!
+!----- compute orographic direction (horizontal orographic aspect ratio)
+!
+ od(i) = olp(i)/max(ol(i),olmin)
+ od(i) = min(od(i),odmax)
+ od(i) = max(od(i),odmin)
+!
+!----- compute length of grid in the along(dxy) and cross(dxyp) wind directions
+!
+ dxy(i) = dxy4(MOD(nwd-1,4)+1)
+ dxyp(i) = dxy4p(MOD(nwd-1,4)+1)
+ enddo
+!
+! END INITIALIZATION; BEGIN GWD CALCULATIONS:
+!
+IF ( ((gwd_opt_ls .EQ. 1).or.(gwd_opt_bl .EQ. 1)).and. &
+ (ls_taper .GT. 1.E-02) ) THEN !====
+!
+!--- saving richardson number in usqj for migwdi
+!
+ do k = kts,km-1
+ do i = its,im
+ ti = 2.0 / (t1(i,k)+t1(i,k+1))
+ rdz = 1./(zl(i,k+1) - zl(i,k))
+ tem1 = u1(i,k) - u1(i,k+1)
+ tem2 = v1(i,k) - v1(i,k+1)
+ dw2 = rcl*(tem1*tem1 + tem2*tem2)
+ shr2 = max(dw2,dw2min) * rdz * rdz
+ bvf2 = g*(g/cp+rdz*(vtj(i,k+1)-vtj(i,k))) * ti
+ usqj(i,k) = max(bvf2/shr2,rimin)
+ bnv2(i,k) = 2.0*g*rdz*(vtk(i,k+1)-vtk(i,k))/(vtk(i,k+1)+vtk(i,k))
+ bnv2(i,k) = max( bnv2(i,k), bnv2min )
+ enddo
+ enddo
+!
+!----compute the "low level" or 1/3 wind magnitude (m/s)
+!
+ do i = its,im
+ ulow(i) = max(sqrt(ubar(i)*ubar(i) + vbar(i)*vbar(i)), 1.0)
+ rulow(i) = 1./ulow(i)
+ enddo
+!
+ do k = kts,km-1
+ do i = its,im
+ velco(i,k) = (0.5*rcs) * ((u1(i,k)+u1(i,k+1)) * ubar(i) &
+ + (v1(i,k)+v1(i,k+1)) * vbar(i))
+ velco(i,k) = velco(i,k) * rulow(i)
+ if ((velco(i,k).lt.veleps) .and. (velco(i,k).gt.0.)) then
+ velco(i,k) = veleps
+ endif
+ enddo
+ enddo
+!
+! no drag when critical level in the base layer
+!
+ do i = its,im
+ ldrag(i) = velco(i,1).le.0.
+ enddo
+!
+! no drag when velco.lt.0
+!
+ do k = kpblmin,kpblmax
+ do i = its,im
+ if (k .lt. kbl(i)) ldrag(i) = ldrag(i).or. velco(i,k).le.0.
+ enddo
+ enddo
+!
+! no drag when bnv2.lt.0
+!
+ do k = kts,kpblmax
+ do i = its,im
+ if (k .lt. kbl(i)) ldrag(i) = ldrag(i).or. bnv2(i,k).lt.0.
+ enddo
+ enddo
+!
+!-----the low level weighted average ri is stored in usqj(1,1; im)
+!-----the low level weighted average n**2 is stored in bnv2(1,1; im)
+!---- this is called bnvl2 in phys_gwd_alpert_sub not bnv2
+!---- rdelks (del(k)/delks) vert ave factor so we can * instead of /
+!
+ do i = its,im
+ wtkbj = (prsl(i,1)-prsl(i,2)) * delks1(i)
+ bnv2(i,1) = wtkbj * bnv2(i,1)
+ usqj(i,1) = wtkbj * usqj(i,1)
+ enddo
+!
+ do k = kpblmin,kpblmax
+ do i = its,im
+ if (k .lt. kbl(i)) then
+ rdelks = (prsl(i,k)-prsl(i,k+1)) * delks1(i)
+ bnv2(i,1) = bnv2(i,1) + bnv2(i,k) * rdelks
+ usqj(i,1) = usqj(i,1) + usqj(i,k) * rdelks
+ endif
+ enddo
+ enddo
+!
+ do i = its,im
+ ldrag(i) = ldrag(i) .or. bnv2(i,1).le.0.0
+ ldrag(i) = ldrag(i) .or. ulow(i).eq.1.0
+ ldrag(i) = ldrag(i) .or. var(i) .le. 0.0
+ enddo
+!
+! set all ri low level values to the low level value
+!
+ do k = kpblmin,kpblmax
+ do i = its,im
+ if (k .lt. kbl(i)) usqj(i,k) = usqj(i,1)
+ enddo
+ enddo
+!
+ do i = its,im
+ if (.not.ldrag(i)) then
+ bnv(i) = sqrt( bnv2(i,1) )
+ fr(i) = bnv(i) * rulow(i) * 2. * var(i) * od(i)
+ fr(i) = min(fr(i),frmax)
+ xn(i) = ubar(i) * rulow(i)
+ yn(i) = vbar(i) * rulow(i)
+ endif
+ enddo
+!
+! compute the base level stress and store it in taub
+! calculate enhancement factor, number of mountains & aspect
+! ratio const. use simplified relationship between standard
+! deviation & critical hgt
+
+ do i = its,im
+ if (.not. ldrag(i)) then
+ efact = (oa(i) + 2.) ** (ce*fr(i)/frc)
+ efact = min( max(efact,efmin), efmax )
+!!!!!!! cleff (effective grid length) is highly tunable parameter
+!!!!!!! the bigger (smaller) value produce weaker (stronger) wave drag
+!WRF cleff = sqrt(dxy(i)**2. + dxyp(i)**2.)
+!WRF cleff = 3. * max(dx(i),cleff)
+ coefm(i) = (1. + ol(i)) ** (oa(i)+1.)
+!WRF xlinv(i) = coefm(i) / cleff
+ xlinv(i) = coefm(i) * cleff
+ tem = fr(i) * fr(i) * oc1(i)
+ gfobnv = gmax * tem / ((tem + cg)*bnv(i))
+ if ( gwd_opt_ls .NE. 0 ) then
+ taub(i) = xlinv(i) * roll(i) * ulow(i) * ulow(i) &
+ * ulow(i) * gfobnv * efact
+ else ! We've gotten what we need for the blocking scheme
+ taub(i) = 0.0
+ end if
+ else
+ taub(i) = 0.0
+ xn(i) = 0.0
+ yn(i) = 0.0
+ endif
+ enddo
+
+ENDIF ! (gwd_opt_ls .EQ. 1).or.(gwd_opt_bl .EQ. 1)
+
+!=========================================================
+! add small-scale wavedrag for stable boundary layer
+!=========================================================
+ XNBV=0.
+ tauwavex0=0.
+ tauwavey0=0.
+ density=1.2
+ utendwave=0.
+ vtendwave=0.
+ zq=0.
+!
+ IF ( (gwd_opt_ss .EQ. 1).and.(ss_taper.GT.1.E-02) ) THEN
+ if (me==master) print *,"in Drag Suite: Running small-scale gravity wave drag"
+!
+! declaring potential temperature
+!
+ do k = kts,km
+ do i = its,im
+ thx(i,k) = t1(i,k)/prslk(i,k)
+ enddo
+ enddo
+!
+ do k = kts,km
+ do i = its,im
+ tvcon = (1.+fv*q1(i,k))
+ thvx(i,k) = thx(i,k)*tvcon
+ enddo
+ enddo
+ ! Calculate mid-layer height (zl), interface height (zq), and layer depth (dz2).
+ do k = kts,km
+ do i = its,im
+ zq(i,k+1) = PHII(i,k+1)*g_inv
+ dz2(i,k) = (PHII(i,k+1)-PHII(i,k))*g_inv
+ zl(i,k) = PHIL(i,k)*g_inv
+ enddo
+ enddo
+
+ do i=its,im
+ hpbl2 = hpbl(i)+10.
+ kpbl2 = kpbl(i)
+ !kvar = MIN(kpbl, k-level of var)
+ kvar = 1
+ do k=kts+1,MAX(kpbl(i),kts+1)
+! IF (zl(i,k)>2.*var(i) .or. zl(i,k)>2*varmax) then
+ IF (zl(i,k)>300.) then
+ kpbl2 = k
+ IF (k == kpbl(i)) then
+ hpbl2 = hpbl(i)+10.
+ ELSE
+ hpbl2 = zl(i,k)+10.
+ ENDIF
+ exit
+ ENDIF
+ enddo
+ if((xland(i)-1.5).le.0. .and. 2.*varss(i).le.hpbl(i))then
+ if(br1(i).gt.0. .and. thvx(i,kpbl2)-thvx(i,kts) > 0.)then
+!WRF cleff_ss = sqrt(dxy(i)**2 + dxyp(i)**2)
+! cleff_ss = 3. * max(dx(i),cleff_ss)
+! cleff_ss = 10. * max(dxmax_ss,cleff_ss)
+!WRF cleff_ss = 0.1 * max(dxmax_ss,cleff_ss)
+ cleff_ss = 0.1 * 12000.
+ coefm_ss(i) = (1. + olss(i)) ** (oass(i)+1.)
+ xlinv(i) = coefm_ss(i) / cleff_ss
+ !govrth(i)=g/(0.5*(thvx(i,kpbl(i))+thvx(i,kts)))
+ govrth(i)=g/(0.5*(thvx(i,kpbl2)+thvx(i,kts)))
+ !XNBV=sqrt(govrth(i)*(thvx(i,kpbl(i))-thvx(i,kts))/hpbl(i))
+ XNBV=sqrt(govrth(i)*(thvx(i,kpbl2)-thvx(i,kts))/hpbl2)
+!
+ !if(abs(XNBV/u1(i,kpbl(i))).gt.xlinv(i))then
+ if(abs(XNBV/u1(i,kpbl2)).gt.xlinv(i))then
+ !tauwavex0=0.5*XNBV*xlinv(i)*(2*MIN(varss(i),75.))**2*ro(i,kts)*u1(i,kpbl(i))
+ !tauwavex0=0.5*XNBV*xlinv(i)*(2.*MIN(varss(i),40.))**2*ro(i,kts)*u1(i,kpbl2)
+ !tauwavex0=0.5*XNBV*xlinv(i)*(2.*MIN(varss(i),40.))**2*ro(i,kts)*u1(i,3)
+ var_temp = MIN(varss(i),varmax_ss) + &
+ MAX(0.,beta_ss*(varss(i)-varmax_ss))
+ ! Note: This is a semi-implicit treatment of the time differencing
+ var_temp2 = 0.5*XNBV*xlinv(i)*(2.*var_temp)**2*ro(i,kvar) ! this is greater than zero
+ tauwavex0=-var_temp2*u1(i,kvar)/(1.+var_temp2*deltim)
+ tauwavex0=tauwavex0*ss_taper
+ else
+ tauwavex0=0.
+ endif
+!
+ !if(abs(XNBV/v1(i,kpbl(i))).gt.xlinv(i))then
+ if(abs(XNBV/v1(i,kpbl2)).gt.xlinv(i))then
+ !tauwavey0=0.5*XNBV*xlinv(i)*(2*MIN(varss(i),75.))**2*ro(i,kts)*v1(i,kpbl(i))
+ !tauwavey0=0.5*XNBV*xlinv(i)*(2.*MIN(varss(i),40.))**2*ro(i,kts)*v1(i,kpbl2)
+ !tauwavey0=0.5*XNBV*xlinv(i)*(2.*MIN(varss(i),40.))**2*ro(i,kts)*v1(i,3)
+ var_temp = MIN(varss(i),varmax_ss) + &
+ MAX(0.,beta_ss*(varss(i)-varmax_ss))
+ ! Note: This is a semi-implicit treatment of the time differencing
+ tauwavey0=-var_temp2*v1(i,kvar)/(1.+var_temp2*deltim)
+ tauwavey0=tauwavey0*ss_taper
+ else
+ tauwavey0=0.
+ endif
+
+ do k=kts,kpbl(i) !MIN(kpbl2+1,km-1)
+!original
+ !utendwave(i,k)=-1.*tauwavex0*2.*max((1.-zl(i,k)/hpbl(i)),0.)/hpbl(i)
+ !vtendwave(i,k)=-1.*tauwavey0*2.*max((1.-zl(i,k)/hpbl(i)),0.)/hpbl(i)
+!new
+ utendwave(i,k)=-1.*tauwavex0*2.*max((1.-zl(i,k)/hpbl2),0.)/hpbl2
+ vtendwave(i,k)=-1.*tauwavey0*2.*max((1.-zl(i,k)/hpbl2),0.)/hpbl2
+!mod-to be used in HRRRv3/RAPv4
+ !utendwave(i,k)=-1.*tauwavex0 * max((1.-zl(i,k)/hpbl2),0.)**2
+ !vtendwave(i,k)=-1.*tauwavey0 * max((1.-zl(i,k)/hpbl2),0.)**2
+ enddo
+ endif
+ endif
+ enddo ! end i loop
+
+ do k = kts,km
+ do i = its,im
+ dudt(i,k) = dudt(i,k) + utendwave(i,k)
+ dvdt(i,k) = dvdt(i,k) + vtendwave(i,k)
+ dusfc(i) = dusfc(i) + utendwave(i,k) * del(i,k)
+ dvsfc(i) = dvsfc(i) + vtendwave(i,k) * del(i,k)
+ enddo
+ enddo
+ if (gwd_opt == 33) then
+ do k = kts,km
+ do i = its,im
+ dusfc_ss(i) = dusfc_ss(i) + utendwave(i,k) * del(i,k)
+ dvsfc_ss(i) = dvsfc_ss(i) + vtendwave(i,k) * del(i,k)
+ dtaux2d_ss(i,k) = utendwave(i,k)
+ dtauy2d_ss(i,k) = vtendwave(i,k)
+ enddo
+ enddo
+ endif
+
+ENDIF ! end if gwd_opt_ss == 1
+
+!================================================================
+! Topographic Form Drag from Beljaars et al. (2004, QJRMS, equ. 16):
+!================================================================
+IF ( (gwd_opt_fd .EQ. 1).and.(ss_taper.GT.1.E-02) ) THEN
+ if (me==master) print *,"in Drag Suite: Running form drag"
+
+ utendform=0.
+ vtendform=0.
+ zq=0.
+
+ IF ( (gwd_opt_ss .NE. 1).and.(ss_taper.GT.1.E-02) ) THEN
+ ! Defining mid-layer height (zl), interface height (zq), and layer depth (dz2).
+ ! This is already done above if the small-scale GWD is activated.
+ do k = kts,km
+ do i = its,im
+ zq(i,k+1) = PHII(i,k+1)*g_inv
+ dz2(i,k) = (PHII(i,k+1)-PHII(i,k))*g_inv
+ zl(i,k) = PHIL(i,k)*g_inv
+ enddo
+ enddo
+ ENDIF
+
+ DO i=its,im
+ IF ((xland(i)-1.5) .le. 0.) then
+ !(IH*kflt**n1)**-1 = (0.00102*0.00035**-1.9)**-1 = 0.00026615161
+ var_temp = MIN(varss(i),varmax_fd) + &
+ MAX(0.,beta_fd*(varss(i)-varmax_fd))
+ var_temp = MIN(var_temp, 250.)
+ a1=0.00026615161*var_temp**2
+! a1=0.00026615161*MIN(varss(i),varmax)**2
+! a1=0.00026615161*(0.5*varss(i))**2
+ ! k1**(n1-n2) = 0.003**(-1.9 - -2.8) = 0.003**0.9 = 0.005363
+ a2=a1*0.005363
+ ! Revise e-folding height based on PBL height and topographic std. dev. -- M. Toy 3/12/2018
+ H_efold = max(2*varss(i),hpbl(i))
+ H_efold = min(H_efold,1500.)
+ DO k=kts,km
+ wsp=SQRT(u1(i,k)**2 + v1(i,k)**2)
+ ! alpha*beta*Cmd*Ccorr*2.109 = 12.*1.*0.005*0.6*2.109 = 0.0759
+ var_temp = 0.0759*EXP(-(zl(i,k)/H_efold)**1.5)*a2* &
+ zl(i,k)**(-1.2)*ss_taper ! this is greater than zero
+ ! Note: This is a semi-implicit treatment of the time differencing
+ ! per Beljaars et al. (2004, QJRMS)
+ utendform(i,k) = - var_temp*wsp*u1(i,k)/(1. + var_temp*deltim*wsp)
+ vtendform(i,k) = - var_temp*wsp*v1(i,k)/(1. + var_temp*deltim*wsp)
+ !IF(zl(i,k) > 4000.) exit
+ ENDDO
+ ENDIF
+ ENDDO
+
+ do k = kts,km
+ do i = its,im
+ dudt(i,k) = dudt(i,k) + utendform(i,k)
+ dvdt(i,k) = dvdt(i,k) + vtendform(i,k)
+ dusfc(i) = dusfc(i) + utendform(i,k) * del(i,k)
+ dvsfc(i) = dvsfc(i) + vtendform(i,k) * del(i,k)
+ enddo
+ enddo
+ if (gwd_opt == 33) then
+ do k = kts,km
+ do i = its,im
+ dtaux2d_fd(i,k) = utendform(i,k)
+ dtauy2d_fd(i,k) = vtendform(i,k)
+ dusfc_fd(i) = dusfc_fd(i) + utendform(i,k) * del(i,k)
+ dvsfc_fd(i) = dvsfc_fd(i) + vtendform(i,k) * del(i,k)
+ enddo
+ enddo
+ endif
+
+ENDIF ! end if gwd_opt_fd == 1
+!=======================================================
+! More for the large-scale gwd component
+IF ( (gwd_opt_ls .EQ. 1).and.(ls_taper.GT.1.E-02) ) THEN
+ if (me==master) print *,"in Drag Suite: Running large-scale gravity wave drag"
+!
+! now compute vertical structure of the stress.
+ do k = kts,kpblmax
+ do i = its,im
+ if (k .le. kbl(i)) taup(i,k) = taub(i)
+ enddo
+ enddo
+!
+ do k = kpblmin, km-1 ! vertical level k loop!
+ kp1 = k + 1
+ do i = its,im
+!
+! unstablelayer if ri < ric
+! unstable layer if upper air vel comp along surf vel <=0 (crit lay)
+! at (u-c)=0. crit layer exists and bit vector should be set (.le.)
+!
+ if (k .ge. kbl(i)) then
+ icrilv(i) = icrilv(i) .or. ( usqj(i,k) .lt. ric) &
+ .or. (velco(i,k) .le. 0.0)
+ brvf(i) = max(bnv2(i,k),bnv2min) ! brunt-vaisala frequency squared
+ brvf(i) = sqrt(brvf(i)) ! brunt-vaisala frequency
+ endif
+ enddo
+!
+ do i = its,im
+ if (k .ge. kbl(i) .and. (.not. ldrag(i))) then
+ if (.not.icrilv(i) .and. taup(i,k) .gt. 0.0 ) then
+ temv = 1.0 / velco(i,k)
+ tem1 = coefm(i)/dxy(i)*(ro(i,kp1)+ro(i,k))*brvf(i)*velco(i,k)*0.5
+ hd = sqrt(taup(i,k) / tem1)
+ fro = brvf(i) * hd * temv
+!
+! rim is the minimum-richardson number by shutts (1985)
+ tem2 = sqrt(usqj(i,k))
+ tem = 1. + tem2 * fro
+ rim = usqj(i,k) * (1.-fro) / (tem * tem)
+!
+! check stability to employ the 'saturation hypothesis'
+! of lindzen (1981) except at tropospheric downstream regions
+!
+ if (rim .le. ric) then ! saturation hypothesis!
+ if ((oa(i) .le. 0.).or.(kp1 .ge. kpblmin )) then
+ temc = 2.0 + 1.0 / tem2
+ hd = velco(i,k) * (2.*sqrt(temc)-temc) / brvf(i)
+ taup(i,kp1) = tem1 * hd * hd
+ endif
+ else ! no wavebreaking!
+ taup(i,kp1) = taup(i,k)
+ endif
+ endif
+ endif
+ enddo
+ enddo
+!
+ if(lcap.lt.km) then
+ do klcap = lcapp1,km
+ do i = its,im
+ taup(i,klcap) = prsi(i,klcap) / prsi(i,lcap) * taup(i,lcap)
+ enddo
+ enddo
+ endif
+
+ENDIF !END LARGE-SCALE TAU CALCULATION
+!===============================================================
+!COMPUTE BLOCKING COMPONENT
+!===============================================================
+IF ( (gwd_opt_bl .EQ. 1) .and. (ls_taper .GT. 1.E-02) ) THEN
+ if (me==master) print *,"in Drag Suite: Running blocking drag"
+
+ do i = its,im
+ if(.not.ldrag(i)) then
+!
+!------- determine the height of flow-blocking layer
+!
+ kblk = 0
+ pe = 0.0
+ do k = km, kpblmin, -1
+ if(kblk.eq.0 .and. k.le.komax(i)) then
+ pe = pe + bnv2(i,k)*(zl(i,komax(i))-zl(i,k))*del(i,k)/g/ro(i,k)
+ ke = 0.5*((rcs*u1(i,k))**2.+(rcs*v1(i,k))**2.)
+!
+!---------- apply flow-blocking drag when pe >= ke
+!
+ if(pe.ge.ke) then
+ kblk = k
+ kblk = min(kblk,kbl(i))
+ zblk = zl(i,kblk)-zl(i,kts)
+ RDXZB(i) = real(k,kind=kind_phys)
+ endif
+ endif
+ enddo
+ if(kblk.ne.0) then
+!
+!--------- compute flow-blocking stress
+!
+ cd = max(2.0-1.0/od(i),0.0)
+ taufb(i,kts) = 0.5 * roll(i) * coefm(i) / max(dxmax_ls,dxy(i))**2 * cd * dxyp(i) &
+ * olp(i) * zblk * ulow(i)**2
+ tautem = taufb(i,kts)/float(kblk-kts)
+ do k = kts+1, kblk
+ taufb(i,k) = taufb(i,k-1) - tautem
+ enddo
+!
+!----------sum orographic GW stress and flow-blocking stress
+!
+ ! taup(i,:) = taup(i,:) + taufb(i,:) ! Keep taup and taufb separate for now
+ endif
+ endif
+ enddo
+
+ENDIF ! end blocking drag
+!===========================================================
+IF ( (gwd_opt_ls .EQ. 1 .OR. gwd_opt_bl .EQ. 1) .and. (ls_taper .GT. 1.E-02) ) THEN
+!
+! calculate - (g)*d(tau)/d(pressure) and deceleration terms dtaux, dtauy
+!
+ do k = kts,km
+ do i = its,im
+ taud_ls(i,k) = 1. * (taup(i,k+1) - taup(i,k)) * csg / del(i,k)
+ taud_bl(i,k) = 1. * (taufb(i,k+1) - taufb(i,k)) * csg / del(i,k)
+ enddo
+ enddo
+!
+! limit de-acceleration (momentum deposition ) at top to 1/2 value
+! the idea is some stuff must go out the 'top'
+ do klcap = lcap,km
+ do i = its,im
+ taud_ls(i,klcap) = taud_ls(i,klcap) * factop
+ taud_bl(i,klcap) = taud_bl(i,klcap) * factop
+ enddo
+ enddo
+!
+! if the gravity wave drag would force a critical line
+! in the lower ksmm1 layers during the next deltim timestep,
+! then only apply drag until that critical line is reached.
+!
+ do k = kts,kpblmax-1
+ do i = its,im
+ if (k .le. kbl(i)) then
+ if((taud_ls(i,k)+taud_bl(i,k)).ne.0.) &
+ dtfac(i) = min(dtfac(i),abs(velco(i,k) &
+ /(deltim*rcs*(taud_ls(i,k)+taud_bl(i,k)))))
+ endif
+ enddo
+ enddo
+!
+ do k = kts,km
+ do i = its,im
+ taud_ls(i,k) = taud_ls(i,k) * dtfac(i) * ls_taper *(1.-rstoch(i))
+ taud_bl(i,k) = taud_bl(i,k) * dtfac(i) * ls_taper *(1.-rstoch(i))
+
+ dtaux = taud_ls(i,k) * xn(i)
+ dtauy = taud_ls(i,k) * yn(i)
+ dtauxb = taud_bl(i,k) * xn(i)
+ dtauyb = taud_bl(i,k) * yn(i)
+
+ !add blocking and large-scale contributions to tendencies
+ dudt(i,k) = dtaux + dtauxb + dudt(i,k)
+ dvdt(i,k) = dtauy + dtauyb + dvdt(i,k)
+
+ if ( gsd_diss_ht_opt .EQ. 1 ) then
+ ! Calculate dissipation heating
+ ! Initial kinetic energy (at t0-dt)
+ eng0 = 0.5*( (rcs*u1(i,k))**2. + (rcs*v1(i,k))**2. )
+ ! Kinetic energy after wave-breaking/flow-blocking
+ eng1 = 0.5*( (rcs*(u1(i,k)+(dtaux+dtauxb)*deltim))**2 + &
+ (rcs*(v1(i,k)+(dtauy+dtauyb)*deltim))**2 )
+ ! Modify theta tendency
+ dtdt(i,k) = dtdt(i,k) + max((eng0-eng1),0.0)/cp/deltim/prslk(i,k)
+ end if
+
+ dusfc(i) = dusfc(i) + taud_ls(i,k)*xn(i)*del(i,k) + taud_bl(i,k)*xn(i)*del(i,k)
+ dvsfc(i) = dvsfc(i) + taud_ls(i,k)*yn(i)*del(i,k) + taud_bl(i,k)*yn(i)*del(i,k)
+ enddo
+ enddo
+
+ ! Finalize dusfc and dvsfc diagnostics
+ do i = its,im
+ dusfc(i) = (-1./g*rcs) * dusfc(i)
+ dvsfc(i) = (-1./g*rcs) * dvsfc(i)
+ enddo
+
+ if (gwd_opt == 33) then
+ do k = kts,km
+ do i = its,im
+ dtaux2d_ls(i,k) = taud_ls(i,k) * xn(i)
+ dtauy2d_ls(i,k) = taud_ls(i,k) * yn(i)
+ dtaux2d_bl(i,k) = taud_bl(i,k) * xn(i)
+ dtauy2d_bl(i,k) = taud_bl(i,k) * yn(i)
+ dusfc_ls(i) = dusfc_ls(i) + dtaux2d_ls(i,k) * del(i,k)
+ dvsfc_ls(i) = dvsfc_ls(i) + dtauy2d_ls(i,k) * del(i,k)
+ dusfc_bl(i) = dusfc_bl(i) + dtaux2d_bl(i,k) * del(i,k)
+ dvsfc_bl(i) = dvsfc_bl(i) + dtauy2d_bl(i,k) * del(i,k)
+ enddo
+ enddo
+ endif
+
+ENDIF
+
+if (gwd_opt == 33) then
+ ! Finalize dusfc and dvsfc diagnostics
+ do i = its,im
+ dusfc_ls(i) = (-1./g*rcs) * dusfc_ls(i)
+ dvsfc_ls(i) = (-1./g*rcs) * dvsfc_ls(i)
+ dusfc_bl(i) = (-1./g*rcs) * dusfc_bl(i)
+ dvsfc_bl(i) = (-1./g*rcs) * dvsfc_bl(i)
+ dusfc_ss(i) = (-1./g*rcs) * dusfc_ss(i)
+ dvsfc_ss(i) = (-1./g*rcs) * dvsfc_ss(i)
+ dusfc_fd(i) = (-1./g*rcs) * dusfc_fd(i)
+ dvsfc_fd(i) = (-1./g*rcs) * dvsfc_fd(i)
+ enddo
+endif
+!
+ return
+ end subroutine drag_suite_run
+!-------------------------------------------------------------------
+!
+!> \section arg_table_drag_suite_finalize Argument Table
+!!
+ subroutine drag_suite_finalize()
+ end subroutine drag_suite_finalize
+
+ end module drag_suite
+
+!> This module contains the CCPP-compliant orographic gravity wave drag post
+!! interstitial codes.
+ module drag_suite_post
+
+ contains
+
+!> \section arg_table_drag_suite_post_init Argument Table
+!!
+ subroutine drag_suite_post_init()
+ end subroutine drag_suite_post_init
+
+!> \section arg_table_drag_suite_post_run Argument Table
+!! \htmlinclude drag_suite_post_run.html
+!!
+ subroutine drag_suite_post_run( &
+ & lssav, ldiag3d, dtf, dusfcg, dvsfcg, dudt, dvdt, dtdt, &
+ & dugwd, dvgwd, du3dt, dv3dt, dt3dt, errmsg, errflg)
+
+ use machine, only : kind_phys
+ implicit none
+
+ logical, intent(in) :: lssav, ldiag3d
+ real(kind=kind_phys), intent(in) :: dtf
+ real(kind=kind_phys), intent(in) :: &
+ & dusfcg(:), dvsfcg(:), dudt(:,:), dvdt(:,:), dtdt(:,:)
+
+ real(kind=kind_phys), intent(inout) :: &
+ & dugwd(:), dvgwd(:), du3dt(:,:), dv3dt(:,:), dt3dt(:,:)
+
+ character(len=*), intent(out) :: errmsg
+ integer, intent(out) :: errflg
+
+ ! Initialize CCPP error handling variables
+ errmsg = ''
+ errflg = 0
+
+ if (lssav) then
+ dugwd(:) = dugwd(:) + dusfcg(:)*dtf
+ dvgwd(:) = dvgwd(:) + dvsfcg(:)*dtf
+
+ if (ldiag3d) then
+ du3dt(:,:) = du3dt(:,:) + dudt(:,:) * dtf
+ dv3dt(:,:) = dv3dt(:,:) + dvdt(:,:) * dtf
+ dt3dt(:,:) = dt3dt(:,:) + dtdt(:,:) * dtf
+ endif
+ endif
+
+ end subroutine drag_suite_post_run
+
+!> \section arg_table_drag_suite_post_finalize Argument Table
+!!
+ subroutine drag_suite_post_finalize()
+ end subroutine drag_suite_post_finalize
+
+ end module drag_suite_post
diff --git a/physics/drag_suite.meta b/physics/drag_suite.meta
new file mode 100644
index 000000000..ab84e937f
--- /dev/null
+++ b/physics/drag_suite.meta
@@ -0,0 +1,862 @@
+[ccpp-arg-table]
+ name = drag_suite_pre_init
+ type = scheme
+
+########################################################################
+[ccpp-arg-table]
+ name = drag_suite_pre_run
+ type = scheme
+[im]
+ standard_name = horizontal_loop_extent
+ long_name = horizontal dimension
+ units = count
+ dimensions = ()
+ type = integer
+ intent = in
+ optional = F
+[nmtvr]
+ standard_name = number_of_statistical_measures_of_subgrid_orography
+ long_name = number of statistical measures of subgrid orography
+ units = count
+ dimensions = ()
+ type = integer
+ intent = in
+ optional = F
+[mntvar]
+ standard_name = statistical_measures_of_subgrid_orography
+ long_name = array of statistical measures of subgrid orography
+ units = various
+ dimensions = (horizontal_dimension,number_of_statistical_measures_of_subgrid_orography)
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
+[hprime]
+ standard_name = standard_deviation_of_subgrid_orography
+ long_name = standard deviation of subgrid orography
+ units = m
+ dimensions = (horizontal_dimension)
+ type = real
+ kind = kind_phys
+ intent = out
+ optional = F
+[oc]
+ standard_name = convexity_of_subgrid_orography
+ long_name = convexity of subgrid orography
+ units = none
+ dimensions = (horizontal_dimension)
+ type = real
+ kind = kind_phys
+ intent = out
+ optional = F
+[oa4]
+ standard_name = asymmetry_of_subgrid_orography
+ long_name = asymmetry of subgrid orography
+ units = none
+ dimensions = (horizontal_dimension,4)
+ type = real
+ kind = kind_phys
+ intent = out
+ optional = F
+[clx]
+ standard_name = fraction_of_grid_box_with_subgrid_orography_higher_than_critical_height
+ long_name = horizontal fraction of grid box covered by subgrid orography higher than critical height
+ units = frac
+ dimensions = (horizontal_dimension,4)
+ type = real
+ kind = kind_phys
+ intent = out
+ optional = F
+[theta]
+ standard_name = angle_from_east_of_maximum_subgrid_orographic_variations
+ long_name = angle with_respect to east of maximum subgrid orographic variations
+ units = degrees
+ dimensions = (horizontal_dimension)
+ type = real
+ kind = kind_phys
+ intent = out
+ optional = F
+[sigma]
+ standard_name = slope_of_subgrid_orography
+ long_name = slope of subgrid orography
+ units = none
+ dimensions = (horizontal_dimension)
+ type = real
+ kind = kind_phys
+ intent = out
+ optional = F
+[gamma]
+ standard_name = anisotropy_of_subgrid_orography
+ long_name = anisotropy of subgrid orography
+ units = none
+ dimensions = (horizontal_dimension)
+ type = real
+ kind = kind_phys
+ intent = out
+ optional = F
+[elvmax]
+ standard_name = maximum_subgrid_orography
+ long_name = maximum of subgrid orography
+ units = m
+ dimensions = (horizontal_dimension)
+ type = real
+ kind = kind_phys
+ intent = out
+ optional = F
+[errmsg]
+ standard_name = ccpp_error_message
+ long_name = error message for error handling in CCPP
+ units = none
+ dimensions = ()
+ type = character
+ kind = len=*
+ intent = out
+ optional = F
+[errflg]
+ standard_name = ccpp_error_flag
+ long_name = error flag for error handling in CCPP
+ units = flag
+ dimensions = ()
+ type = integer
+ intent = out
+ optional = F
+
+########################################################################
+[ccpp-arg-table]
+ name = drag_suite_pre_finalize
+ type = scheme
+
+########################################################################
+[ccpp-arg-table]
+ name = drag_suite_init
+ type = scheme
+
+########################################################################
+[ccpp-arg-table]
+ name = drag_suite_run
+ type = scheme
+[im]
+ standard_name = horizontal_loop_extent
+ long_name = horizontal loop extent
+ units = count
+ dimensions = ()
+ type = integer
+ intent = in
+ optional = F
+[ix]
+ standard_name = horizontal_dimension
+ long_name = horizontal dimension
+ units = count
+ dimensions = ()
+ type = integer
+ intent = in
+ optional = F
+[km]
+ standard_name = vertical_dimension
+ long_name = number of vertical layers
+ units = count
+ dimensions = ()
+ type = integer
+ intent = in
+ optional = F
+[dvdt]
+ standard_name = tendency_of_y_wind_due_to_model_physics
+ long_name = meridional wind tendency due to model physics
+ units = m s-2
+ dimensions = (horizontal_dimension,vertical_dimension)
+ type = real
+ kind = kind_phys
+ intent = inout
+ optional = F
+[dudt]
+ standard_name = tendency_of_x_wind_due_to_model_physics
+ long_name = zonal wind tendency due to model physics
+ units = m s-2
+ dimensions = (horizontal_dimension,vertical_dimension)
+ type = real
+ kind = kind_phys
+ intent = inout
+ optional = F
+[dtdt]
+ standard_name = tendency_of_air_temperature_due_to_model_physics
+ long_name = air temperature tendency due to model physics
+ units = K s-1
+ dimensions = (horizontal_dimension,vertical_dimension)
+ type = real
+ kind = kind_phys
+ intent = inout
+ optional = F
+[u1]
+ standard_name = x_wind
+ long_name = zonal wind
+ units = m s-1
+ dimensions = (horizontal_dimension,vertical_dimension)
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
+[v1]
+ standard_name = y_wind
+ long_name = meridional wind
+ units = m s-1
+ dimensions = (horizontal_dimension,vertical_dimension)
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
+[t1]
+ standard_name = air_temperature
+ long_name = mid-layer temperature
+ units = K
+ dimensions = (horizontal_dimension,vertical_dimension)
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
+[q1]
+ standard_name = water_vapor_specific_humidity
+ long_name = mid-layer specific humidity of water vapor
+ units = kg kg-1
+ dimensions = (horizontal_dimension,vertical_dimension)
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
+[kpbl]
+ standard_name = vertical_index_at_top_of_atmosphere_boundary_layer
+ long_name = vertical index at top atmospheric boundary layer
+ units = index
+ dimensions = (horizontal_dimension)
+ type = integer
+ intent = in
+ optional = F
+[prsi]
+ standard_name = air_pressure_at_interface
+ long_name = interface pressure
+ units = Pa
+ dimensions = (horizontal_dimension,vertical_dimension_plus_one)
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
+[del]
+ standard_name = air_pressure_difference_between_midlayers
+ long_name = difference between mid-layer pressures
+ units = Pa
+ dimensions = (horizontal_dimension,vertical_dimension)
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
+[prsl]
+ standard_name = air_pressure
+ long_name = mid-layer pressure
+ units = Pa
+ dimensions = (horizontal_dimension,vertical_dimension)
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
+[prslk]
+ standard_name = dimensionless_exner_function_at_model_layers
+ long_name = mid-layer Exner function
+ units = none
+ dimensions = (horizontal_dimension,vertical_dimension)
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
+[phii]
+ standard_name = geopotential_at_interface
+ long_name = interface geopotential
+ units = m2 s-2
+ dimensions = (horizontal_dimension,vertical_dimension_plus_one)
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
+[phil]
+ standard_name = geopotential
+ long_name = mid-layer geopotential
+ units = m2 s-2
+ dimensions = (horizontal_dimension,vertical_dimension)
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
+[deltim]
+ standard_name = time_step_for_physics
+ long_name = physics time step
+ units = s
+ dimensions = ()
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
+[kdt]
+ standard_name = index_of_time_step
+ long_name = current time step index
+ units = index
+ dimensions = ()
+ type = integer
+ intent = in
+ optional = F
+[var]
+ standard_name = standard_deviation_of_subgrid_orography
+ long_name = standard deviation of subgrid orography
+ units = m
+ dimensions = (horizontal_dimension)
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
+[oc1]
+ standard_name = convexity_of_subgrid_orography
+ long_name = convexity of subgrid orography
+ units = none
+ dimensions = (horizontal_dimension)
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
+[oa4]
+ standard_name = asymmetry_of_subgrid_orography
+ long_name = asymmetry of subgrid orography
+ units = none
+ dimensions = (horizontal_dimension,4)
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
+[ol4]
+ standard_name = fraction_of_grid_box_with_subgrid_orography_higher_than_critical_height
+ long_name = horizontal fraction of grid box covered by subgrid orography higher than critical height
+ units = frac
+ dimensions = (horizontal_dimension,4)
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
+[theta]
+ standard_name = angle_from_east_of_maximum_subgrid_orographic_variations
+ long_name = angle with respect to east of maximum subgrid orographic variations
+ units = degrees
+ dimensions = (horizontal_dimension)
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
+[sigma]
+ standard_name = slope_of_subgrid_orography
+ long_name = slope of subgrid orography
+ units = none
+ dimensions = (horizontal_dimension)
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
+[gamma]
+ standard_name = anisotropy_of_subgrid_orography
+ long_name = anisotropy of subgrid orography
+ units = none
+ dimensions = (horizontal_dimension)
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
+[elvmax]
+ standard_name = maximum_subgrid_orography
+ long_name = maximum of subgrid orography
+ units = m
+ dimensions = (horizontal_dimension)
+ type = real
+ kind = kind_phys
+ intent = inout
+ optional = F
+[dtaux2d_ls]
+ standard_name = x_momentum_tendency_from_large_scale_gwd
+ long_name = x momentum tendency from large scale gwd
+ units = m s-2
+ dimensions = (horizontal_dimension,vertical_dimension)
+ type = real
+ kind = kind_phys
+ intent = out
+ optional = F
+[dtauy2d_ls]
+ standard_name = y_momentum_tendency_from_large_scale_gwd
+ long_name = y momentum tendency from large scale gwd
+ units = m s-2
+ dimensions = (horizontal_dimension,vertical_dimension)
+ type = real
+ kind = kind_phys
+ intent = out
+ optional = F
+[dtaux2d_bl]
+ standard_name = x_momentum_tendency_from_blocking_drag
+ long_name = x momentum tendency from blocking drag
+ units = m s-2
+ dimensions = (horizontal_dimension,vertical_dimension)
+ type = real
+ kind = kind_phys
+ intent = out
+ optional = F
+[dtauy2d_bl]
+ standard_name = y_momentum_tendency_from_blocking_drag
+ long_name = y momentum tendency from blocking drag
+ units = m s-2
+ dimensions = (horizontal_dimension,vertical_dimension)
+ type = real
+ kind = kind_phys
+ intent = out
+ optional = F
+[dtaux2d_ss]
+ standard_name = x_momentum_tendency_from_small_scale_gwd
+ long_name = x momentum tendency from small scale gwd
+ units = m s-2
+ dimensions = (horizontal_dimension,vertical_dimension)
+ type = real
+ kind = kind_phys
+ intent = out
+ optional = F
+[dtauy2d_ss]
+ standard_name = y_momentum_tendency_from_small_scale_gwd
+ long_name = y momentum tendency from small scale gwd
+ units = m s-2
+ dimensions = (horizontal_dimension,vertical_dimension)
+ type = real
+ kind = kind_phys
+ intent = out
+ optional = F
+[dtaux2d_fd]
+ standard_name = x_momentum_tendency_from_form_drag
+ long_name = x momentum tendency from form drag
+ units = m s-2
+ dimensions = (horizontal_dimension,vertical_dimension)
+ type = real
+ kind = kind_phys
+ intent = out
+ optional = F
+[dtauy2d_fd]
+ standard_name = y_momentum_tendency_from_form_drag
+ long_name = y momentum tendency from form drag
+ units = m s-2
+ dimensions = (horizontal_dimension,vertical_dimension)
+ type = real
+ kind = kind_phys
+ intent = out
+ optional = F
+[dusfc]
+ standard_name = instantaneous_x_stress_due_to_gravity_wave_drag
+ long_name = zonal surface stress due to orographic gravity wave drag
+ units = Pa
+ dimensions = (horizontal_dimension)
+ type = real
+ kind = kind_phys
+ intent = out
+ optional = F
+[dvsfc]
+ standard_name = instantaneous_y_stress_due_to_gravity_wave_drag
+ long_name = meridional surface stress due to orographic gravity wave drag
+ units = Pa
+ dimensions = (horizontal_dimension)
+ type = real
+ kind = kind_phys
+ intent = out
+ optional = F
+[dusfc_ls]
+ standard_name = integrated_x_momentum_flux_from_large_scale_gwd
+ long_name = integrated x momentum flux from large scale gwd
+ units = Pa s
+ dimensions = (horizontal_dimension)
+ type = real
+ kind = kind_phys
+ intent = out
+ optional = F
+[dvsfc_ls]
+ standard_name = integrated_y_momentum_flux_from_large_scale_gwd
+ long_name = integrated y momentum flux from large scale gwd
+ units = Pa s
+ dimensions = (horizontal_dimension)
+ type = real
+ kind = kind_phys
+ intent = out
+ optional = F
+[dusfc_bl]
+ standard_name = integrated_x_momentum_flux_from_blocking_drag
+ long_name = integrated x momentum flux from blocking drag
+ units = Pa s
+ dimensions = (horizontal_dimension)
+ type = real
+ kind = kind_phys
+ intent = out
+ optional = F
+[dvsfc_bl]
+ standard_name = integrated_y_momentum_flux_from_blocking_drag
+ long_name = integrated y momentum flux from blocking drag
+ units = Pa s
+ dimensions = (horizontal_dimension)
+ type = real
+ kind = kind_phys
+ intent = out
+ optional = F
+[dusfc_ss]
+ standard_name = integrated_x_momentum_flux_from_small_scale_gwd
+ long_name = integrated x momentum flux from small scale gwd
+ units = Pa s
+ dimensions = (horizontal_dimension)
+ type = real
+ kind = kind_phys
+ intent = out
+ optional = F
+[dvsfc_ss]
+ standard_name = integrated_y_momentum_flux_from_small_scale_gwd
+ long_name = integrated y momentum flux from small scale gwd
+ units = Pa s
+ dimensions = (horizontal_dimension)
+ type = real
+ kind = kind_phys
+ intent = out
+ optional = F
+[dusfc_fd]
+ standard_name = integrated_x_momentum_flux_from_form_drag
+ long_name = integrated x momentum flux from form drag
+ units = Pa s
+ dimensions = (horizontal_dimension)
+ type = real
+ kind = kind_phys
+ intent = out
+ optional = F
+[dvsfc_fd]
+ standard_name = integrated_y_momentum_flux_from_form_drag
+ long_name = integrated y momentum flux from form drag
+ units = Pa s
+ dimensions = (horizontal_dimension)
+ type = real
+ kind = kind_phys
+ intent = out
+ optional = F
+[slmsk]
+ standard_name = sea_land_ice_mask_real
+ long_name = landmask: sea/land/ice=0/1/2
+ units = flag
+ dimensions = (horizontal_dimension)
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
+[br1]
+ standard_name = bulk_richardson_number_at_lowest_model_level
+ long_name = bulk Richardson number at the surface
+ units = none
+ dimensions = (horizontal_dimension)
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
+[hpbl]
+ standard_name = atmosphere_boundary_layer_thickness
+ long_name = PBL thickness
+ units = m
+ dimensions = (horizontal_dimension)
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
+[g]
+ standard_name = gravitational_acceleration
+ long_name = gravitational acceleration
+ units = m s-2
+ dimensions = ()
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
+[cp]
+ standard_name = specific_heat_of_dry_air_at_constant_pressure
+ long_name = specific heat of dry air at constant pressure
+ units = J kg-1 K-1
+ dimensions = ()
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
+[rd]
+ standard_name = gas_constant_dry_air
+ long_name = ideal gas constant for dry air
+ units = J kg-1 K-1
+ dimensions = ()
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
+[rv]
+ standard_name = gas_constant_water_vapor
+ long_name = ideal gas constant for water vapor
+ units = J kg-1 K-1
+ dimensions = ()
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
+[fv]
+ standard_name = ratio_of_vapor_to_dry_air_gas_constants_minus_one
+ long_name = (rv/rd) - 1 (rv = ideal gas constant for water vapor)
+ units = none
+ dimensions = ()
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
+[pi]
+ standard_name = pi
+ long_name = ratio of a circle's circumference to its diameter
+ units = radians
+ dimensions = ()
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
+[imx]
+ standard_name = number_of_equatorial_longitude_points
+ long_name = number of longitude points along the equator
+ units = count
+ dimensions = ()
+ type = integer
+ intent = in
+ optional = F
+[cdmbgwd]
+ standard_name = multiplication_factors_for_mountain_blocking_and_orographic_gravity_wave_drag
+ long_name = multiplic. factors for (1) mountain blocking drag coeff. and (2) ref. level orographic gravity wave drag
+ units = none
+ dimensions = (2)
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
+[me]
+ standard_name = mpi_rank
+ long_name = rank of the current MPI task
+ units = index
+ dimensions = ()
+ type = integer
+ intent = in
+ optional = F
+[master]
+ standard_name = mpi_root
+ long_name = master MPI-rank
+ units = index
+ dimensions = ()
+ type = integer
+ intent = in
+ optional = F
+[lprnt]
+ standard_name = flag_print
+ long_name = flag for debugging printouts
+ units = flag
+ dimensions = ()
+ type = logical
+ intent = in
+ optional = F
+[ipr]
+ standard_name = horizontal_index_of_printed_column
+ long_name = horizontal index of column used in debugging printouts
+ units = index
+ dimensions = ()
+ type = integer
+ intent = in
+ optional = F
+[rdxzb]
+ standard_name = level_of_dividing_streamline
+ long_name = level of the dividing streamline
+ units = none
+ dimensions = (horizontal_dimension)
+ type = real
+ kind = kind_phys
+ intent = out
+ optional = F
+[dx]
+ standard_name = cell_size
+ long_name = size of the grid cell
+ units = m
+ dimensions = (horizontal_dimension)
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
+[gwd_opt]
+ standard_name = gwd_opt
+ long_name = flag to choose gwd scheme
+ units = flag
+ dimensions = ()
+ type = integer
+ intent = in
+ optional = F
+[errmsg]
+ standard_name = ccpp_error_message
+ long_name = error message for error handling in CCPP
+ units = none
+ dimensions = ()
+ type = character
+ kind = len=*
+ intent = out
+ optional = F
+[errflg]
+ standard_name = ccpp_error_flag
+ long_name = error flag for error handling in CCPP
+ units = flag
+ dimensions = ()
+ type = integer
+ intent = out
+ optional = F
+
+########################################################################
+[ccpp-arg-table]
+ name = drag_suite_finalize
+ type = scheme
+
+########################################################################
+[ccpp-arg-table]
+ name = drag_suite_post_init
+ type = scheme
+
+########################################################################
+[ccpp-arg-table]
+ name = drag_suite_post_run
+ type = scheme
+[lssav]
+ standard_name = flag_diagnostics
+ long_name = flag for calculating diagnostic fields
+ units = flag
+ dimensions = ()
+ type = logical
+ intent = in
+ optional = F
+[ldiag3d]
+ standard_name = flag_diagnostics_3D
+ long_name = flag for calculating 3-D diagnostic fields
+ units = flag
+ dimensions = ()
+ type = logical
+ intent = in
+ optional = F
+[dtf]
+ standard_name = time_step_for_dynamics
+ long_name = dynamics time step
+ units = s
+ dimensions = ()
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
+[dusfcg]
+ standard_name = instantaneous_x_stress_due_to_gravity_wave_drag
+ long_name = zonal surface stress due to orographic gravity wave drag
+ units = Pa
+ dimensions = (horizontal_dimension)
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
+[dvsfcg]
+ standard_name = instantaneous_y_stress_due_to_gravity_wave_drag
+ long_name = meridional surface stress due to orographic gravity wave drag
+ units = Pa
+ dimensions = (horizontal_dimension)
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
+[dudt]
+ standard_name = tendency_of_x_wind_due_to_model_physics
+ long_name = zonal wind tendency due to model physics
+ units = m s-2
+ dimensions = (horizontal_dimension,vertical_dimension)
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
+[dvdt]
+ standard_name = tendency_of_y_wind_due_to_model_physics
+ long_name = meridional wind tendency due to model physics
+ units = m s-2
+ dimensions = (horizontal_dimension,vertical_dimension)
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
+[dtdt]
+ standard_name = tendency_of_air_temperature_due_to_model_physics
+ long_name = air temperature tendency due to model physics
+ units = K s-1
+ dimensions = (horizontal_dimension,vertical_dimension)
+ type = real
+ kind = kind_phys
+ intent = in
+ optional = F
+[dugwd]
+ standard_name = time_integral_of_x_stress_due_to_gravity_wave_drag
+ long_name = integral over time of zonal stress due to gravity wave drag
+ units = Pa s
+ dimensions = (horizontal_dimension)
+ type = real
+ kind = kind_phys
+ intent = inout
+ optional = F
+[dvgwd]
+ standard_name = time_integral_of_y_stress_due_to_gravity_wave_drag
+ long_name = integral over time of meridional stress due to gravity wave drag
+ units = Pa s
+ dimensions = (horizontal_dimension)
+ type = real
+ kind = kind_phys
+ intent = inout
+ optional = F
+[du3dt]
+ standard_name = cumulative_change_in_x_wind_due_to_orographic_gravity_wave_drag
+ long_name = cumulative change in zonal wind due to orographic gravity wave drag
+ units = m s-1
+ dimensions = (horizontal_dimension,vertical_dimension)
+ type = real
+ kind = kind_phys
+ intent = inout
+ optional = F
+[dv3dt]
+ standard_name = cumulative_change_in_y_wind_due_to_orographic_gravity_wave_drag
+ long_name = cumulative change in meridional wind due to orographic gravity wave drag
+ units = m s-1
+ dimensions = (horizontal_dimension,vertical_dimension)
+ type = real
+ kind = kind_phys
+ intent = inout
+ optional = F
+[dt3dt]
+ standard_name = cumulative_change_in_temperature_due_to_orographic_gravity_wave_drag
+ long_name = cumulative change in temperature due to orographic gravity wave drag
+ units = K
+ dimensions = (horizontal_dimension,vertical_dimension)
+ type = real
+ kind = kind_phys
+ intent = inout
+ optional = F
+[errmsg]
+ standard_name = ccpp_error_message
+ long_name = error message for error handling in CCPP
+ units = none
+ dimensions = ()
+ type = character
+ kind = len=*
+ intent = out
+ optional = F
+[errflg]
+ standard_name = ccpp_error_flag
+ long_name = error flag for error handling in CCPP
+ units = flag
+ dimensions = ()
+ type = integer
+ intent = out
+ optional = F
+
+########################################################################
+[ccpp-arg-table]
+ name = drag_suite_post_finalize
+ type = scheme
diff --git a/physics/gfdl_cloud_microphys.F90 b/physics/gfdl_cloud_microphys.F90
index 6d907e40a..1ccedb956 100644
--- a/physics/gfdl_cloud_microphys.F90
+++ b/physics/gfdl_cloud_microphys.F90
@@ -113,7 +113,7 @@ end subroutine gfdl_cloud_microphys_finalize
!! \htmlinclude gfdl_cloud_microphys_run.html
!!
subroutine gfdl_cloud_microphys_run( &
- levs, im, con_g, con_fvirt, con_rd, frland, garea, &
+ levs, im, con_g, con_fvirt, con_rd, frland, garea, islmsk, &
gq0, gq0_ntcw, gq0_ntrw, gq0_ntiw, gq0_ntsw, gq0_ntgl, gq0_ntclamt, &
gt0, gu0, gv0, vvl, prsl, phii, del, &
rain0, ice0, snow0, graupel0, prcp0, sr, &
@@ -136,6 +136,7 @@ subroutine gfdl_cloud_microphys_run( &
integer, intent(in ) :: levs, im
real(kind=kind_phys), intent(in ) :: con_g, con_fvirt, con_rd
real(kind=kind_phys), intent(in ), dimension(1:im) :: frland, garea
+ integer, intent(in ), dimension(1:im) :: islmsk
real(kind=kind_phys), intent(inout), dimension(1:im,1:levs) :: gq0, gq0_ntcw, gq0_ntrw, gq0_ntiw, &
gq0_ntsw, gq0_ntgl, gq0_ntclamt
real(kind=kind_phys), intent(inout), dimension(1:im,1:levs) :: gt0, gu0, gv0
@@ -170,9 +171,6 @@ subroutine gfdl_cloud_microphys_run( &
real(kind=kind_phys), dimension(:,:), allocatable :: den
real(kind=kind_phys) :: onebg
real(kind=kind_phys) :: tem
-#ifdef TRANSITION
- real(kind=kind_phys), volatile :: volatile_var1, volatile_var2
-#endif
! Initialize CCPP error handling variables
errmsg = ''
@@ -260,18 +258,10 @@ subroutine gfdl_cloud_microphys_run( &
! calculate fraction of frozen precipitation using unscaled
! values of rain0, ice0, snow0, graupel0 (for bit-for-bit)
do i=1,im
-#ifdef TRANSITION
- volatile_var1 = rain0(i)+snow0(i)+ice0(i)+graupel0(i)
- volatile_var2 = snow0(i)+ice0(i)+graupel0(i)
- prcp0(i) = volatile_var1 * tem
- if ( volatile_var1 * tem > rainmin ) then
- sr(i) = volatile_var2 / volatile_var1
-#else
prcp0(i) = (rain0(i)+snow0(i)+ice0(i)+graupel0(i)) * tem
if ( prcp0(i) > rainmin ) then
sr(i) = (snow0(i) + ice0(i) + graupel0(i)) &
/ (rain0(i) + snow0(i) + ice0(i) + graupel0(i))
-#endif
else
sr(i) = 0.0
endif
@@ -309,9 +299,11 @@ subroutine gfdl_cloud_microphys_run( &
enddo
enddo
call cloud_diagnosis (1, im, 1, levs, den(1:im,1:levs), &
+ del(1:im,1:levs), islmsk(1:im), &
gq0_ntcw(1:im,1:levs), gq0_ntiw(1:im,1:levs), &
- gq0_ntrw(1:im,1:levs), gq0_ntsw(1:im,1:levs), &
- gq0_ntgl(1:im,1:levs), gt0(1:im,1:levs), &
+ gq0_ntrw(1:im,1:levs), &
+ gq0_ntsw(1:im,1:levs) + gq0_ntgl(1:im,1:levs), &
+ gq0_ntgl(1:im,1:levs)*0.0, gt0(1:im,1:levs), &
rew(1:im,1:levs), rei(1:im,1:levs), rer(1:im,1:levs),&
res(1:im,1:levs), reg(1:im,1:levs))
deallocate(den)
diff --git a/physics/gfdl_cloud_microphys.meta b/physics/gfdl_cloud_microphys.meta
index c2ce3f8f5..7f31637bf 100644
--- a/physics/gfdl_cloud_microphys.meta
+++ b/physics/gfdl_cloud_microphys.meta
@@ -180,6 +180,14 @@
kind = kind_phys
intent = in
optional = F
+[islmsk]
+ standard_name = sea_land_ice_mask
+ long_name = sea/land/ice mask (=0/1/2)
+ units = flag
+ dimensions = (horizontal_dimension)
+ type = integer
+ intent = in
+ optional = F
[gq0]
standard_name = water_vapor_specific_humidity_updated_by_physics
long_name = water vapor specific humidity updated by physics
diff --git a/physics/gfdl_fv_sat_adj.F90 b/physics/gfdl_fv_sat_adj.F90
index d535ebc91..f5c84cd99 100644
--- a/physics/gfdl_fv_sat_adj.F90
+++ b/physics/gfdl_fv_sat_adj.F90
@@ -49,7 +49,7 @@ module fv_sat_adj
! | gfdl_cloud_microphys_mod |
! ql_gen, qi_gen, qi0_max, ql_mlt, ql0_max, qi_lim, qs_mlt,
! tau_r2g, tau_smlt, tau_i2s, tau_v2l, tau_l2v, tau_imlt, tau_l2r,
-! rad_rain, rad_snow, rad_graupel, dw_ocean, dw_land |
+! rad_rain, rad_snow, rad_graupel, dw_ocean, dw_land, tintqs
!
!
! DH* TODO - MAKE THIS INPUT ARGUMENTS *DH
@@ -64,8 +64,7 @@ module fv_sat_adj
use gfdl_cloud_microphys_mod, only: ql_gen, qi_gen, qi0_max, ql_mlt, ql0_max, qi_lim, qs_mlt
use gfdl_cloud_microphys_mod, only: icloud_f, sat_adj0, t_sub, cld_min
use gfdl_cloud_microphys_mod, only: tau_r2g, tau_smlt, tau_i2s, tau_v2l, tau_l2v, tau_imlt, tau_l2r
- use gfdl_cloud_microphys_mod, only: rad_rain, rad_snow, rad_graupel, dw_ocean, dw_land
-
+ use gfdl_cloud_microphys_mod, only: rad_rain, rad_snow, rad_graupel, dw_ocean, dw_land, tintqs
#ifdef MULTI_GASES
use ccpp_multi_gases_mod, only: multi_gases_init, &
multi_gases_finalize, &
@@ -296,10 +295,6 @@ subroutine fv_sat_adj_run(mdt, zvir, is, ie, isd, ied, kmp, km, kmdelz, js, je,
! Local variables
real(kind=kind_dyn), dimension(is:ie,js:je) :: dpln
-#ifdef TRANSITION
- ! For bit-for-bit reproducibility
- real(kind=kind_dyn), volatile :: volatile_var
-#endif
integer :: kdelz
integer :: k, j, i
@@ -317,9 +312,6 @@ subroutine fv_sat_adj_run(mdt, zvir, is, ie, isd, ied, kmp, km, kmdelz, js, je,
!$OMP ql,qv,te0,fast_mp_consv, &
!$OMP hydrostatic,ng,zvir,pkz, &
!$OMP akap,te0_2d,ngas,qvi) &
-#ifdef TRANSITION
-!$OMP private(volatile_var) &
-#endif
!$OMP private(k,j,i,kdelz,dpln)
#endif
@@ -351,28 +343,13 @@ subroutine fv_sat_adj_run(mdt, zvir, is, ie, isd, ied, kmp, km, kmdelz, js, je,
do j=js,je
do i=is,ie
#ifdef MOIST_CAPPA
-#ifdef TRANSITION
- volatile_var = log(rrg*delp(i,j,k)/delz(i,j,k)*pt(i,j,k))
- pkz(i,j,k) = exp(cappa(i,j,k)*volatile_var)
-#else
pkz(i,j,k) = exp(cappa(i,j,k)*log(rrg*delp(i,j,k)/delz(i,j,k)*pt(i,j,k)))
-#endif
-#else
-#ifdef TRANSITION
-#ifdef MULTI_GASES
- volatile_var = log(rrg*delp(i,j,k)/delz(i,j,k)*pt(i,j,k))
- pkz(i,j,k) = exp(akap*(virqd(q(i,j,k,1:num_gas))/vicpqd(q(i,j,k,1:num_gas))*volatile_var)
-#else
- volatile_var = log(rrg*delp(i,j,k)/delz(i,j,k)*pt(i,j,k))
- pkz(i,j,k) = exp(akap*volatile_var)
-#endif
#else
#ifdef MULTI_GASES
pkz(i,j,k) = exp(akap*(virqd(q(i,j,k,1:num_gas))/vicpqd(q(i,j,k,1:num_gas))*log(rrg*delp(i,j,k)/delz(i,j,k)*pt(i,j,k)))
#else
pkz(i,j,k) = exp(akap*log(rrg*delp(i,j,k)/delz(i,j,k)*pt(i,j,k)))
#endif
-#endif
#endif
enddo
enddo
@@ -1052,9 +1029,13 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te,
do i = is, ie
+ if(tintqs) then
+ tin = pt1(i)
+ else
tin = pt1 (i) - (lcp2 (i) * q_cond (i) + icp2 (i) * q_sol (i)) ! minimum temperature
! tin = pt1 (i) - ((lv00 + d0_vap * pt1 (i)) * q_cond (i) + &
! (li00 + dc_ice * pt1 (i)) * q_sol (i)) / (mc_air (i) + qpz (i) * c_vap)
+ endif
! -----------------------------------------------------------------------
! determine saturated specific humidity
@@ -1097,14 +1078,14 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te,
! icloud_f = 2: binary cloud scheme (0 / 1)
! -----------------------------------------------------------------------
- if (rh > 0.75 .and. qpz (i) > 1.e-6) then
+ if (rh > 0.75 .and. qpz (i) > 1.e-8) then
dq = hvar (i) * qpz (i)
q_plus = qpz (i) + dq
q_minus = qpz (i) - dq
if (icloud_f == 2) then
if (qpz (i) > qstar (i)) then
qa (i, j) = 1.
- elseif (qstar (i) < q_plus .and. q_cond (i) > 1.e-6) then
+ elseif (qstar (i) < q_plus .and. q_cond (i) > 1.e-8) then
qa (i, j) = ((q_plus - qstar (i)) / dq) ** 2
qa (i, j) = min (1., qa (i, j))
else
@@ -1124,7 +1105,7 @@ subroutine fv_sat_adj_work(mdt, zvir, is, ie, js, je, ng, hydrostatic, consv_te,
qa (i, j) = 0.
endif
! impose minimum cloudiness if substantial q_cond (i) exist
- if (q_cond (i) > 1.e-6) then
+ if (q_cond (i) > 1.e-8) then
qa (i, j) = max (cld_min, qa (i, j))
endif
qa (i, j) = min (1., qa (i, j))
diff --git a/physics/gwdc.f b/physics/gwdc.f
index 80898c47b..9909a3100 100644
--- a/physics/gwdc.f
+++ b/physics/gwdc.f
@@ -22,7 +22,7 @@ end subroutine gwdc_pre_init
subroutine gwdc_pre_run ( &
& im, cgwf, dx, work1, work2, dlength, cldf, &
& levs, kbot, ktop, dtp, gt0, gt0_init, del, cumabs, &
- & errmsg, errflg )
+ & do_cnvgwd, errmsg, errflg )
use machine, only : kind_phys
implicit none
@@ -38,6 +38,7 @@ subroutine gwdc_pre_run ( &
real(kind=kind_phys), intent(out) :: &
& dlength(:), cldf(:), cumabs(:)
+ logical, intent(in) :: do_cnvgwd
character(len=*), intent(out) :: errmsg
integer, intent(out) :: errflg
@@ -49,6 +50,14 @@ subroutine gwdc_pre_run ( &
errmsg = ''
errflg = 0
+ ! DH*
+ if (.not. do_cnvgwd) then
+ write(0,*) "ERROR: , GWDC_PRE CALLED BUT DO_CNVGWD FALSE"
+ call sleep(5)
+ stop
+ end if
+ ! *DH
+
do i = 1, im
tem1 = dx(i)
tem2 = tem1
diff --git a/physics/gwdc.meta b/physics/gwdc.meta
index b87529aec..2151cc5f7 100644
--- a/physics/gwdc.meta
+++ b/physics/gwdc.meta
@@ -137,6 +137,14 @@
kind = kind_phys
intent = out
optional = F
+[do_cnvgwd]
+ standard_name = flag_for_convective_gravity_wave_drag
+ long_name = flag for convective gravity wave drag (gwd)
+ units = flag
+ dimensions = ()
+ type = logical
+ intent = in
+ optional = F
[errmsg]
standard_name = ccpp_error_message
long_name = error message for error handling in CCPP
diff --git a/physics/gwdps.f b/physics/gwdps.f
index 366a8b974..0ea2c8754 100644
--- a/physics/gwdps.f
+++ b/physics/gwdps.f
@@ -299,12 +299,8 @@ subroutine gwdps_run( &
! Interface variables
integer, intent(in) :: im, ix, km, imx, kdt, ipr, me
integer, intent(in) :: KPBL(IM) ! Index for the PBL top layer!
- ! DH* adding intent(in) information for the following variables
- ! changes the results on Theia/Intel - skip for bit-for-bit results *DH
-! real(kind=kind_phys), intent(in) :: &
-! & deltim, G, CP, RD, RV, cdmbgwd(2)
- real(kind=kind_phys) deltim, G, CP, RD, RV, cdmbgwd(2)
- ! *DH
+ real(kind=kind_phys), intent(in) :: &
+ & deltim, G, CP, RD, RV, cdmbgwd(4)
real(kind=kind_phys), intent(inout) :: &
& A(IX,KM), B(IX,KM), C(IX,KM)
real(kind=kind_phys), intent(in) :: &
@@ -382,7 +378,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)
+ &, DTFAC(IM), XLINV(IM), DELKS(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) &
@@ -392,7 +389,8 @@ subroutine gwdps_run( &
! real(kind=kind_phys) VELKO(KM-1)
integer kref(IM), kint(im), iwk(im), ipt(im)
! for lm mtn blocking
- integer kreflm(IM), iwklm(im)
+ integer iwklm(im)
+! integer kreflm(IM), iwklm(im)
integer idxzb(im), ktrial, klevm1
!
real(kind=kind_phys) gor, gocp, fv, gr2, bnv, fr &
@@ -470,7 +468,7 @@ subroutine gwdps_run( &
do i=1,npt
iwklm(i) = 2
IDXZB(i) = 0
- kreflm(i) = 0
+! kreflm(i) = 0
enddo
! if (lprnt)
! & print *,' in gwdps_lm.f npt,IM,IX,IY,km,me=',npt,IM,IX,IY,km,me
@@ -552,14 +550,14 @@ subroutine gwdps_run( &
!
DO I = 1, npt
J = ipt(i)
- DELKS(I) = 1.0 / (PRSI(J,1) - PRSI(J,iwklm(i)))
- DELKS1(I) = 1.0 / (PRSL(J,1) - PRSL(J,iwklm(i)))
- UBAR (I) = 0.0
- VBAR (I) = 0.0
- ROLL (I) = 0.0
- PE (I) = 0.0
- EK (I) = 0.0
- BNV2bar(I) = (PRSL(J,1)-PRSL(J,2)) * DELKS1(I) * BNV2LM(I,1)
+ DELKS(I) = 1.0 / (PRSI(J,1) - PRSI(J,iwklm(i)))
+! DELKS1(I) = 1.0 / (PRSI(J,1) - PRSL(J,iwklm(i)))
+ UBAR (I) = 0.0
+ VBAR (I) = 0.0
+ ROLL (I) = 0.0
+ PE (I) = 0.0
+ EK (I) = 0.0
+ BNV2bar(I) = (PRSI(J,1)-PRSL(J,1)) * DELKS(I) * BNV2LM(I,1)
ENDDO
! --- find the dividing stream line height
@@ -567,13 +565,13 @@ subroutine gwdps_run( &
! --- iwklm(i) is the k-index of mtn elvmax elevation
!> - Find the dividing streamline height starting from the level above
!! the maximum mountain height and processing downward.
- DO Ktrial = KMLL, 1, -1
- DO I = 1, npt
- IF ( Ktrial < iwklm(I) .and. kreflm(I) == 0 ) then
- kreflm(I) = Ktrial
- ENDIF
- ENDDO
- ENDDO
+! DO Ktrial = KMLL, 1, -1
+! DO I = 1, npt
+! IF ( Ktrial < iwklm(I) .and. kreflm(I) == 0 ) then
+! kreflm(I) = Ktrial
+! ENDIF
+! ENDDO
+! ENDDO
! print *,' in gwdps_lm.f 4 npt=',npt,kreflm(npt),me
!
! --- in the layer kreflm(I) to 1 find PE (which needs N, ELVMAX)
@@ -582,13 +580,17 @@ subroutine gwdps_run( &
! --- is the vert ave of quantities from the surface to mtn top.
!
DO I = 1, npt
- DO K = 1, Kreflm(I)
+ DO K = 1, iwklm(i)-1
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
- RDELKS = (PRSL(J,K)-PRSL(J,K+1)) * DELKS1(I)
+ if (k < iwklm(I)-1) then
+ RDELKS = (PRSL(J,K)-PRSL(J,K+1)) * DELKS(I)
+ else
+ RDELKS = (PRSL(J,K)-PRSI(J,K+1)) * DELKS(I)
+ endif
BNV2bar(I) = BNV2bar(I) + BNV2lm(I,K) * RDELKS
! --- these vert ave are for diags, testing and GWD to follow (*j*).
ENDDO
@@ -862,14 +864,14 @@ subroutine gwdps_run( &
J = ipt(i)
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)))
+! DELKS1(I) = 1.0 / (PRSI(J,1) - PRSL(J,kref(I)))
UBAR (I) = 0.0
VBAR (I) = 0.0
ROLL (I) = 0.0
KBPS = MAX(KBPS, kref(I))
KMPS = MIN(KMPS, kref(I))
!
- BNV2bar(I) = (PRSL(J,1)-PRSL(J,2)) * DELKS1(I) * BNV2(I,1)
+ BNV2bar(I) = (PRSI(J,1)-PRSL(J,1)) * DELKS(I) * BNV2(I,1)
ENDDO
! print *,' in gwdps_lm.f GWD:15 =',KBPS,KMPS
KBPSP1 = KBPS + 1
@@ -883,7 +885,11 @@ subroutine gwdps_run( &
VBAR(I) = VBAR(I) + RDELKS * V1(J,K) ! Mean V below kref
!
ROLL(I) = ROLL(I) + RDELKS * RO(I,K) ! Mean RO below kref
- RDELKS = (PRSL(J,K)-PRSL(J,K+1)) * DELKS1(I)
+ if (k < kref(i)-1) then
+ RDELKS = (PRSL(J,K)-PRSL(J,K+1)) * DELKS(I)
+ else
+ RDELKS = (PRSL(J,K)-PRSI(J,K+1)) * DELKS(I)
+ endif
BNV2bar(I) = BNV2bar(I) + BNV2(I,K) * RDELKS
ENDIF
ENDDO
diff --git a/physics/gwdps.meta b/physics/gwdps.meta
index 97b6abae3..0a141b208 100644
--- a/physics/gwdps.meta
+++ b/physics/gwdps.meta
@@ -318,7 +318,7 @@
standard_name = multiplication_factors_for_mountain_blocking_and_orographic_gravity_wave_drag
long_name = multiplic. factors for (1) mountain blocking drag coeff. and (2) ref. level orographic gravity wave drag
units = none
- dimensions = (2)
+ dimensions = (4)
type = real
kind = kind_phys
intent = in
diff --git a/physics/machine.F b/physics/machine.F
index f50a950d7..896b665da 100644
--- a/physics/machine.F
+++ b/physics/machine.F
@@ -1,10 +1,8 @@
module machine
-#if 0
!! \section arg_table_machine
!! \htmlinclude machine.html
!!
-#endif
implicit none
@@ -20,6 +18,7 @@ module machine
&, kind_phys = 8 ,kind_taum=8 &
&, kind_grid = 8 &
&, kind_REAL = 8 &! used in cmp_comm
+ &, kind_LOGICAL = 4 &
&, kind_INTEGER = 4 ! -,,-
#else
@@ -34,6 +33,7 @@ module machine
&, kind_phys = 4 ,kind_taum=4 &
&, kind_grid = 4 &
&, kind_REAL = 4 &! used in cmp_comm
+ &, kind_LOGICAL = 4 &
&, kind_INTEGER = 4 ! -,,-
#endif
diff --git a/physics/machine.meta b/physics/machine.meta
index 43fc8770a..d93f50e09 100644
--- a/physics/machine.meta
+++ b/physics/machine.meta
@@ -19,3 +19,15 @@
units = none
dimensions = ()
type = integer
+[kind_LOGICAL]
+ standard_name = kind_LOGICAL
+ long_name = definition of kind_LOGICAL
+ units = none
+ dimensions = ()
+ type = integer
+[kind_INTEGER]
+ standard_name = kind_INTEGER
+ long_name = definition of kind_INTEGER
+ units = none
+ dimensions = ()
+ type = integer
diff --git a/physics/module_BL_MYJPBL.F90 b/physics/module_BL_MYJPBL.F90
new file mode 100755
index 000000000..b23e67cb5
--- /dev/null
+++ b/physics/module_BL_MYJPBL.F90
@@ -0,0 +1,2188 @@
+!-----------------------------------------------------------------------
+!
+ MODULE MODULE_BL_MYJPBL
+!
+!-----------------------------------------------------------------------
+!
+!*** THE MYJ PBL SCHEME
+!
+!-----------------------------------------------------------------------
+!
+! USE MODULE_INCLUDE
+!
+! USE MODULE_CONSTANTS,ONLY : A2,A3,A4,CP,ELIV,ELWV,ELIWV &
+! ,EP_1,EPSQ &
+! ,G,P608,PI,PQ0,R_D,R_V,RHOWATER &
+! ,STBOLT,CAPPA
+
+ USE machine, only: kfpt => kind_phys, &
+ kint => kind_INTEGER, &
+ klog => kind_LOGICAL
+
+!-----------------------------------------------------------------------
+!
+ IMPLICIT NONE
+!
+!-----------------------------------------------------------------------
+! integer,parameter :: isingle=selected_int_kind(r=9)
+! integer,parameter :: idouble=selected_int_kind(r=18)
+! integer,parameter :: single=selected_real_kind(p=6,r=37)
+! integer,parameter :: double=selected_real_kind(p=13,r=200)
+
+! integer,parameter:: &
+! klog=4 &
+! ,kint=isingle &
+! ,kdin=idouble &
+! ,kfpt=single &
+! ,kdbl=double
+
+! real (kind=kfpt),parameter :: r4_in=x'ffbfffff'
+! real (kind=kdbl),parameter :: r8_in=x'fff7ffffffffffff'
+! integer(kind=kint),parameter :: i4_in=-999 ! -huge(1)
+
+ ! integer,parameter:: &
+ ! klog=4 & ! logical variables
+ ! ,kint=4 & ! integer variables
+ ! !,kfpt=4 & ! floating point variables
+ ! ,kfpt=8 & ! floating point variables
+ ! ,kdbl=8 ! double precision
+
+ REAL(kind=kfpt),PARAMETER :: A2=17.2693882,A3=273.15,A4=35.86,CP=1004.6 &
+ ,ELIV=2.850e6,ELWV=2.501e6,R_V=461.6 &
+! ,EPSQ=1.e-12,EPSQ2=0.02,G=9.8060226 &
+ ,EPSQ=1.e-12,G=9.8060226 &
+ ,PQ0=379.90516,R_D=287.04,EP_1=R_V/R_D-1. &
+ ,P608=R_V/R_D-1.,PI=3.141592653589793 &
+ ,RHOWATER=1000.,STBOLT=5.67051E-8,CAPPA=R_D/CP
+ REAL(kind=kfpt),PARAMETER :: eliwv=2.683e6
+!
+ REAL(kind=kfpt),PARAMETER :: CONW=1./G,CONT=CP/G,CONQ=ELWV/G
+
+!-----------------------------------------------------------------------
+!
+ PRIVATE
+!
+ PUBLIC:: MYJPBL_INIT, MYJPBL
+!
+!-----------------------------------------------------------------------
+!-----------------------------------------------------------------------
+!*** FOR MYJ TURBULENCE
+!-----------------------------------------------------------------------
+!-----------------------------------------------------------------------
+!
+ REAL(KIND=KFPT),PARAMETER:: &
+ ELEVFC=0.6
+!
+ REAL(KIND=KFPT),PARAMETER:: &
+ VKARMAN=0.4 &
+!
+ ,XLS=ELIV,XLV=ELWV &
+ ,RLIVWV=XLS/XLV,ELOCP=2.72E6/CP &
+!
+ ,EPS1=1.E-12,EPS2=0. &
+ ,EPSRU=1.E-7,EPSRS=1.E-7 &
+ ,EPSTRB=1.E-24 &
+ ,FH=1.10 &
+!
+ ,ALPH=0.30,BETA=1./273.,EL0MAX=1000.,EL0MIN=1. &
+! ,ELFC=0.5,GAM1=0.2222222222222222222 &
+! ,ELFC=0.23*0.25,GAM1=0.2222222222222222222 &
+ ,ELFC=1.,GAM1=0.2222222222222222222 &
+!
+ ,A1=0.659888514560862645 &
+ ,A2X=0.6574209922667784586 &
+ ,B1=11.87799326209552761 &
+ ,B2=7.226971804046074028 &
+ ,C1=0.000830955950095854396 &
+ ,ELZ0=0.,ESQ=5.0 &
+!
+ ,SEAFC=0.98,PQ0SEA=PQ0*SEAFC &
+!
+ ,BTG=BETA*G &
+ ,ESQHF=0.5*5.0 &
+ ,RB1=1./B1
+!
+ REAL(KIND=KFPT),PARAMETER:: &
+ ADNH= 9.*A1*A2X*A2X*(12.*A1+3.*B2)*BTG*BTG &
+ ,ADNM=18.*A1*A1*A2X*(B2-3.*A2X)*BTG &
+ ,ANMH=-9.*A1*A2X*A2X*BTG*BTG &
+ ,ANMM=-3.*A1*A2X*(3.*A2X+3.*B2*C1+18.*A1*C1-B2)*BTG &
+ ,BDNH= 3.*A2X*(7.*A1+B2)*BTG &
+ ,BDNM= 6.*A1*A1 &
+ ,BEQH= A2X*B1*BTG+3.*A2X*(7.*A1+B2)*BTG &
+ ,BEQM=-A1*B1*(1.-3.*C1)+6.*A1*A1 &
+ ,BNMH=-A2X*BTG &
+ ,BNMM=A1*(1.-3.*C1) &
+ ,BSHH=9.*A1*A2X*A2X*BTG &
+ ,BSHM=18.*A1*A1*A2X*C1 &
+ ,BSMH=-3.*A1*A2X*(3.*A2X+3.*B2*C1+12.*A1*C1-B2)*BTG &
+ ,CESH=A2X &
+ ,CESM=A1*(1.-3.*C1) &
+ ,CNV=EP_1*G/BTG
+!
+!-----------------------------------------------------------------------
+!*** FREE TERM IN THE EQUILIBRIUM EQUATION FOR (L/Q)**2
+!-----------------------------------------------------------------------
+!
+ REAL(KIND=KFPT),PARAMETER:: &
+ AEQH=9.*A1*A2X*A2X*B1*BTG*BTG &
+ +9.*A1*A2X*A2X*(12.*A1+3.*B2)*BTG*BTG &
+ ,AEQM=3.*A1*A2X*B1*(3.*A2X+3.*B2*C1+18.*A1*C1-B2) &
+ *BTG+18.*A1*A1*A2X*(B2-3.*A2X)*BTG
+!
+!-----------------------------------------------------------------------
+!*** FORBIDDEN TURBULENCE AREA
+!-----------------------------------------------------------------------
+!
+ REAL(KIND=KFPT),PARAMETER:: &
+ REQU=-AEQH/AEQM &
+ ,EPSGH=1.E-9,EPSGM=REQU*EPSGH
+!
+!-----------------------------------------------------------------------
+!*** NEAR ISOTROPY FOR SHEAR TURBULENCE, WW/Q2 LOWER LIMIT
+!-----------------------------------------------------------------------
+!
+ REAL(KIND=KFPT),PARAMETER:: &
+ UBRYL=(18.*REQU*A1*A1*A2X*B2*C1*BTG &
+ +9.*A1*A2X*A2X*B2*BTG*BTG) &
+ /(REQU*ADNM+ADNH) &
+ ,UBRY=(1.+EPSRS)*UBRYL,UBRY3=3.*UBRY
+!
+ REAL(KIND=KFPT),PARAMETER:: &
+ AUBH=27.*A1*A2X*A2X*B2*BTG*BTG-ADNH*UBRY3 &
+ ,AUBM=54.*A1*A1*A2X*B2*C1*BTG -ADNM*UBRY3 &
+ ,BUBH=(9.*A1*A2X+3.*A2X*B2)*BTG-BDNH*UBRY3 &
+ ,BUBM=18.*A1*A1*C1 -BDNM*UBRY3 &
+ ,CUBR=1. - UBRY3 &
+ ,RCUBR=1./CUBR
+!
+!-----------------------------------------------------------------------
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+!---LOOK-UP TABLES------------------------------------------------------
+INTEGER(KIND=KINT),PARAMETER:: &
+ ITBL=401 & ! CONVECTION TABLES, DIMENSION 1
+,JTBL=1201 & ! CONVECTION TABLES, DIMENSION 2
+,KERFM=301 & ! SIZE OF ERF HALF TABLE
+,KERFM2=KERFM-2 ! INTERNAL POINTS OF ERF HALF TABLE
+
+REAL(KIND=KFPT),PARAMETER:: &
+ PL=2500. & ! LOWER BOUND OF PRESSURE RANGE
+,PH=105000. & ! UPPER BOUND OF PRESSURE RANGE
+,THL=210. & ! LOWER BOUND OF POTENTIAL TEMPERATURE RANGE
+,THH=365. & ! UPPER BOUND OF POTENTIAL TEMPERATURE RANGE
+,XEMIN=0. & ! LOWER BOUND OF ERF HALF TABLE
+,XEMAX=3. ! UPPER BOUND OF ERF HALF TABLE
+
+REAL(KIND=KFPT),PRIVATE,SAVE:: &
+ RDP & ! SCALING FACTOR FOR PRESSURE
+,RDQ & ! SCALING FACTOR FOR HUMIDITY
+,RDTH & ! SCALING FACTOR FOR POTENTIAL TEMPERATURE
+,RDTHE & ! SCALING FACTOR FOR EQUIVALENT POT. TEMPERATURE
+,RDXE ! ERF HALF TABLE SCALING FACTOR
+
+REAL(KIND=KFPT),DIMENSION(1:ITBL),PRIVATE,SAVE:: &
+ STHE & ! RANGE FOR EQUIVALENT POTENTIAL TEMPERATURE
+,THE0 ! BASE FOR EQUIVALENT POTENTIAL TEMPERATURE
+
+REAL(KIND=KFPT),DIMENSION(1:JTBL),PRIVATE,SAVE:: &
+ QS0 & ! BASE FOR SATURATION SPECIFIC HUMIDITY
+,SQS ! RANGE FOR SATURATION SPECIFIC HUMIDITY
+
+REAL(KIND=KFPT),DIMENSION(1:KERFM),PRIVATE,SAVE:: &
+ HERFF ! HALF ERF TABLE
+
+REAL(KIND=KFPT),DIMENSION(1:ITBL,1:JTBL),PRIVATE,SAVE:: &
+ PTBL ! SATURATION PRESSURE TABLE
+
+REAL(KIND=KFPT),DIMENSION(1:JTBL,1:ITBL),PRIVATE,SAVE:: &
+ TTBL ! TEMPERATURE TABLE
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+!-----------------------------------------------------------------------
+!
+ CONTAINS
+!
+!-----------------------------------------------------------------------
+!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
+!-----------------------------------------------------------------------
+!
+! REFERENCES: JANJIC (2001), NCEP OFFICE NOTE 437
+!
+! ABSTRACT:
+! MYJ UPDATES THE TURBULENT KINETIC ENERGY WITH THE PRODUCTION/
+! DISSIPATION TERM AND THE VERTICAL DIFFUSION TERM
+! (USING AN IMPLICIT FORMULATION) FROM MELLOR-YAMADA
+! LEVEL 2.5 AS EXTENDED BY JANJIC. EXCHANGE COEFFICIENTS FOR
+! THE SURFACE LAYER ARE COMPUTED FROM THE MONIN-OBUKHOV THEORY.
+! THE TURBULENT VERTICAL EXCHANGE IS THEN EXECUTED.
+!
+!-----------------------------------------------------------------------
+ SUBROUTINE MYJPBL(NTSD,ME,DT_PHS,EPSL,EPSQ2,HT,STDH,DZ,DEL &
+ ,PMID,PINH,TH,T,EXNER,Q,CWM,U,V &
+ ,TSK,QSFC,CHKLOWQ,THZ0,QZ0,UZ0,VZ0 &
+ ,XLAND,SICE,SNOW &
+ ,Q2,EXCH_H,USTAR,Z0,EL_MYJ,PBLH,KPBL,CT &
+ ,AKHS,AKMS,ELFLX,MIXHT,THLM,QLM &
+ ,RUBLTEN,RVBLTEN,RTHBLTEN,RQBLTEN,RQCBLTEN &
+ ,DUSFC,DVSFC,DTSFC,DQSFC,xkzo,xkzmo,ICT &
+ ,IDS,IDE,JDS,JDE &
+ ,IMS,IME,JMS,JME &
+ ,ITS,ITE,JTS,JTE,LM)
+
+! SUBROUTINE MYJPBL(DT,NPHS,EPSL,EPSQ2,HT,STDH,DZ &
+! ,PMID,PINH,TH,T,EXNER,Q,CWM,U,V &
+! ,TSK,QSFC,CHKLOWQ,THZ0,QZ0,UZ0,VZ0 &
+! ,XLAND,SICE,SNOW &
+! ,Q2,EXCH_H,USTAR,Z0,EL_MYJ,PBLH,KPBL,CT &
+! ,AKHS,AKMS,ELFLX,MIXHT &
+! ,RUBLTEN,RVBLTEN,RTHBLTEN,RQBLTEN,RQCBLTEN &
+! ,IDS,IDE,JDS,JDE &
+! ,IMS,IME,JMS,JME &
+! ,ITS,ITE,JTS,JTE,LM)
+
+!----------------------------------------------------------------------
+!
+ IMPLICIT NONE
+!
+ logical(kind=klog),save:: &
+ reinit
+!----------------------------------------------------------------------
+ INTEGER(KIND=KINT),INTENT(IN):: &
+ IDS,IDE,JDS,JDE &
+ ,IMS,IME,JMS,JME &
+ ,ITS,ITE,JTS,JTE,LM
+!
+ INTEGER,INTENT(IN) :: ICT,ME,NTSD
+
+! INTEGER(KIND=KINT),INTENT(IN):: &
+! NPHS
+!
+ INTEGER(KIND=KINT),DIMENSION(IMS:IME,JMS:JME),INTENT(OUT):: &
+ KPBL
+!
+ REAL(KIND=KFPT),INTENT(IN):: &
+ DT_PHS
+! DT
+!
+ real(kind=kfpt),dimension(1:lm-1),intent(inout):: EPSL
+ real(kind=kfpt),dimension(1:lm),intent(in):: EPSQ2
+!
+ REAL(KIND=KFPT),DIMENSION(IMS:IME,JMS:JME),INTENT(IN):: &
+ HT,SICE,SNOW,STDH &
+ ,TSK,XLAND &
+ ,CHKLOWQ,ELFLX,THLM,QLM
+!
+ REAL(KIND=KFPT),DIMENSION(IMS:IME,JMS:JME,1:LM),INTENT(IN):: &
+ DZ,EXNER,PMID,Q,CWM,U,V,T,TH,DEL,xkzo,xkzmo
+!
+ REAL(KIND=KFPT),DIMENSION(IMS:IME,JMS:JME,1:LM+1),INTENT(IN):: &
+ PINH
+!
+ REAL(KIND=KFPT),DIMENSION(IMS:IME,JMS:JME),INTENT(OUT):: &
+ MIXHT &
+ ,PBLH
+!
+ REAL(KIND=KFPT),DIMENSION(IMS:IME,JMS:JME,1:LM),INTENT(OUT):: &
+ EL_MYJ
+!
+ REAL(KIND=KFPT),DIMENSION(IMS:IME,JMS:JME,1:LM),INTENT(OUT):: &
+ RQCBLTEN &
+ ,RUBLTEN,RVBLTEN &
+ ,RTHBLTEN,RQBLTEN
+!
+ REAL(kind=KFPT),DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: &
+ DUSFC,DVSFC &
+ ,DTSFC,DQSFC
+!
+ REAL(KIND=KFPT),DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT):: &
+ AKHS,AKMS
+!
+ REAL(KIND=KFPT),DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT):: &
+ CT,QSFC,QZ0 &
+ ,THZ0,USTAR &
+ ,UZ0,VZ0,Z0
+!
+ REAL(KIND=KFPT),DIMENSION(IMS:IME,JMS:JME,1:LM),INTENT(INOUT):: &
+ EXCH_H &
+ ,Q2
+!
+!----------------------------------------------------------------------
+!***
+!*** LOCAL VARIABLES
+!***
+ INTEGER(KIND=KINT):: &
+ I,IQTB,ITTB,J,K,LLOW,LMH,LMXL
+!
+ INTEGER(KIND=KINT),DIMENSION(IMS:IME,JMS:JME):: &
+ LPBL
+!
+ REAL(KIND=KFPT):: &
+ AKHS_DENS,AKMS_DENS,BQ,BQS00K,BQS10K &
+ ,DCDT,DELTAZ,DQDT,DTDIF,DTDT,DTTURBL &
+ ,P00K,P01K,P10K,P11K,PELEVFC,PP1,PSFC,PSP,PTOP &
+ ,QBT,QFC1,QLOW,QQ1,QX &
+ ,RDTTURBL,RG,RSQDT,RXNERS,RXNSFC &
+ ,SEAMASK,SQ,SQS00K,SQS10K &
+ ,THBT,THNEW,THOLD,TQ,TTH &
+ ,ULOW,VLOW,RSTDH,STDFAC,ZSF,ZSX,ZSY,ZUV
+!
+ REAL(KIND=KFPT),DIMENSION(1:LM):: &
+ CWMK,PK,PSK,Q2K,QK,RHOK,RXNERK,THEK,THK,THVK,TK,UK,VK
+!
+ REAL(KIND=KFPT),DIMENSION(1:LM-1):: &
+ AKHK,AKMK,DCOL,EL,GH,GM
+!
+ REAL(KIND=KFPT),DIMENSION(1:LM+1):: &
+ ZHK
+!
+ REAL(KIND=KFPT),DIMENSION(IMS:IME,JMS:JME):: &
+ THSK
+!
+ REAL(KIND=KFPT),DIMENSION(IMS:IME,JMS:JME,1:LM):: &
+ RXNER,THV
+!
+ REAL(KIND=KFPT),DIMENSION(IMS:IME,JMS:JME,1:LM-1):: &
+ AKH,AKM
+!
+ REAL(KIND=KFPT),DIMENSION(IMS:IME,JMS:JME,1:LM+1):: &
+ ZINT
+!
+!*** Begin debugging
+ REAL(KIND=KFPT):: ZSL_DIAG
+ INTEGER(KIND=KINT):: IMD,JMD,PRINT_DIAG
+!*** End debugging
+!-----------------------------------------------------------------------
+!***********************************************************************
+ data reinit/.false./
+!-----------------------------------------------------------------------
+! if(reinit) then
+! call MYJPBL_INIT( &
+! 1,IDE,1,1,LM, &
+! 1,IDE,1,1, &
+! 1,IDE,1,1)
+! reinit=.false.
+! endif
+!
+!----------------------------------------------------------------------
+!**********************************************************************
+!----------------------------------------------------------------------
+!
+!*** Begin debugging
+ IMD=(IMS+IME)/2
+ JMD=(JMS+JME)/2
+!*** End debugging
+!
+!*** MAKE PREPARATIONS
+!
+!----------------------------------------------------------------------
+ STDFAC=1.
+!----------------------------------------------------------------------
+! DTTURBL=DT*NPHS
+ DTTURBL=DT_PHS
+ RDTTURBL=1./DTTURBL
+ RSQDT=SQRT(RDTTURBL)
+ DTDIF=DTTURBL
+ RG=1./G
+!
+ DO K=1,LM-1
+ DO J=JTS,JTE
+ DO I=ITS,ITE
+ AKM(I,J,K)=0.
+ ENDDO
+ ENDDO
+ ENDDO
+!
+ DO K=1,LM+1
+ DO J=JTS,JTE
+ DO I=ITS,ITE
+ ZINT(I,J,K)=0.
+ ENDDO
+ ENDDO
+ ENDDO
+!
+ DO J=JTS,JTE
+ DO I=ITS,ITE
+ ZINT(I,J,LM+1)=HT(I,J) ! Z AT BOTTOM OF LOWEST SIGMA LAYER
+ ENDDO
+ ENDDO
+!
+ DO K=LM,1,-1
+ DO J=JTS,JTE
+ DO I=ITS,ITE
+ ZINT(I,J,K)=ZINT(I,J,K+1)+DZ(I,J,K)
+ RXNER(I,J,K)=1./EXNER(I,J,K)
+ THV(I,J,K)=(Q(I,J,K)*0.608+(1.-CWM(I,J,K)))*TH(I,J,K)
+ ENDDO
+ ENDDO
+ ENDDO
+!
+ DO J=JTS,JTE
+ DO I=ITS,ITE
+ EL_MYJ(I,J,LM)=0.
+ ENDDO
+ ENDDO
+ DO J=JTS,JTE
+ DO I=ITS,ITE
+ DUSFC(I,J)=0.
+ DVSFC(I,J)=0.
+ DTSFC(I,J)=0.
+ DQSFC(I,J)=0.
+ ENDDO
+ ENDDO
+
+!
+!----------------------------------------------------------------------
+!.......................................................................
+!ZJ$OMP PARALLEL DO &
+!ZJ$OMP PRIVATE(J,I,LMH,PTOP,PSFC,SEAMASK,K,TK,THVK,QK,Q2K,RXNERK, &
+!ZJ$OMP PK,UK,VK,Q2K,ZHK,LMXL,GM,GH,EL,AKMK,AKHK,DELTAZ), &
+!ZJ$OMP SCHEDULE(DYNAMIC)
+!.......................................................................
+!----------------------------------------------------------------------
+ setup_integration: DO J=JTS,JTE
+!----------------------------------------------------------------------
+!
+ DO I=ITS,ITE
+!
+ LMH=LM
+!
+ PTOP=PINH(I,J,1)
+ PSFC=PINH(I,J,LMH+1)
+!
+!*** CONVERT LAND MASK (1 FOR SEA; 0 FOR LAND)
+!
+ SEAMASK=XLAND(I,J)-1.
+!
+!*** FILL 1-D VERTICAL ARRAYS
+!
+ DO K=LM,1,-1
+ PK(K)=PMID(I,J,K)
+ TK(K)=T(I,J,K)
+ QK(K)=Q(I,J,K)
+ THVK(K)=THV(I,J,K)
+ RXNERK(K)=RXNER(I,J,K)
+ UK(K)=U(I,J,K)
+ VK(K)=V(I,J,K)
+ Q2K(K)=Q2(I,J,K)
+!
+!*** COMPUTE THE HEIGHTS OF THE LAYER INTERFACES
+!
+ ZHK(K)=ZINT(I,J,K)
+!
+ ENDDO
+ ZHK(LM+1)=HT(I,J) ! Z AT BOTTOM OF LOWEST SIGMA LAYER
+!
+!*** POTENTIAL INSTABILITY
+!
+ PELEVFC=PMID(I,J,LMH)*ELEVFC
+!
+ DO K=LMH,1,-1
+!-----------------------------------------------------------------------
+ IF(K==LMH .OR. PMID(I,J,K)>PELEVFC) THEN
+!---PREPARATION FOR SEARCH FOR MAX CAPE---------------------------------
+ QBT=QK(K)
+ THBT=TH(I,J,K)
+ TTH=(THBT-THL)*RDTH
+ QQ1=TTH-AINT(TTH)
+ ITTB=INT(TTH)+1
+!---KEEPING INDICES WITHIN THE TABLE------------------------------------
+ IF(ITTB.LT.1)THEN
+ ITTB=1
+ QQ1=0.
+ ELSE IF(ITTB.GE.JTBL)THEN
+ ITTB=JTBL-1
+ QQ1=0.
+ ENDIF
+!---BASE AND SCALING FACTOR FOR SPEC. HUMIDITY--------------------------
+ BQS00K=QS0(ITTB)
+ SQS00K=SQS(ITTB)
+ BQS10K=QS0(ITTB+1)
+ SQS10K=SQS(ITTB+1)
+!--------------SCALING SPEC. HUMIDITY & TABLE INDEX---------------------
+ BQ=(BQS10K-BQS00K)*QQ1+BQS00K
+ SQ=(SQS10K-SQS00K)*QQ1+SQS00K
+ TQ=(QBT-BQ)/SQ*RDQ
+ PP1=TQ-AINT(TQ)
+ IQTB=INT(TQ)+1
+!----------------KEEPING INDICES WITHIN THE TABLE-----------------------
+ IF(IQTB.LT.1)THEN
+ IQTB=1
+ PP1=0.
+ ELSEIF(IQTB.GE.ITBL)THEN
+ IQTB=ITBL-1
+ PP1=0.
+ ENDIF
+!--------------SATURATION PRESSURE AT FOUR SURROUNDING TABLE PTS.-------
+ P00K=PTBL(IQTB ,ITTB )
+ P10K=PTBL(IQTB+1,ITTB )
+ P01K=PTBL(IQTB ,ITTB+1)
+ P11K=PTBL(IQTB+1,ITTB+1)
+!--------------SATURATION POINT VARIABLES AT THE BOTTOM-----------------
+ PSP=P00K+(P10K-P00K)*PP1+(P01K-P00K)*QQ1 &
+ +(P00K-P10K-P01K+P11K)*PP1*QQ1
+ RXNERS=(1.E5/PSP)**CAPPA
+ THEK(K)=THBT*EXP(ELOCP*QBT*RXNERS/THBT)
+ PSK (K)=PSP
+!-----------------------------------------------------------------------
+ ELSE
+!-----------------------------------------------------------------------
+ THEK(K)=THEK(K+1)
+ PSK (K)=PINH(I,J,1)
+!-----------------------------------------------------------------------
+ ENDIF
+!-----------------------------------------------------------------------
+ ENDDO
+!
+!*** Begin debugging
+! IF(I==IMD.AND.J==JMD)THEN
+! PRINT_DIAG=1
+! ELSE
+! PRINT_DIAG=0
+! ENDIF
+! IF(I==227.AND.J==363)PRINT_DIAG=2
+!*** End debugging
+!
+!----------------------------------------------------------------------
+!***
+!*** FIND THE MIXING LENGTH
+!***
+ CALL MIXLEN(LMH,RSQDT,UK,VK,THVK,THEK &
+ ,Q2K,EPSL,EPSQ2,ZHK,PK,PSK,RXNERK,GM,GH,EL &
+ ,PBLH(I,J),LPBL(I,J),LMXL,CT(I,J),MIXHT(I,J) &
+ ,I,J,LM)
+!
+!----------------------------------------------------------------------
+!***
+!*** SOLVE FOR THE PRODUCTION/DISSIPATION OF
+!*** THE TURBULENT KINETIC ENERGY
+!***
+!
+ CALL PRODQ2(NTSD,ME,LMH,DTTURBL,USTAR(I,J),GM,GH,EL,Q2K &
+ ,EPSL,EPSQ2,I,J,LM)
+
+! if(i.eq.4)print*,'11ql test Q2(LMH)=',Q2K(LMH),B1,USTAR(I,J)
+!
+!----------------------------------------------------------------------
+!*** THE MODEL LAYER (COUNTING UPWARD) CONTAINING THE TOP OF THE PBL
+!----------------------------------------------------------------------
+!
+ KPBL(I,J)=LPBL(I,J)
+!
+!----------------------------------------------------------------------
+!***
+!*** FIND THE EXCHANGE COEFFICIENTS IN THE FREE ATMOSPHERE
+!***
+ CALL DIFCOF(NTSD,ME,LMH,LMXL,GM,GH,EL,TK,Q2K,ZHK,AKMK,AKHK,I,J,LM &
+ ,PRINT_DIAG,KPBL(I,J))
+!
+!*** COUNTING DOWNWARD FROM THE TOP, THE EXCHANGE COEFFICIENTS AKH
+!*** ARE DEFINED ON THE BOTTOMS OF THE LAYERS 1 TO LM-1. COUNTING
+!*** COUNTING UPWARD FROM THE BOTTOM, THOSE SAME COEFFICIENTS EXCH_H
+!*** ARE DEFINED ON THE TOPS OF THE LAYERS 1 TO LM-1.
+!
+ DO K=1,LM-1
+
+ DELTAZ=0.5*(ZHK(K)-ZHK(K+2))
+ AKHK(K)=max(AKHK(K),xkzo(I,J,K)/DELTAZ) ! add minimum background diffusion
+ AKMK(K)=max(AKMK(K),xkzmo(I,J,K)/DELTAZ)
+ if((THVK(LM)-THVK(K)).GT.0.) then
+ AKHK(K)=max(AKHK(K),3./DELTAZ) ! add minimum background diffusion
+ AKMK(K)=max(AKMK(K),3./DELTAZ)
+ end if
+ AKH(I,J,K)=AKHK(K)
+ AKM(I,J,K)=AKMK(K)
+ EXCH_H(I,J,K)=AKHK(K)*DELTAZ
+ ENDDO
+!
+!----------------------------------------------------------------------
+!***
+!*** CARRY OUT THE VERTICAL DIFFUSION OF
+!*** TURBULENT KINETIC ENERGY
+!***
+!
+ CALL VDIFQ(LMH,DTDIF,Q2K,EL,ZHK,I,J,LM)
+!
+!*** SAVE THE NEW Q2 AND MIXING LENGTH.
+!
+ DO K=1,LM
+ Q2(I,J,K)=MAX(Q2K(K),EPSQ2(K))
+ IF(K