From b7a35311940736efe39de9c62f22e3a28b024f4e Mon Sep 17 00:00:00 2001 From: Grant Firl Date: Fri, 13 Dec 2019 11:44:17 -0700 Subject: [PATCH 1/6] add preliminary satmedmfvdifq documentation --- physics/docs/ccpp_doxyfile | 6 + physics/docs/pdftxt/GFS_SATMEDMFVDIFQ.txt | 35 ++++ physics/docs/pdftxt/all_shemes_list.txt | 1 + physics/mfpbltq.f | 2 +- physics/mfscuq.f | 2 +- physics/satmedmfvdifq.F | 239 +++++++++++++++------- physics/tridi.f | 3 + 7 files changed, 213 insertions(+), 75 deletions(-) create mode 100644 physics/docs/pdftxt/GFS_SATMEDMFVDIFQ.txt diff --git a/physics/docs/ccpp_doxyfile b/physics/docs/ccpp_doxyfile index e4b2e0501..339ddb3f8 100644 --- a/physics/docs/ccpp_doxyfile +++ b/physics/docs/ccpp_doxyfile @@ -113,6 +113,7 @@ INPUT = pdftxt/mainpage.txt \ pdftxt/GFS_SFCSICE.txt \ pdftxt/GFS_HEDMF.txt \ pdftxt/GFS_SATMEDMF.txt \ + pdftxt/GFS_SATMEDMFVDIFQ.txt \ pdftxt/GFS_GWDPS.txt \ pdftxt/GFS_OZPHYS.txt \ pdftxt/GFS_H2OPHYS.txt \ @@ -189,6 +190,11 @@ INPUT = pdftxt/mainpage.txt \ ../mfpblt.f \ ../mfscu.f \ ../tridi.f \ +### satmedmfvdifq + ../satmedmfvdifq.F \ + ../mfpbltq.f \ + ../mfscuq.f \ + ../tridi.f \ ### Orographic Gravity Wave ../gwdps.f \ ### Rayleigh Dampling diff --git a/physics/docs/pdftxt/GFS_SATMEDMFVDIFQ.txt b/physics/docs/pdftxt/GFS_SATMEDMFVDIFQ.txt new file mode 100644 index 000000000..de543fe6c --- /dev/null +++ b/physics/docs/pdftxt/GFS_SATMEDMFVDIFQ.txt @@ -0,0 +1,35 @@ +/** +\page GFS_SATMEDMFVDIFQ GFS Scale-aware TKE-based Moist Eddy-Diffusion Mass-Flux (EDMF) PBL and Free Atmospheric Turbulence Scheme +\section des_satmedmfvdifq Description + +The current operational \ref GFS_HEDMF uses a hybrid EDMF parameterization for the convective PBL (Han et al. 2016 \cite Han_2016; +Han et al. 2017 \cite han_et_al_2017), where the EDMF scheme is applied only for the strongly unstable PBL, while the eddy-diffusivity +counter-gradient(EDCG) scheme is used for the weakly unstable PBL. The new TKE-EDMF is an extended version of \ref GFS_HEDMF with below enhancement: + +-# Eddy diffusivity (K) is now a function of TKE which is prognostically predicted + +-# EDMF approach is applied for all the unstable PBL + +-# EDMF approach is also applied to the stratocumulus-top-driven turbulence mixing + +-# It includes a moist-adiabatic process when updraft thermal becomes saturated + +-# Scale-aware capability + +-# It includes interaction between TKE and cumulus convection + +The CCPP-compliant subroutine satmedmfvdifq_run() computes subgrid vertical turbulence mixing using scale-aware +TKE-based moist eddy-diffusion mass-flux paramterization (Han et al. 2019 \cite Han_2019) +- For the convective boundary layer, the scheme adopts EDMF parameterization (Siebesma et al. (2007)\cite Siebesma_2007) +to take into account nonlocal transport by large eddies(mfpbltq.f) +- A new mass-flux paramterization for stratocumulus-top-induced turbulence mixing has been introduced (mfscuq.f; previously, +it was an eddy diffusion form) +- For local turbulence mixing, a TKE closure model is used. + +\section intra_satmedmfvdifq Intraphysics Communication +\ref arg_table_satmedmfvdifq_run + +\section gen_pbl_satmedmfvdifq General Algorithm +\ref gen_satmedmfvdifq + +*/ diff --git a/physics/docs/pdftxt/all_shemes_list.txt b/physics/docs/pdftxt/all_shemes_list.txt index 3f2290d7b..7e5e3298e 100644 --- a/physics/docs/pdftxt/all_shemes_list.txt +++ b/physics/docs/pdftxt/all_shemes_list.txt @@ -14,6 +14,7 @@ parameterizations in suites. - \b PBL \b and \b Turbulence - \subpage GFS_HEDMF - \subpage GFS_SATMEDMF + - \subpage GFS_SATMEDMFVDIFQ - \subpage GSD_MYNNEDMF - \b Land \b Surface \b Model diff --git a/physics/mfpbltq.f b/physics/mfpbltq.f index 0f4004444..a6fc22cef 100644 --- a/physics/mfpbltq.f +++ b/physics/mfpbltq.f @@ -3,7 +3,7 @@ !! updraft parcel properties for thermals driven by surface heating !! for use in the TKE-EDMF PBL scheme (updated version). -!>\ingroup satmedmfq +!>\ingroup satmedmfvdifq !! This subroutine computes mass flux and updraft parcel properties for !! thermals driven by surface heating. !!\section mfpbltq_gen GFS mfpblt General Algorithm diff --git a/physics/mfscuq.f b/physics/mfscuq.f index c6f66b74b..3390c3e58 100644 --- a/physics/mfscuq.f +++ b/physics/mfscuq.f @@ -2,7 +2,7 @@ !! This file contains the mass flux and downdraft parcel preperties !! parameterization for stratocumulus-top-driven turbulence (updated version). -!>\ingroup satmedmfq +!>\ingroup satmedmfvdifq !! This subroutine computes mass flux and downdraft parcel properties !! for stratocumulus-top-driven turbulence. !! \section mfscuq GFS mfscu General Algorithm diff --git a/physics/satmedmfvdifq.F b/physics/satmedmfvdifq.F index 546cefca6..8a93cc5fa 100644 --- a/physics/satmedmfvdifq.F +++ b/physics/satmedmfvdifq.F @@ -7,6 +7,15 @@ module satmedmfvdifq contains +!> \defgroup satmedmfvdifq GFS Scale-aware TKE-based Moist Eddy-Diffusivity Mass-flux (TKE-EDMF, updated version) Scheme Module +!! @{ +!! \brief This subroutine contains all of the logic for the +!! scale-aware TKE-based moist eddy-diffusion mass-flux (TKE-EDMF, updated version) scheme. +!! For local turbulence mixing, a TKE closure model is used. +!! Updated version of satmedmfvdif.f (May 2019) to have better low level +!! inversion, to reduce the cold bias in lower troposphere, +!! and to reduce the negative wind speed bias in upper troposphere + !> \section arg_table_satmedmfvdifq_init Argument Table !! \htmlinclude satmedmfvdifq_init.html !! @@ -33,30 +42,21 @@ end subroutine satmedmfvdifq_init subroutine satmedmfvdifq_finalize () end subroutine satmedmfvdifq_finalize -!> \defgroup satmedmfq GFS Scale-aware TKE-based Moist Eddy-Diffusivity Mass-flux (TKE-EDMF, updated version) Scheme Module -!! @{ -!! \brief This subroutine contains all of the logic for the -!! scale-aware TKE-based moist eddy-diffusion mass-flux (TKE-EDMF, updated version) scheme. -!! !> \section arg_table_satmedmfvdifq_run Argument Table !! \htmlinclude satmedmfvdifq_run.html !! -!!\section gen_satmedmfvdif GFS satmedmfvdif General Algorithm -!! satmedmfvdif_run() computes subgrid vertical turbulence mixing +!!\section gen_satmedmfvdifq GFS satmedmfvdifq General Algorithm +!! satmedmfvdifq_run() computes subgrid vertical turbulence mixing !! using the scale-aware TKE-based moist eddy-diffusion mass-flux (EDMF) parameterization of !! Han and Bretherton (2019) \cite Han_2019 . !! -# The local turbulent mixing is represented by an eddy-diffusivity scheme which !! is a function of a prognostic TKE. !! -# For the convective boundary layer, nonlocal transport by large eddies -!! (mfpblt.f), is represented using a mass flux approach (Siebesma et al.(2007) \cite Siebesma_2007 ). +!! (mfpbltq.f), is represented using a mass flux approach (Siebesma et al.(2007) \cite Siebesma_2007 ). !! -# A mass-flux approach is also used to represent the stratocumulus-top-induced turbulence -!! (mfscu.f). -!! For local turbulence mixing, a TKE closure model is used. -!! Updated version of satmedmfvdif.f (May 2019) to have better low level -!! inversion, to reduce the cold bias in lower troposphere, -!! and to reduce the negative wind speed bias in upper troposphere +!! (mfscuq.f). !! \section detail_satmedmfvidfq GFS satmedmfvdifq Detailed Algorithm -!> @{ +!! @{ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & & grav,rd,cp,rv,hvap,hfus,fv,eps,epsm1, & & dv,du,tdt,rtg,u1,v1,t1,q1,swh,hlw,xmu,garea, & @@ -241,6 +241,7 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & errmsg = '' errflg = 0 +!> ## Compute preliminary variables from input arguments dt2 = delt rdt = 1. / dt2 ! @@ -251,7 +252,8 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & km1 = km - 1 kmpbl = km / 2 kmscu = km / 2 -! +!> - Compute physical height of the layer centers and interfaces from +!! the geopotential height (\p zi and \p zl) do k=1,km do i=1,im zi(i,k) = phii(i,k) * gravi @@ -276,11 +278,12 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & zm(i,k) = zi(i,k+1) enddo enddo -! horizontal grid size +!> - Compute horizontal grid size (\p gdx) do i=1,im gdx(i) = sqrt(garea(i)) enddo -! +!> - Initialize tke value at vertical layer centers and interfaces +!! from tracer (\p tke and \p tkeh) do k=1,km do i=1,im tke(i,k) = max(q1(i,k,ntke), tkmin) @@ -291,7 +294,7 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & tkeh(i,k) = 0.5 * (tke(i,k) + tke(i,k+1)) enddo enddo -! +!> - Compute reciprocal of \f$ \Delta z \f$ (rdzt) do k = 1,km1 do i=1,im rdzt(i,k) = 1.0 / (zl(i,k+1) - zl(i,k)) @@ -299,12 +302,18 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & enddo enddo ! -! set background diffusivities as a function of -! horizontal grid size with xkzm_h & xkzm_m for gdx >= 25km -! and 0.01 for gdx=5m, i.e., -! xkzm_hx = 0.01 + (xkzm_h - 0.01)/(xkgdx-5.) * (gdx-5.) -! xkzm_mx = 0.01 + (xkzm_h - 0.01)/(xkgdx-5.) * (gdx-5.) -! +!> - Compute reciprocal of pressure (tx1, tx2) + +!> - Compute minimum turbulent mixing length (rlmnz) + +!> - Compute background vertical diffusivities for scalars and momentum (xkzo and xkzmo) + +!> - set background diffusivities as a function of +!! horizontal grid size with xkzm_h & xkzm_m for gdx >= 25km +!! and 0.01 for gdx=5m, i.e., +!! \n xkzm_hx = 0.01 + (xkzm_h - 0.01)/(xkgdx-5.) * (gdx-5.) +!! \n xkzm_mx = 0.01 + (xkzm_h - 0.01)/(xkgdx-5.) * (gdx-5.) + do i=1,im kx1(i) = 1 tx1(i) = 1.0 / prsi(i,1) @@ -352,7 +361,8 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & endif enddo enddo -! + +!> - Some output variables and logical flags are initialized do i = 1,im z0(i) = 0.01 * zorl(i) dusfc(i) = 0. @@ -376,7 +386,9 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & kcld(i) = km1 endif enddo -! + +!> - Compute \f$\theta\f$(theta), and \f$q_l\f$(qlx), \f$\theta_e\f$(thetae), +!! \f$\theta_v\f$(thvx),\f$\theta_{l,v}\f$ (thlvx) including ice water do k=1,km do i=1,im pix(i,k) = psk(i) / prslk(i,k) @@ -403,10 +415,9 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & gotvx(i,k) = g / tvx(i,k) enddo enddo -! -! compute an empirical cloud fraction based on -! Xu & Randall's (1996,JAS) study -! + +!> - Compute an empirical cloud fraction based on +!! Xu and Randall (1996) \cite xu_and_randall_1996 do k = 1, km do i = 1, im plyr(i,k) = 0.01 * prsl(i,k) ! pa to mb (hpa) @@ -433,7 +444,7 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & enddo enddo ! -! compute buoyancy modified by clouds +!> - Compute buoyancy modified by clouds ! do k = 1, km1 do i = 1, im @@ -456,6 +467,8 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! +!> - Initialize diffusion coefficients to 0 and calculate the total +!! radiative heating rate (dku, dkt, radx) do k=1,km1 do i=1,im dku(i,k) = 0. @@ -467,14 +480,31 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & radx(i,k) = tem*(swh(i,k)*xmu(i)+hlw(i,k)) enddo enddo -! +!> - Compute stable/unstable PBL flag (pblflg) based on the total +!! surface energy flux (\e false if the total surface energy flux +!! is into the surface) do i = 1,im sflux(i) = heat(i) + evap(i)*fv*theta(i,1) if(.not.sfcflg(i) .or. sflux(i) <= 0.) pblflg(i)=.false. enddo ! -! compute critical bulk richardson number -! +!> ## Calculate the PBL height +!! The calculation of the boundary layer height follows Troen and Mahrt (1986) \cite troen_and_mahrt_1986 section 3. The approach is to find the level in the column where a modified bulk Richardson number exceeds a critical value. +!! - Compute critical bulk Richardson number (\f$Rb_{cr}\f$) (crb) +!! - For the unstable PBL, crb is a constant (0.25) +!! - For the stable boundary layer (SBL), \f$Rb_{cr}\f$ varies +!! with the surface Rossby number, \f$R_{0}\f$, as given by +!! Vickers and Mahrt (2004) \cite Vickers_2004 +!! \f[ +!! Rb_{cr}=0.16(10^{-7}R_{0})^{-0.18} +!! \f] +!! \f[ +!! R_{0}=\frac{U_{10}}{f_{0}z_{0}} +!! \f] +!! where \f$U_{10}\f$ is the wind speed at 10m above the ground surface, +!! \f$f_0\f$ the Coriolis parameter, and \f$z_{0}\f$ the surface roughness +!! length. To avoid too much variation, we restrict \f$Rb_{cr}\f$ to vary +!! within the range of 0.15~0.35 do i = 1,im if(pblflg(i)) then ! thermal(i) = thvx(i,1) @@ -490,7 +520,7 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & crb(i) = max(min(crb(i), crbmax), crbmin) endif enddo -! +!> - Compute \f$\frac{\Delta t}{\Delta z}\f$ , \f$u_*\f$ do i=1,im dtdz1(i) = dt2 / (zi(i,2)-zi(i,1)) enddo @@ -499,7 +529,8 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & ustar(i) = sqrt(stress(i)) enddo ! -! compute buoyancy (bf) and winshear square +!> - Compute buoyancy \f$\frac{\partial \theta_v}{\partial z}\f$ (bf) +!! and the wind shear squared (shr2) ! do k = 1, km1 do i = 1, im @@ -511,14 +542,18 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & enddo enddo ! -! find pbl height based on bulk richardson number (mrf pbl scheme) +! Find pbl height based on bulk richardson number (mrf pbl scheme) ! and also for diagnostic purpose ! do i=1,im flg(i) = .false. rbup(i) = rbsoil(i) enddo -! +!> - Given the thermal's properties and the critical Richardson number, +!! a loop is executed to find the first level above the surface (kpblx) where +!! the modified Richardson number is greater than the critical Richardson +!! number, using equation 10a from Troen and Mahrt (1996) \cite troen_and_mahrt_1986 +!! (also equation 8 from Hong and Pan (1996) \cite hong_and_pan_1996): do k = 1, kmpbl do i = 1, im if(.not.flg(i)) then @@ -533,6 +568,9 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & endif enddo enddo +!> - Once the level is found, some linear interpolation is performed to find +!! the exact height of the boundary layer top (where \f$R_{i} > Rb_{cr}\f$) +!! and the PBL height (hpbl and kpbl) and the PBL top index are saved. do i = 1,im if(kpblx(i) > 1) then k = kpblx(i) @@ -554,8 +592,15 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & if(kpbl(i) <= 1) pblflg(i)=.false. enddo ! -! compute similarity parameters -! +!> ## Compute Monin-Obukhov similarity parameters +!! - Calculate the Monin-Obukhov nondimensional stability paramter, commonly +!! referred to as \f$\zeta\f$ using the following equation from Businger et al.(1971) \cite businger_et_al_1971 +!! (eqn 28): +!! \f[ +!! \zeta = Ri_{sfc}\frac{F_m^2}{F_h} = \frac{z}{L} +!! \f] +!! where \f$F_m\f$ and \f$F_h\f$ are surface Monin-Obukhov stability functions calculated in sfc_diff.f and +!! \f$L\f$ is the Obukhov length. do i=1,im zol(i) = max(rbsoil(i)*fm(i)*fm(i)/fh(i),rimin) if(sfcflg(i)) then @@ -563,7 +608,17 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & else zol(i) = max(zol(i),zfmin) endif -! +!> - Calculate the nondimensional gradients of momentum and temperature (\f$\phi_m\f$ (phim) and \f$\phi_h\f$(phih)) are calculated using +!! eqns 5 and 6 from Hong and Pan (1996) \cite hong_and_pan_1996 depending on the surface layer stability: +!! - For the unstable and neutral conditions: +!! \f[ +!! \phi_m=(1-16\frac{0.1h}{L})^{-1/4} +!! \phi_h=(1-16\frac{0.1h}{L})^{-1/2} +!! \f] +!! - For the stable regime +!! \f[ +!! \phi_m=\phi_t=(1+5\frac{0.1h}{L}) +!! \f] zol1 = zol(i)*sfcfrac*hpbl(i)/zl(i,1) if(sfcflg(i)) then tem = 1.0 / (1. - aphi16*zol1) @@ -575,6 +630,21 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & endif enddo ! +!> - The \f$z/L\f$ (zol) is used as the stability criterion for the PBL.Currently, +!! strong unstable (convective) PBL for \f$z/L < -0.02\f$ and weakly and moderately +!! unstable PBL for \f$0>z/L>-0.02\f$ +!> - Compute the velocity scale \f$w_s\f$ (wscale) (eqn 22 of Han et al. 2019). It +!! is represented by the value scaled at the top of the surface layer: +!! \f[ +!! w_s=(u_*^3+7\alpha\kappa w_*^3)^{1/3} +!! \f] +!! where \f$u_*\f$ (ustar) is the surface friction velocity,\f$\alpha\f$ is the ratio +!! of the surface layer height to the PBL height (specified as sfcfrac =0.1), +!! \f$\kappa =0.4\f$ is the von Karman constant, and \f$w_*\f$ is the convective velocity +!! scale defined as eqn23 of Han et al.(2019): +!! \f[ +!! w_{*}=[(g/T)\overline{(w'\theta_v^{'})}_0h]^{1/3} +!! \f] do i=1,im if(pblflg(i)) then if(zol(i) < zolcru) then @@ -589,7 +659,8 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & endif enddo ! -! compute a thermal excess +!> ## The counter-gradient terms for temperature and humidity are calculated. +!! - Equation 4 of Hong and Pan (1996) \cite hong_and_pan_1996 and are used to calculate the "scaled virtual temperature excess near the surface" (equation 9 in Hong and Pan (1996) \cite hong_and_pan_1996) for use in the mass-flux algorithm. ! do i = 1,im if(pcnvflg(i)) then @@ -603,7 +674,10 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! look for stratocumulus -! +!> ## Determine whether stratocumulus layers exist and compute quantities +!! - Starting at the PBL top and going downward, if the level is less than 2.5 km +!! and \f$q_l\geq q_{lcr}\f$ then set kcld = k (find the cloud top index in the PBL. +!! If no cloud water above the threshold is hound, \e scuflg is set to F. do i=1,im flg(i) = scuflg(i) enddo @@ -631,7 +705,11 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & do i = 1, im if(scuflg(i) .and. kcld(i)==km1) scuflg(i)=.false. enddo -! +!> - Starting at the PBL top and going downward, if the level is less +!! than the cloud top, find the level of the minimum radiative heating +!! rate wihin the cloud. If the level of the minimum is the lowest model +!! level or the minimum radiative heating rate is positive, then set +!! scuflg to F. do i = 1, im flg(i)=scuflg(i) enddo @@ -655,9 +733,10 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & enddo ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! compute components for mass flux mixing by large thermals +!> ## Compute components for mass flux mixing by large thermals !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! +!> - If the PBL is convective, the updraft properties are initialized +!! to be the same as the state variables. do k = 1, km do i = 1, im if(pcnvflg(i)) then @@ -684,12 +763,14 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & enddo enddo enddo -! +!> - Call mfpbltq(), which is an EDMF parameterization (Siebesma et al.(2007) \cite Siebesma_2007) +!! to take into account nonlocal transport by large eddies. For details of the mfpbltq subroutine, step into its documentation ::mfpbltq call mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,dt2, & pcnvflg,zl,zm,q1,t1,u1,v1,plyr,pix,thlx,thvx, & gdx,hpbl,kpbl,vpert,buou,xmf, & tcko,qcko,ucko,vcko,xlamue,bl_upfr) -! +!> - Call mfscuq(), which is a new mass-flux parameterization for +!! stratocumulus-top-induced turbulence mixing. For details of the mfscuq subroutine, step into its documentation ::mfscuq call mfscuq(im,ix,km,kmscu,ntcw,ntrac1,dt2, & scuflg,zl,zm,q1,t1,u1,v1,plyr,pix, & thlx,thvx,thlvx,gdx,thetae, @@ -697,8 +778,8 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & & tcdo,qcdo,ucdo,vcdo,xlamde,bl_dnfr) ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! compute prandtl number and exchange coefficient varying with height -! + +!> ## Compute Prandtl number \f$P_r\f$ (prn) and exchange coefficient varying with height do k = 1, kmpbl do i = 1, im if(k < kpbl(i)) then @@ -742,8 +823,8 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & ! enddo ! enddo ! -! The background vertical diffusivities in the inversion layers are limited -! to be less than or equal to xkzminv +!> ## The background vertical diffusivities in the inversion layers are limited +!! to be less than or equal to xkzinv ! do k = 1,km1 do i=1,im @@ -758,7 +839,7 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & enddo ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! compute an asymtotic mixing length +!> ## Compute an asymtotic mixing length ! do k = 1, km1 do i = 1, im @@ -818,7 +899,18 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & ! tem = 0.5 * (zi(i,k+1)-zi(i,k)) tem1 = min(tem, rlmnz(i,k)) -! +!> - Following Bougeault and Lacarrere(1989), the characteristic length +!! scale (\f$l_2\f$) (eqn 10 in Han et al.(2019) \cite Han_2019) is given by: +!!\f[ +!! l_2=min(l_{up},l_{down}) +!!\f] +!! and dissipation length scale \f$l_d\f$ is given by: +!!\f[ +!! l_d=(l_{up}l_{down})^{1/2} +!!\f] +!! where \f$l_{up}\f$ and \f$l_{down}\f$ are the distances that a parcel +!! having an initial TKE can travel upward and downward before being stopped +!! by buoyancy effects. ptem2 = min(zlup,zldn) rlam(i,k) = elmfac * ptem2 rlam(i,k) = max(rlam(i,k), tem1) @@ -831,7 +923,8 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & ! enddo enddo -! +!> - Compute the surface layer length scale (\f$l_1\f$) following +!! Nakanishi (2001) \cite Nakanish_2001 (eqn 9 of Han et al.(2019) \cite Han_2019) do k = 1, km1 do i = 1, im tem = vk * zl(i,k) @@ -860,7 +953,7 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & enddo ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! compute eddy diffusivities +!> ## Compute eddy diffusivities !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! do k = 1, km1 @@ -913,8 +1006,8 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & ! enddo enddo -! -! compute a minimum TKE deduced from background diffusivity for momentum. +!> ## Compute TKE. +!! - Compute a minimum TKE deduced from background diffusivity for momentum. ! do k = 1, km1 do i = 1, im @@ -933,7 +1026,7 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & enddo ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! compute buoyancy and shear productions of tke +!> - Compute buoyancy and shear productions of TKE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! do k = 1, km1 @@ -1057,7 +1150,7 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & enddo ! !---------------------------------------------------------------------- -! first predict tke due to tke production & dissipation(diss) +!> - First predict tke due to tke production & dissipation(diss) ! dtn = dt2 / float(ndt) do n = 1, ndt @@ -1075,7 +1168,7 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & enddo enddo ! -! compute updraft & downdraft properties for tke +!> - Compute updraft & downdraft properties for TKE ! do k = 1, km do i = 1, im @@ -1113,7 +1206,7 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & enddo ! !---------------------------------------------------------------------- -! compute tridiagonal matrix elements for turbulent kinetic energy +!> - Compute tridiagonal matrix elements for turbulent kinetic energy ! do i=1,im ad(i,1) = 1.0 @@ -1161,11 +1254,11 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & enddo enddo c -c solve tridiagonal problem for tke +!> - Call tridit() to solve tridiagonal problem for TKE c call tridit(im,km,1,al,ad,au,f1,au,f1) c -c recover tendency of tke +!> - Recover the tendency of tke c do k = 1,km do i = 1,im @@ -1175,7 +1268,7 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & enddo enddo c -c compute tridiagonal matrix elements for heat and moisture +!> ## Compute tridiagonal matrix elements for heat and moisture c do i=1,im ad(i,1) = 1. @@ -1284,11 +1377,11 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & enddo endif c -c solve tridiagonal problem for heat and moisture +!> - Call tridin() to solve tridiagonal problem for heat and moisture c call tridin(im,km,ntrac1,al,ad,au,f1,f2,au,f1,f2) c -c recover tendencies of heat and moisture +!> - Recover the tendencies of heat and moisture c do k = 1,km do i = 1,im @@ -1313,7 +1406,7 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & enddo endif ! -! add tke dissipative heating to temperature tendency +!> ## Add TKE dissipative heating to temperature tendency ! if(dspheat) then do k = 1,km1 @@ -1326,7 +1419,7 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & enddo endif c -c compute tridiagonal matrix elements for momentum +!> ## Compute tridiagonal matrix elements for momentum c do i=1,im ad(i,1) = 1.0 + dtdz1(i) * stress(i) / spd1(i) @@ -1384,11 +1477,11 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & enddo enddo c -c solve tridiagonal problem for momentum +!> - Call tridi2() to solve tridiagonal problem for momentum c call tridi2(im,km,al,ad,au,f1,f2,au,f1,f2) c -c recover tendencies of momentum +!> - Recover the tendencies of momentum c do k = 1,km do i = 1,im @@ -1402,7 +1495,7 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & enddo ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! pbl height for diagnostic purpose +!> ## Save PBL height for diagnostic purpose ! do i = 1, im hpbl(i) = hpblx(i) @@ -1413,5 +1506,5 @@ subroutine satmedmfvdifq_run(ix,im,km,ntrac,ntcw,ntiw,ntke, & return end subroutine satmedmfvdifq_run !> @} - +!! @} end module satmedmfvdifq diff --git a/physics/tridi.f b/physics/tridi.f index 22a35ea9c..bd44bcc86 100644 --- a/physics/tridi.f +++ b/physics/tridi.f @@ -42,6 +42,7 @@ end subroutine tridi1 c----------------------------------------------------------------------- !>\ingroup satmedmf +!>\ingroup satmedmfvdifq !> This subroutine .. subroutine tridi2(l,n,cl,cm,cu,r1,r2,au,a1,a2) cc @@ -84,6 +85,7 @@ end subroutine tridi2 c----------------------------------------------------------------------- !>\ingroup satmedmf +!>\ingroup satmedmfvdifq !> Routine to solve the tridiagonal system to calculate u- and !! v-momentum at \f$ t + \Delta t \f$; part of two-part process to !! calculate time tendencies due to vertical diffusion. @@ -154,6 +156,7 @@ end subroutine tridin c----------------------------------------------------------------------- !>\ingroup satmedmf +!>\ingroup satmedmfvdifq !! This subroutine solves tridiagonal problem for TKE. subroutine tridit(l,n,nt,cl,cm,cu,rt,au,at) !----------------------------------------------------------------------- From ee1065bae62e21319ff55375af1221bebc9dcfca Mon Sep 17 00:00:00 2001 From: "Bin.Liu" Date: Wed, 18 Dec 2019 02:59:06 +0000 Subject: [PATCH 2/6] Connect HAFS version of GFS EDMF PBL scheme with CCPP (Qingfu, Bin, Chunxi, and Weiguo). --- physics/moninedmf_hafs.f | 1555 +++++++++++++++++++++++++++++++++++ physics/moninedmf_hafs.meta | 526 ++++++++++++ 2 files changed, 2081 insertions(+) create mode 100644 physics/moninedmf_hafs.f create mode 100644 physics/moninedmf_hafs.meta diff --git a/physics/moninedmf_hafs.f b/physics/moninedmf_hafs.f new file mode 100644 index 000000000..5c6ff85a8 --- /dev/null +++ b/physics/moninedmf_hafs.f @@ -0,0 +1,1555 @@ +!> \file moninedmf_hafs.f +!! Contains most of the hybrid eddy-diffusivity mass-flux scheme except for the +!! subroutine that calculates the mass flux and updraft properties. + +!> This module contains the CCPP-compliant hybrid eddy-diffusivity mass-flux +!! scheme. + module hedmf_hafs + + contains + +!> \section arg_table_hedmf_hafs_init Argument Table +!! \htmlinclude hedmf_hafs_init.html +!! + subroutine hedmf_hafs_init (moninq_fac,errmsg,errflg) + use machine, only : kind_phys + implicit none + real(kind=kind_phys), intent(in ) :: moninq_fac + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + ! Initialize CCPP error handling variables + errmsg = '' + errflg = 0 + + if (moninq_fac == 0) then + errflg = 1 + write(errmsg,'(*(a))') 'Logic error: moninq_fac == 0', & + & ' is incompatible with moninedmf_hafs' + end if + end subroutine hedmf_hafs_init + + subroutine hedmf_hafs_finalize () + end subroutine hedmf_hafs_finalize + + +!> \defgroup HEDMF GFS Hybrid Eddy-Diffusivity Mass-Flux (HEDMF) Scheme Module +!! @{ +!! \brief This subroutine contains all of logic for the +!! Hybrid EDMF PBL scheme except for the calculation of +!! the updraft properties and mass flux. +!! +!> \section arg_table_hedmf_hafs_run Argument Table +!! \htmlinclude hedmf_hafs_run.html +!! +!! \section general_edmf GFS Hybrid EDMF General Algorithm +!! -# Compute preliminary variables from input arguments. +!! -# Calculate the first estimate of the PBL height ("Predictor step"). +!! -# Calculate Monin-Obukhov similarity parameters. +!! -# Update thermal properties of surface parcel and recompute PBL height ("Corrector step"). +!! -# Determine whether stratocumulus layers exist and compute quantities needed for enhanced diffusion. +!! -# Calculate the inverse Prandtl number. +!! -# Compute diffusion coefficients below the PBL top. +!! -# Compute diffusion coefficients above the PBL top. +!! -# If the PBL is convective, call the mass flux scheme to replace the countergradient terms. +!! -# Compute enhanced diffusion coefficients related to stratocumulus-topped PBLs. +!! -# Solve for the temperature and moisture tendencies due to vertical mixing. +!! -# Calculate heating due to TKE dissipation and add to the tendency for temperature. +!! -# Solve for the horizontal momentum tendencies and add them to output tendency terms. +!! \section detailed_hedmf GFS Hybrid HEDMF Detailed Algorithm +!! @{ + subroutine hedmf_hafs_run(ix,im,km,ntrac,ntcw,dv,du,tau,rtg, & + & u1,v1,t1,q1,swh,hlw,xmu, & + & psk,rbsoil,zorl,u10m,v10m,fm,fh, & + & tsea,heat,evap,stress,spd1,kpbl, & + & prsi,del,prsl,prslk,phii,phil,delt,dspheat, & + & dusfc,dvsfc,dtsfc,dqsfc,hpbl,hgamt,hgamq,dkt, & + & kinver,xkzm_m,xkzm_h,xkzm_s,lprnt,ipr, & + & xkzminv,moninq_fac,islimsk,errmsg,errflg) +! + use machine , only : kind_phys + use funcphys , only : fpvs + use physcons, grav => con_g, rd => con_rd, cp => con_cp & + &, hvap => con_hvap, fv => con_fvirt + implicit none +! +! arguments +! + logical, intent(in) :: lprnt + integer, intent(in) :: ipr + integer, intent(in) :: ix, im, km, ntrac, ntcw, kinver(im) + integer, intent(in) :: islimsk(1:im) + integer, intent(out) :: kpbl(im) + +! + real(kind=kind_phys), intent(in) :: delt, xkzm_m, xkzm_h, xkzm_s + real(kind=kind_phys), intent(in) :: xkzminv, moninq_fac + real(kind=kind_phys), intent(inout) :: dv(im,km), du(im,km), & + & tau(im,km), rtg(im,km,ntrac) + real(kind=kind_phys), intent(in) :: & + & u1(ix,km), v1(ix,km), & + & t1(ix,km), q1(ix,km,ntrac), & + & swh(ix,km), hlw(ix,km), & + & xmu(im), psk(im), & + & rbsoil(im), zorl(im), & + & u10m(im), v10m(im), & + & fm(im), fh(im), & + & tsea(im), & + & heat(im), evap(im), & + & stress(im), spd1(im) + real(kind=kind_phys), intent(in) :: & + & prsi(ix,km+1), del(ix,km), & + & prsl(ix,km), prslk(ix,km), & + & phii(ix,km+1), phil(ix,km) + real(kind=kind_phys), intent(out) :: & + & dusfc(im), dvsfc(im), & + & dtsfc(im), dqsfc(im), & + & hpbl(im), dkt(im,km-1) + + real(kind=kind_phys), intent(inout) :: & + & hgamt(im), hgamq(im) +! + logical, intent(in) :: dspheat +! flag for tke dissipative heating + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + +! +! locals +! + integer i,iprt,is,iun,k,kk,km1,kmpbl,latd,lond + integer lcld(im),icld(im),kcld(im),krad(im) + integer kx1(im), kpblx(im) +! +! real(kind=kind_phys) betaq(im), betat(im), betaw(im), + real(kind=kind_phys) phih(im), phim(im), hpblx(im), & + & rbdn(im), rbup(im), & + & beta(im), sflux(im), & + & z0(im), crb(im), wstar(im), & + & zol(im), ustmin(im), ustar(im), & + & thermal(im),wscale(im), wscaleu(im) +! + real(kind=kind_phys) theta(im,km),thvx(im,km), thlvx(im,km), & + & qlx(im,km), thetae(im,km), & + & qtx(im,km), bf(im,km-1), diss(im,km), & + & radx(im,km-1), & + & govrth(im), hrad(im), & +! & hradm(im), radmin(im), vrad(im), & + & radmin(im), vrad(im), & + & zd(im), zdd(im), thlvx1(im) +! + real(kind=kind_phys) rdzt(im,km-1),dktx(im,km-1), & + & zi(im,km+1), zl(im,km), xkzo(im,km-1), & + & dku(im,km-1), xkzmo(im,km-1), & + & cku(im,km-1), ckt(im,km-1), & + & ti(im,km-1), shr2(im,km-1), & + & al(im,km-1), ad(im,km), & + & au(im,km-1), a1(im,km), & + & a2(im,km*ntrac) +! + real(kind=kind_phys) tcko(im,km), qcko(im,km,ntrac), & + & ucko(im,km), vcko(im,km), xmf(im,km) +! + real(kind=kind_phys) prinv(im), rent(im) +! + logical pblflg(im), sfcflg(im), scuflg(im), flg(im) + logical ublflg(im), pcnvflg(im) +! +! pcnvflg: true for convective(strongly unstable) pbl +! ublflg: true for unstable but not convective(strongly unstable) pbl +! + real(kind=kind_phys) aphi16, aphi5, bvf2, wfac, + & cfac, conq, cont, conw, + & dk, dkmax, dkmin, + & dq1, dsdz2, dsdzq, dsdzt, + & dsdzu, dsdzv, + & dsig, dt2, dthe1, dtodsd, + & dtodsu, dw2, dw2min, g, + & gamcrq, gamcrt, gocp, + & gravi, f0, + & prnum, prmax, prmin, pfac, crbcon, + & qmin, tdzmin, qtend, crbmin,crbmax, + & rbint, rdt, rdz, qlmin, + & ri, rimin, rl2, rlam, rlamun, + & rone, rzero, sfcfrac, + & spdk2, sri, zol1, zolcr, zolcru, + & robn, ttend, + & utend, vk, vk2, + & ust3, wst3, + & vtend, zfac, vpert, cteit, + & rentf1, rentf2, radfac, + & zfmin, zk, tem, tem1, tem2, + & xkzm, xkzmu, + & ptem, ptem1, ptem2, tx1(im), tx2(im) +! + real(kind=kind_phys) zstblmax,h1, h2, qlcr, actei, + & cldtime + +!! for aplha + real(kind=kind_phys) WSPM(IM,KM-1) + integer kLOC ! RGF + real :: xDKU, ALPHA ! RGF + + integer :: useshape + real :: smax,ashape,sz2h, sksfc,skmax,ashape1,skminusk0, hmax + + +!cc + parameter(gravi=1.0/grav) + parameter(g=grav) + parameter(gocp=g/cp) + parameter(cont=cp/g,conq=hvap/g,conw=1.0/g) ! for del in pa +! parameter(cont=1000.*cp/g,conq=1000.*hvap/g,conw=1000./g) ! for del in kpa + parameter(rlam=30.0,vk=0.4,vk2=vk*vk) + parameter(prmin=0.25,prmax=4.,zolcr=0.2,zolcru=-0.5) + parameter(dw2min=0.0001,dkmin=0.0,dkmax=1000.,rimin=-100.) + parameter(crbcon=0.25,crbmin=0.15,crbmax=0.35) + parameter(wfac=7.0,cfac=6.5,pfac=2.0,sfcfrac=0.1) +! parameter(qmin=1.e-8,xkzm=1.0,zfmin=1.e-8,aphi5=5.,aphi16=16.) + parameter(qmin=1.e-8, zfmin=1.e-8,aphi5=5.,aphi16=16.) + parameter(tdzmin=1.e-3,qlmin=1.e-12,f0=1.e-4) + parameter(h1=0.33333333,h2=0.66666667) +! parameter(cldtime=500.,xkzminv=0.3) + parameter(cldtime=500.) +! parameter(cldtime=500.,xkzmu=3.0,xkzminv=0.3) +! parameter(gamcrt=3.,gamcrq=2.e-3,rlamun=150.0) + parameter(gamcrt=3.,gamcrq=0.,rlamun=150.0) + parameter(rentf1=0.2,rentf2=1.0,radfac=0.85) + parameter(iun=84) +! +! parameter (zstblmax = 2500., qlcr=1.0e-5) +! parameter (zstblmax = 2500., qlcr=3.0e-5) +! parameter (zstblmax = 2500., qlcr=3.5e-5) +! parameter (zstblmax = 2500., qlcr=1.0e-4) + parameter (zstblmax = 2500., qlcr=3.5e-5) +! parameter (actei = 0.23) + parameter (actei = 0.7) + +! HAFS PBL: height-dependent ALPHA + useshape=2 !0-- no change, origincal ALPHA adjustment,1-- shape1, 2-- shape2(adjust above sfc) + alpha=moninq_fac + + ! write(0,*)'in PBL,alpha=',alpha + + ! write(0,*)'islimsk=',(islimsk(i),i=1,im) + +c +c----------------------------------------------------------------------- +c + 601 format(1x,' moninp lat lon step hour ',3i6,f6.1) + 602 format(1x,' k',' z',' t',' th', + 1 ' tvh',' q',' u',' v', + 2 ' sp') + 603 format(1x,i5,8f9.1) + 604 format(1x,' sfc',9x,f9.1,18x,f9.1) + 605 format(1x,' k zl spd2 thekv the1v' + 1 ,' thermal rbup') + 606 format(1x,i5,6f8.2) + 607 format(1x,' kpbl hpbl fm fh hgamt', + 1 ' hgamq ws ustar cd ch') + 608 format(1x,i5,9f8.2) + 609 format(1x,' k pr dkt dku ',i5,3f8.2) + 610 format(1x,' k pr dkt dku ',i5,3f8.2,' l2 ri t2', + 1 ' sr2 ',2f8.2,2e10.2) +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +! Initialize CCPP error handling variables + errmsg = '' + errflg = 0 +!> ## Compute preliminary variables from input arguments + +! compute preliminary variables +! + if (ix .lt. im) stop +! +! iprt = 0 +! if(iprt.eq.1) then +!cc latd = 0 +! lond = 0 +! else +!cc latd = 0 +! lond = 0 +! endif +! + dt2 = delt + rdt = 1. / dt2 + km1 = km - 1 + kmpbl = km / 2 +!> - Compute physical height of the layer centers and interfaces from the geopotential height (zi and zl) + do k=1,km + do i=1,im + zi(i,k) = phii(i,k) * gravi + zl(i,k) = phil(i,k) * gravi + enddo + enddo + do i=1,im + zi(i,km+1) = phii(i,km+1) * gravi + enddo +!> - Compute reciprocal of \f$ \Delta z \f$ (rdzt) + do k = 1,km1 + do i=1,im + rdzt(i,k) = 1.0 / (zl(i,k+1) - zl(i,k)) + enddo + enddo +!> - Compute reciprocal of pressure (tx1, tx2) + do i=1,im + kx1(i) = 1 + tx1(i) = 1.0 / prsi(i,1) + tx2(i) = tx1(i) + enddo +!> - Compute background vertical diffusivities for scalars and momentum (xkzo and xkzmo) + do k = 1,km1 + do i=1,im + xkzo(i,k) = 0.0 + xkzmo(i,k) = 0.0 + if (k < kinver(i)) then +! vertical background diffusivity + ptem = prsi(i,k+1) * tx1(i) + tem1 = 1.0 - ptem + tem1 = tem1 * tem1 * 10.0 + xkzo(i,k) = xkzm_h * min(1.0, exp(-tem1)) + +! vertical background diffusivity for momentum + if (ptem >= xkzm_s) then + xkzmo(i,k) = xkzm_m + kx1(i) = k + 1 + else + if (k == kx1(i) .and. k > 1) tx2(i) = 1.0 / prsi(i,k) + tem1 = 1.0 - prsi(i,k+1) * tx2(i) + tem1 = tem1 * tem1 * 5.0 + xkzmo(i,k) = xkzm_m * min(1.0, exp(-tem1)) + endif + endif + enddo + enddo + +! if (lprnt) then +! print *,' xkzo=',(xkzo(ipr,k),k=1,km1) +! print *,' xkzmo=',(xkzmo(ipr,k),k=1,km1) +! endif +! +! diffusivity in the inversion layer is set to be xkzminv (m^2/s) +!> - The background scalar vertical diffusivity is limited to be less than or equal to xkzminv + do k = 1,kmpbl + do i=1,im +! if(zi(i,k+1) > 200..and.zi(i,k+1) < zstblmax) then + if(zi(i,k+1) > 250.) then + tem1 = (t1(i,k+1)-t1(i,k)) * rdzt(i,k) + if(tem1 > 1.e-5) then + xkzo(i,k) = min(xkzo(i,k),xkzminv) + endif + endif + enddo + enddo +!> - Some output variables and logical flags are initialized + do i = 1,im + z0(i) = 0.01 * zorl(i) + dusfc(i) = 0. + dvsfc(i) = 0. + dtsfc(i) = 0. + dqsfc(i) = 0. + wscale(i)= 0. + wscaleu(i)= 0. + kpbl(i) = 1 + hpbl(i) = zi(i,1) + hpblx(i) = zi(i,1) + pblflg(i)= .true. + sfcflg(i)= .true. + if(rbsoil(i) > 0.) sfcflg(i) = .false. + ublflg(i)= .false. + pcnvflg(i)= .false. + scuflg(i)= .true. + if(scuflg(i)) then + radmin(i)= 0. + rent(i) = rentf1 + hrad(i) = zi(i,1) +! hradm(i) = zi(i,1) + krad(i) = 1 + icld(i) = 0 + lcld(i) = km1 + kcld(i) = km1 + zd(i) = 0. + endif + enddo +!> - Compute \f$\theta\f$ (theta), \f$q_l\f$ (qlx), \f$q_t\f$ (qtx), \f$\theta_e\f$ (thetae), \f$\theta_v\f$ (thvx), \f$\theta_{l,v}\f$ (thlvx) + do k = 1,km + do i = 1,im + theta(i,k) = t1(i,k) * psk(i) / prslk(i,k) + qlx(i,k) = max(q1(i,k,ntcw),qlmin) + qtx(i,k) = max(q1(i,k,1),qmin)+qlx(i,k) + ptem = qlx(i,k) + ptem1 = hvap*max(q1(i,k,1),qmin)/(cp*t1(i,k)) + thetae(i,k)= theta(i,k)*(1.+ptem1) + thvx(i,k) = theta(i,k)*(1.+fv*max(q1(i,k,1),qmin)-ptem) + ptem2 = theta(i,k)-(hvap/cp)*ptem + thlvx(i,k) = ptem2*(1.+fv*qtx(i,k)) + enddo + enddo +!> - Initialize diffusion coefficients to 0 and calculate the total radiative heating rate (dku, dkt, radx) + do k = 1,km1 + do i = 1,im + dku(i,k) = 0. + dkt(i,k) = 0. + dktx(i,k) = 0. + cku(i,k) = 0. + ckt(i,k) = 0. + tem = zi(i,k+1)-zi(i,k) + radx(i,k) = tem*(swh(i,k)*xmu(i)+hlw(i,k)) + enddo + enddo +!> - Set lcld to first index above 2.5km + do i=1,im + flg(i) = scuflg(i) + enddo + do k = 1, km1 + do i=1,im + if(flg(i).and.zl(i,k) >= zstblmax) then + lcld(i)=k + flg(i)=.false. + endif + enddo + enddo +! +! compute virtual potential temp gradient (bf) and winshear square +!> - Compute \f$\frac{\partial \theta_v}{\partial z}\f$ (bf) and the wind shear squared (shr2) + do k = 1, km1 + do i = 1, im + rdz = rdzt(i,k) + bf(i,k) = (thvx(i,k+1)-thvx(i,k))*rdz + ti(i,k) = 2./(t1(i,k)+t1(i,k+1)) + dw2 = (u1(i,k)-u1(i,k+1))**2 + & + (v1(i,k)-v1(i,k+1))**2 + shr2(i,k) = max(dw2,dw2min)*rdz*rdz + enddo + enddo +!> - Calculate \f$\frac{g}{\theta}\f$ (govrth), \f$\beta = \frac{\Delta t}{\Delta z}\f$ (beta), \f$u_*\f$ (ustar), total surface flux (sflux), and set pblflag to false if the total surface energy flux is into the surface + do i = 1,im + govrth(i) = g/theta(i,1) + enddo +! + do i=1,im + beta(i) = dt2 / (zi(i,2)-zi(i,1)) + enddo +! + do i=1,im + ustar(i) = sqrt(stress(i)) + enddo +! + do i = 1,im + sflux(i) = heat(i) + evap(i)*fv*theta(i,1) + if(.not.sfcflg(i) .or. sflux(i) <= 0.) pblflg(i)=.false. + enddo +!> ## Calculate the first estimate of the PBL height (``Predictor step") +!! The calculation of the boundary layer height follows Troen and Mahrt (1986) \cite troen_and_mahrt_1986 section 3. The approach is to find the level in the column where a modified bulk Richardson number exceeds a critical value. +!! +!! The temperature of the thermal is of primary importance. For the initial estimate of the PBL height, the thermal is assumed to have one of two temperatures. If the boundary layer is stable, the thermal is assumed to have a temperature equal to the surface virtual temperature. Otherwise, the thermal is assumed to have the same virtual potential temperature as the lowest model level. For the stable case, the critical bulk Richardson number becomes a function of the wind speed and roughness length, otherwise it is set to a tunable constant. +! compute the pbl height +! + do i=1,im + flg(i) = .false. + rbup(i) = rbsoil(i) + + IF ( ALPHA .GT. 0.0) THEN ! ALPHA + + if(pblflg(i)) then + thermal(i) = thvx(i,1) + crb(i) = crbcon + else + thermal(i) = tsea(i)*(1.+fv*max(q1(i,1,1),qmin)) + tem = sqrt(u10m(i)**2+v10m(i)**2) + tem = max(tem, 1.) + robn = tem / (f0 * z0(i)) + tem1 = 1.e-7 * robn + crb(i) = 0.16 * (tem1 ** (-0.18)) + crb(i) = max(min(crb(i), crbmax), crbmin) + endif + + ELSE +! use variable Ri for all conditions + if(pblflg(i)) then + thermal(i) = thvx(i,1) + else + thermal(i) = tsea(i)*(1.+fv*max(q1(i,1,1),qmin)) + endif + tem = sqrt(u10m(i)**2+v10m(i)**2) + tem = max(tem, 1.) + robn = tem / (f0 * z0(i)) + tem1 = 1.e-7 * robn +! crb(i) = 0.16 * (tem1 ** (-0.18)) + crb(i) = crbcon + IF(islimsk(i).ne.0) crb(I) = 0.16*(tem1)**(-0.18) + IF(islimsk(i).eq.0) crb(I) = 0.25*(tem1)**(-0.18) + crb(i) = max(min(crb(i), crbmax), crbmin) + ENDIF ! ALPHA + + enddo + +!> Given the thermal's properties and the critical Richardson number, a loop is executed to find the first level above the surface where the modified Richardson number is greater than the critical Richardson number, using equation 10a from Troen and Mahrt (1986) \cite troen_and_mahrt_1986 (also equation 8 from Hong and Pan (1996) \cite hong_and_pan_1996): +!! \f[ +!! h = Ri\frac{T_0\left|\vec{v}(h)\right|^2}{g\left(\theta_v(h) - \theta_s\right)} +!! \f] +!! where \f$h\f$ is the PBL height, \f$Ri\f$ is the Richardson number, \f$T_0\f$ is the virtual potential temperature near the surface, \f$\left|\vec{v}\right|\f$ is the wind speed, and \f$\theta_s\f$ is for the thermal. Rearranging this equation to calculate the modified Richardson number at each level, k, for comparison with the critical value yields: +!! \f[ +!! Ri_k = gz(k)\frac{\left(\theta_v(k) - \theta_s\right)}{\theta_v(1)*\vec{v}(k)} +!! \f] + do k = 1, kmpbl + do i = 1, im + if(.not.flg(i)) then + rbdn(i) = rbup(i) + spdk2 = max((u1(i,k)**2+v1(i,k)**2),1.) + rbup(i) = (thvx(i,k)-thermal(i))* + & (g*zl(i,k)/thvx(i,1))/spdk2 + kpbl(i) = k + flg(i) = rbup(i) > crb(i) + endif + enddo + enddo + +!> Once the level is found, some linear interpolation is performed to find the exact height of the boundary layer top (where \f$Ri = Ri_{cr}\f$) and the PBL height and the PBL top index are saved (hpblx and kpblx, respectively) + do i = 1,im + if(kpbl(i) > 1) then + k = kpbl(i) + if(rbdn(i) >= crb(i)) then + rbint = 0. + elseif(rbup(i) <= crb(i)) then + rbint = 1. + else + rbint = (crb(i)-rbdn(i))/(rbup(i)-rbdn(i)) + endif + hpbl(i) = zl(i,k-1) + rbint*(zl(i,k)-zl(i,k-1)) + if(hpbl(i) < zi(i,kpbl(i))) kpbl(i) = kpbl(i) - 1 + else + hpbl(i) = zl(i,1) + kpbl(i) = 1 + endif + kpblx(i) = kpbl(i) + hpblx(i) = hpbl(i) + enddo +! +! compute similarity parameters +!> ## Calculate Monin-Obukhov similarity parameters +!! Using the initial guess for the PBL height, Monin-Obukhov similarity parameters are calculated. They are needed to refine the PBL height calculation and for calculating diffusion coefficients. +!! +!! First, calculate the Monin-Obukhov nondimensional stability parameter, commonly referred to as \f$\zeta\f$ using the following equation from Businger et al. (1971) \cite businger_et_al_1971 (equation 28): +!! \f[ +!! \zeta = Ri_{sfc}\frac{F_m^2}{F_h} = \frac{z}{L} +!! \f] +!! where \f$F_m\f$ and \f$F_h\f$ are surface Monin-Obukhov stability functions calculated in sfc_diff.f and \f$L\f$ is the Obukhov length. Then, the nondimensional gradients of momentum and temperature (phim and phih) are calculated using equations 5 and 6 from Hong and Pan (1996) \cite hong_and_pan_1996 depending on the surface layer stability. Then, the velocity scale valid for the surface layer (\f$w_s\f$, wscale) is calculated using equation 3 from Hong and Pan (1996) \cite hong_and_pan_1996. For the neutral and unstable PBL above the surface layer, the convective velocity scale, \f$w_*\f$, is calculated according to: +!! \f[ +!! w_* = \left(\frac{g}{\theta_0}h\overline{w'\theta_0'}\right)^{1/3} +!! \f] +!! and the mixed layer velocity scale is then calculated with equation 6 from Troen and Mahrt (1986) \cite troen_and_mahrt_1986 +!! \f[ +!! w_s = (u_*^3 + 7\epsilon k w_*^3)^{1/3} +!! \f] + do i=1,im + zol(i) = max(rbsoil(i)*fm(i)*fm(i)/fh(i),rimin) + if(sfcflg(i)) then + zol(i) = min(zol(i),-zfmin) + else + zol(i) = max(zol(i),zfmin) + endif + zol1 = zol(i)*sfcfrac*hpbl(i)/zl(i,1) + if(sfcflg(i)) then +! phim(i) = (1.-aphi16*zol1)**(-1./4.) +! phih(i) = (1.-aphi16*zol1)**(-1./2.) + tem = 1.0 / (1. - aphi16*zol1) + phih(i) = sqrt(tem) + phim(i) = sqrt(phih(i)) + else + phim(i) = 1. + aphi5*zol1 + phih(i) = phim(i) + endif + wscale(i) = ustar(i)/phim(i) + ustmin(i) = ustar(i)/aphi5 + wscale(i) = max(wscale(i),ustmin(i)) + enddo + do i=1,im + if(pblflg(i)) then + if(zol(i) < zolcru .and. kpbl(i) > 1) then + pcnvflg(i) = .true. + else + ublflg(i) = .true. + endif + wst3 = govrth(i)*sflux(i)*hpbl(i) + wstar(i)= wst3**h1 + ust3 = ustar(i)**3. + wscaleu(i) = (ust3+wfac*vk*wst3*sfcfrac)**h1 + wscaleu(i) = max(wscaleu(i),ustmin(i)) + endif + enddo +! +! compute counter-gradient mixing term for heat and moisture +!> ## Update thermal properties of surface parcel and recompute PBL height ("Corrector step"). +!! Next, the counter-gradient terms for temperature and humidity are calculated using equation 4 of Hong and Pan (1996) \cite hong_and_pan_1996 and are used to calculate the "scaled virtual temperature excess near the surface" (equation 9 in Hong and Pan (1996) \cite hong_and_pan_1996) so that the properties of the thermal are updated to recalculate the PBL height. + do i = 1,im + if(ublflg(i)) then + hgamt(i) = min(cfac*heat(i)/wscaleu(i),gamcrt) + hgamq(i) = min(cfac*evap(i)/wscaleu(i),gamcrq) + vpert = hgamt(i) + hgamq(i)*fv*theta(i,1) + vpert = min(vpert,gamcrt) + thermal(i)= thermal(i)+max(vpert,0.) + hgamt(i) = max(hgamt(i),0.0) + hgamq(i) = max(hgamq(i),0.0) + endif + enddo +! +! enhance the pbl height by considering the thermal excess +!> The PBL height calculation follows the same procedure as the predictor step, except that it uses an updated virtual potential temperature for the thermal. + do i=1,im + flg(i) = .true. + if(ublflg(i)) then + flg(i) = .false. + rbup(i) = rbsoil(i) + endif + enddo + do k = 2, kmpbl + do i = 1, im + if(.not.flg(i)) then + rbdn(i) = rbup(i) + spdk2 = max((u1(i,k)**2+v1(i,k)**2),1.) + rbup(i) = (thvx(i,k)-thermal(i))* + & (g*zl(i,k)/thvx(i,1))/spdk2 + kpbl(i) = k + flg(i) = rbup(i) > crb(i) + endif + enddo + enddo + do i = 1,im + if(ublflg(i)) then + k = kpbl(i) + if(rbdn(i) >= crb(i)) then + rbint = 0. + elseif(rbup(i) <= crb(i)) then + rbint = 1. + else + rbint = (crb(i)-rbdn(i))/(rbup(i)-rbdn(i)) + endif + hpbl(i) = zl(i,k-1) + rbint*(zl(i,k)-zl(i,k-1)) + if(hpbl(i) < zi(i,kpbl(i))) kpbl(i) = kpbl(i) - 1 + if(kpbl(i) <= 1) then + ublflg(i) = .false. + pblflg(i) = .false. + endif + endif + enddo +! +! look for stratocumulus +!> ## Determine whether stratocumulus layers exist and compute quantities needed for enhanced diffusion +!! - Starting at the PBL top and going downward, if the level is less than 2.5 km and \f$q_l>q_{l,cr}\f$ then set kcld = k (find the cloud top index in the PBL). If no cloud water above the threshold is found, scuflg is set to F. + do i = 1, im + flg(i)=scuflg(i) + enddo + do k = kmpbl,1,-1 + do i = 1, im + if(flg(i) .and. k <= lcld(i)) then + if(qlx(i,k).ge.qlcr) then + kcld(i)=k + flg(i)=.false. + endif + endif + enddo + enddo + do i = 1, im + if(scuflg(i) .and. kcld(i)==km1) scuflg(i)=.false. + enddo +!> - Starting at the PBL top and going downward, if the level is less than the cloud top, find the level of the minimum radiative heating rate within the cloud. If the level of the minimum is the lowest model level or the minimum radiative heating rate is positive, then set scuflg to F. + do i = 1, im + flg(i)=scuflg(i) + enddo + do k = kmpbl,1,-1 + do i = 1, im + if(flg(i) .and. k <= kcld(i)) then + if(qlx(i,k) >= qlcr) then + if(radx(i,k) < radmin(i)) then + radmin(i)=radx(i,k) + krad(i)=k + endif + else + flg(i)=.false. + endif + endif + enddo + enddo + do i = 1, im + if(scuflg(i) .and. krad(i) <= 1) scuflg(i)=.false. + if(scuflg(i) .and. radmin(i)>=0.) scuflg(i)=.false. + enddo +!> - Starting at the PBL top and going downward, count the number of levels below the minimum radiative heating rate level that have cloud water above the threshold. If there are none, then set the scuflg to F. + do i = 1, im + flg(i)=scuflg(i) + enddo + do k = kmpbl,2,-1 + do i = 1, im + if(flg(i) .and. k <= krad(i)) then + if(qlx(i,k) >= qlcr) then + icld(i)=icld(i)+1 + else + flg(i)=.false. + endif + endif + enddo + enddo + do i = 1, im + if(scuflg(i) .and. icld(i) < 1) scuflg(i)=.false. + enddo +!> - Find the height of the interface where the minimum in radiative heating rate is located. If this height is less than the second model interface height, then set the scuflg to F. + do i = 1, im + if(scuflg(i)) then + hrad(i) = zi(i,krad(i)+1) +! hradm(i)= zl(i,krad(i)) + endif + enddo +! + do i = 1, im + if(scuflg(i) .and. hrad(i) - Calculate the hypothetical \f$\theta_v\f$ at the minimum radiative heating level that a parcel would reach due to radiative cooling after a typical cloud turnover time spent at that level. + do i = 1, im + if(scuflg(i)) then + k = krad(i) + tem = zi(i,k+1)-zi(i,k) + tem1 = cldtime*radmin(i)/tem + thlvx1(i) = thlvx(i,k)+tem1 +! if(thlvx1(i) > thlvx(i,k-1)) scuflg(i)=.false. + endif + enddo +!> - Determine the distance that a parcel would sink downwards starting from the level of minimum radiative heating rate by comparing the hypothetical minimum \f$\theta_v\f$ calculated above with the environmental \f$\theta_v\f$. + do i = 1, im + flg(i)=scuflg(i) + enddo + do k = kmpbl,1,-1 + do i = 1, im + if(flg(i) .and. k <= krad(i))then + if(thlvx1(i) <= thlvx(i,k))then + tem=zi(i,k+1)-zi(i,k) + zd(i)=zd(i)+tem + else + flg(i)=.false. + endif + endif + enddo + enddo +!> - Calculate the cloud thickness, where the cloud top is the in-cloud minimum radiative heating level and the bottom is determined previously. + do i = 1, im + if(scuflg(i))then + kk = max(1, krad(i)+1-icld(i)) + zdd(i) = hrad(i)-zi(i,kk) + endif + enddo +!> - Find the largest between the cloud thickness and the distance of a sinking parcel, then determine the smallest of that number and the height of the minimum in radiative heating rate. Set this number to \f$zd\f$. Using \f$zd\f$, calculate the characteristic velocity scale of cloud-top radiative cooling-driven turbulence. + do i = 1, im + if(scuflg(i))then + zd(i) = max(zd(i),zdd(i)) + zd(i) = min(zd(i),hrad(i)) + tem = govrth(i)*zd(i)*(-radmin(i)) + vrad(i)= tem**h1 + endif + enddo +! +! compute inverse prandtl number +!> ## Calculate the inverse Prandtl number +!! For an unstable PBL, the Prandtl number is calculated according to Hong and Pan (1996) \cite hong_and_pan_1996, equation 10, whereas for a stable boundary layer, the Prandtl number is simply \f$Pr = \frac{\phi_h}{\phi_m}\f$. + do i = 1, im + if(ublflg(i)) then + tem = phih(i)/phim(i)+cfac*vk*sfcfrac + else + tem = phih(i)/phim(i) + endif + prinv(i) = 1.0 / tem + prinv(i) = min(prinv(i),prmax) + prinv(i) = max(prinv(i),prmin) + enddo + do i = 1, im + if(zol(i) > zolcr) then + kpbl(i) = 1 + endif + enddo + +!!! HAFS PBL, Bgin adjustment +! RGF determine wspd at roughly 500 m above surface, or as close as possible, +! reuse SPDK2 +! zi(i,k) is AGL, right? May not matter if applied only to water grid points + if(moninq_fac.lt.0)then + + DO I=1,IM + SPDK2 = 0. + WSPM(i,1) = 0. + DO K = 1, KMPBL ! kmpbl is like a max possible pbl height + if(zi(i,k).le.500.and.zi(i,k+1).gt.500.)then ! find level bracketing 500 m + SPDK2 = SQRT(U1(i,k)*U1(i,k)+V1(i,k)*V1(i,k)) ! wspd near 500 m + WSPM(i,1) = SPDK2/0.6 ! now the Km limit for 500 m. just store in K=1 + WSPM(i,2) = float(k) ! height of level at gridpoint i. store in K=2 +! if(i.eq.25) print *,' IK ',i,k,' ZI ',zi(i,k), ' WSPM1 ',wspm(i,1),' +! KMPBL ',kmpbl,' KPBL ',kpbl(i) + endif + ENDDO + ENDDO ! i + + endif ! moninq_fac < 0 + + +! +! compute diffusion coefficients below pbl +!> ## Compute diffusion coefficients below the PBL top +!! Below the PBL top, the diffusion coefficients (\f$K_m\f$ and \f$K_h\f$) are calculated according to equation 2 in Hong and Pan (1996) \cite hong_and_pan_1996 where a different value for \f$w_s\f$ (PBL vertical velocity scale) is used depending on the PBL stability. \f$K_h\f$ is calculated from \f$K_m\f$ using the Prandtl number. The calculated diffusion coefficients are checked so that they are bounded by maximum values and the local background diffusion coefficients. + + IF (ALPHA > 0) THEN ! AAAAAAAAAAAAAAAAAAAAAAAAAAA + + do k = 1, kmpbl + do i=1,im + if(k < kpbl(i)) then +! zfac = max((1.-(zi(i,k+1)-zl(i,1))/ +! 1 (hpbl(i)-zl(i,1))), zfmin) + zfac = max((1.-zi(i,k+1)/hpbl(i)), zfmin) + tem = zi(i,k+1) * (zfac**pfac) * moninq_fac ! lmh suggested by kg + if(pblflg(i)) then + tem1 = vk * wscaleu(i) * tem +! dku(i,k) = xkzmo(i,k) + tem1 +! dkt(i,k) = xkzo(i,k) + tem1 * prinv(i) + dku(i,k) = tem1 + dkt(i,k) = tem1 * prinv(i) + else + tem1 = vk * wscale(i) * tem +! dku(i,k) = xkzmo(i,k) + tem1 +! dkt(i,k) = xkzo(i,k) + tem1 * prinv(i) + dku(i,k) = tem1 + dkt(i,k) = tem1 * prinv(i) + endif + dku(i,k) = min(dku(i,k),dkmax) + dku(i,k) = max(dku(i,k),xkzmo(i,k)) + dkt(i,k) = min(dkt(i,k),dkmax) + dkt(i,k) = max(dkt(i,k),xkzo(i,k)) + dktx(i,k)= dkt(i,k) + endif + enddo + enddo + + ELSE ! ALPHA <0 AAAAAAAAAAAAA + + do i=1,im + do k = 1, kmpbl + if(k < kpbl(i)) then +! zfac = max((1.-(zi(i,k+1)-zl(i,1))/ +! 1 (hpbl(i)-zl(i,1))), zfmin) + zfac = max((1.-zi(i,k+1)/hpbl(i)), zfmin) + ! tem = zi(i,k+1) * (zfac**pfac) * moninq_fac ! lmh suggested by kg + tem = zi(i,k+1) * (zfac**pfac) * abs( moninq_fac) + +!!!! CHANGES FOR HEIGHT-DEPENDENT K ADJUSTMENT, WANG W + if(useshape .ge. 1) then + sz2h=(ZI(I,K+1)-ZL(I,1))/(HPBL(I)-ZL(I,1)) + sz2h=max(sz2h,zfmin) + sz2h=min(sz2h,1.0) + zfac=(1.0-sz2h)**pfac +! smax=0.148 !! max value of this shape function + smax=0.148 !! max value of this shape function + hmax=0.333 !! roughly height if max K + skmax=hmax*(1.0-hmax)**pfac + sksfc=min(ZI(I,2)/HPBL(I),0.05) ! surface layer top, 0.05H or ZI(2) (Zi(1)=0) + sksfc=sksfc*(1-sksfc)**pfac + + zfac=max(zfac,zfmin) + ashape=max(ABS(moninq_fac),0.2) ! should not be smaller than 0.2, otherwise too much adjustment(?) + if(useshape ==1) then + ashape=( 1.0 - ((sz2h*zfac/smax)**0.25) + & *( 1.0 - ashape ) ) + tem = zi(i,k+1) * (zfac) * ashape + endif + + if (useshape == 2) then !only adjus K that is > K_surface_top + ashape1=1.0 + if (skmax > sksfc) ashape1=(skmax*ashape-sksfc)/ + & (skmax-sksfc) + skminusk0=ZI(I,K+1)*zfac - HPBL(i)*sksfc + tem = zi(i,k+1) * (zfac) ! no adjustment + if (skminusk0 > 0) then ! only adjust K which is > surface top K + tem = skminusk0*ashape1 + HPBL(i)*sksfc + endif + endif + endif ! endif useshape>1 +!!!! END OF CHAGES , WANG W + + + if(pblflg(i)) then + tem1 = vk * wscaleu(i) * tem +! dku(i,k) = xkzmo(i,k) + tem1 +! dkt(i,k) = xkzo(i,k) + tem1 * prinv(i) + dku(i,k) = tem1 + dkt(i,k) = tem1 * prinv(i) + else + tem1 = vk * wscale(i) * tem +! dku(i,k) = xkzmo(i,k) + tem1 +! dkt(i,k) = xkzo(i,k) + tem1 * prinv(i) + dku(i,k) = tem1 + dkt(i,k) = tem1 * prinv(i) + endif + dku(i,k) = min(dku(i,k),dkmax) + dku(i,k) = max(dku(i,k),xkzmo(i,k)) + dkt(i,k) = min(dkt(i,k),dkmax) + dkt(i,k) = max(dkt(i,k),xkzo(i,k)) + dktx(i,k)= dkt(i,k) + endif + enddo !K loop + +! possible modification of first guess DKU, under certain conditions +! (1) this applies only to columns over water + + IF(islimsk(i).eq.0)then ! sea only + +! (2) alpha test +! if alpha < 0, find alpha for each column and do the loop again +! if alpha > 0, we are finished + + + if(alpha.lt.0)then ! variable alpha test + +! k-level of layer around 500 m + kLOC = INT(WSPM(i,2)) +! print *,' kLOC ',kLOC,' KPBL ',KPBL(I) + +! (3) only do this IF KPBL(I) >= kLOC. Otherwise, we are finished, with DKU as +! if alpha = +1 + + if(KPBL(I).gt.kLOC)then + + xDKU = DKU(i,kLOC) ! Km at k-level +! (4) DKU check. +! WSPM(i,1) is the KM cap for the 500-m level. +! if DKU at 500-m level < WSPM(i,1), do not limit Km ANYWHERE. Alpha = +! abs(alpha). No need to recalc. +! if DKU at 500-m level > WSPM(i,1), then alpha = WSPM(i,1)/xDKU for entire +! column + if(xDKU.ge.WSPM(i,1)) then ! ONLY if DKU at 500-m exceeds cap, otherwise already done + + WSPM(i,3) = WSPM(i,1)/xDKU ! ratio of cap to Km at k-level, store in WSPM(i,3) + !WSPM(i,4) = amin1(WSPM(I,3),1.0) ! this is new column alpha. cap at 1. ! should never be needed + WSPM(i,4) = min(WSPM(I,3),1.0) ! this is new column alpha. cap at 1. ! should never be needed + !! recalculate K capped by WSPM(i,1) + do k = 1, kmpbl + if(k < kpbl(i)) then +! zfac = max((1.-(zi(i,k+1)-zl(i,1))/ +! 1 (hpbl(i)-zl(i,1))), zfmin) + zfac = max((1.-zi(i,k+1)/hpbl(i)), zfmin) + ! tem = zi(i,k+1) * (zfac**pfac) + tem = zi(i,k+1) * (zfac**pfac) * WSPM(i,4) + + +!!!! CHANGES FOR HEIGHT-DEPENDENT K ADJUSTMENT, WANG W + if(useshape .ge. 1) then + sz2h=(ZI(I,K+1)-ZL(I,1))/(HPBL(I)-ZL(I,1)) + sz2h=max(sz2h,zfmin) + sz2h=min(sz2h,1.0) + zfac=(1.0-sz2h)**pfac + smax=0.148 !! max value of this shape function + hmax=0.333 !! roughly height if max K + skmax=hmax*(1.0-hmax)**pfac + sksfc=min(ZI(I,2)/HPBL(I),0.05) ! surface layer top, 0.05H or ZI(2) (Zi(1)=0) + sksfc=sksfc*(1-sksfc)**pfac + + zfac=max(zfac,zfmin) + ashape=max(WSPM(i,4),0.2) !! adjustment coef should not smaller than 0.2 + if(useshape ==1) then + ashape=( 1.0 - ((sz2h*zfac/smax)**0.25) + & *( 1.0 - ashape ) ) + tem = zi(i,k+1) * (zfac) * ashape +! if(k ==5) write(0,*)'min alf, height-depend alf',WSPM(i,4),ashape + endif ! endif useshape=1 + + if (useshape == 2) then !only adjus K that is > K_surface_top + ashape1=1.0 + if (skmax > sksfc) ashape1=(skmax*ashape-sksfc)/ + & (skmax-sksfc) + + skminusk0=ZI(I,K+1)*zfac - HPBL(i)*sksfc + tem = zi(i,k+1) * (zfac) ! no adjustment +! if(k ==5) write(0,*)'before, dku,ashape,ashpe1', +! & tem*wscaleu(i)*vk,ashape,ashape1 + if (skminusk0 > 0) then ! only adjust K which is > surface top K + tem = skminusk0*ashape1 + HPBL(i)*sksfc + endif +! if(k ==5)write(0,*) +! & 'after,dku,k_sfc,skmax,sksfc,zi(2),hpbl' +! & ,tem*wscaleu(i)*vk,WSCALEU(I)*VK*HPBL(i)*sksfc, skmax, +! & sksfc,ZI(I,2),HPBL(I) + + endif ! endif useshape=2 + endif ! endif useshape>1 +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + + if(pblflg(i)) then + tem1 = vk * wscaleu(i) * tem +! dku(i,k) = xkzmo(i,k) + tem1 +! dkt(i,k) = xkzo(i,k) + tem1 * prinv(i) + dku(i,k) = tem1 + dkt(i,k) = tem1 * prinv(i) + else + tem1 = vk * wscale(i) * tem +! dku(i,k) = xkzmo(i,k) + tem1 +! dkt(i,k) = xkzo(i,k) + tem1 * prinv(i) + dku(i,k) = tem1 + dkt(i,k) = tem1 * prinv(i) + endif + dku(i,k) = min(dku(i,k),dkmax) + dku(i,k) = max(dku(i,k),xkzmo(i,k)) + dkt(i,k) = min(dkt(i,k),dkmax) + dkt(i,k) = max(dkt(i,k),xkzo(i,k)) + dktx(i,k)= dkt(i,k) + endif + enddo !K loop + endif ! xDKU.ge.WSPM(i,1) + endif ! KPBL(I).ge.kLOC + endif ! alpha < 0 + endif ! islimsk=0 + + enddo !I loop + ENDIF !AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA + +! +! compute diffusion coefficients based on local scheme above pbl +!> ## Compute diffusion coefficients above the PBL top +!! Diffusion coefficients above the PBL top are computed as a function of local stability (gradient Richardson number), shear, and a length scale from Louis (1979) \cite louis_1979 : +!! \f[ +!! K_{m,h}=l^2f_{m,h}(Ri_g)\left|\frac{\partial U}{\partial z}\right| +!! \f] +!! The functions used (\f$f_{m,h}\f$) depend on the local stability. First, the gradient Richardson number is calculated as +!! \f[ +!! Ri_g=\frac{\frac{g}{T}\frac{\partial \theta_v}{\partial z}}{\frac{\partial U}{\partial z}^2} +!! \f] +!! where \f$U\f$ is the horizontal wind. For the unstable case (\f$Ri_g < 0\f$), the Richardson number-dependent functions are given by +!! \f[ +!! f_h(Ri_g) = 1 + \frac{8\left|Ri_g\right|}{1 + 1.286\sqrt{\left|Ri_g\right|}}\\ +!! \f] +!! \f[ +!! f_m(Ri_g) = 1 + \frac{8\left|Ri_g\right|}{1 + 1.746\sqrt{\left|Ri_g\right|}}\\ +!! \f] +!! For the stable case, the following formulas are used +!! \f[ +!! f_h(Ri_g) = \frac{1}{\left(1 + 5Ri_g\right)^2}\\ +!! \f] +!! \f[ +!! Pr = \frac{K_h}{K_m} = 1 + 2.1Ri_g +!! \f] +!! The source for the formulas used for the Richardson number-dependent functions is unclear. They are different than those used in Hong and Pan (1996) \cite hong_and_pan_1996 as the previous documentation suggests. They follow equation 14 of Louis (1979) \cite louis_1979 for the unstable case, but it is unclear where the values of the coefficients \f$b\f$ and \f$c\f$ from that equation used in this scheme originate. Finally, the length scale, \f$l\f$ is calculated according to the following formula from Hong and Pan (1996) \cite hong_and_pan_1996 +!! \f[ +!! \frac{1}{l} = \frac{1}{kz} + \frac{1}{l_0}\\ +!! \f] +!! \f[ +!! or\\ +!! \f] +!! \f[ +!! l=\frac{l_0kz}{l_0+kz} +!! \f] +!! where \f$l_0\f$ is currently 30 m for stable conditions and 150 m for unstable. Finally, the diffusion coefficients are kept in a range bounded by the background diffusion and the maximum allowable values. + do k = 1, km1 + do i=1,im + if(k >= kpbl(i)) then + bvf2 = g*bf(i,k)*ti(i,k) + ri = max(bvf2/shr2(i,k),rimin) + zk = vk*zi(i,k+1) + if(ri < 0.) then ! unstable regime + rl2 = zk*rlamun/(rlamun+zk) + dk = rl2*rl2*sqrt(shr2(i,k)) + sri = sqrt(-ri) +! dku(i,k) = xkzmo(i,k) + dk*(1+8.*(-ri)/(1+1.746*sri)) +! dkt(i,k) = xkzo(i,k) + dk*(1+8.*(-ri)/(1+1.286*sri)) + dku(i,k) = dk*(1+8.*(-ri)/(1+1.746*sri)) + dkt(i,k) = dk*(1+8.*(-ri)/(1+1.286*sri)) + else ! stable regime + rl2 = zk*rlam/(rlam+zk) +!! tem = rlam * sqrt(0.01*prsi(i,k)) +!! rl2 = zk*tem/(tem+zk) + dk = rl2*rl2*sqrt(shr2(i,k)) + tem1 = dk/(1+5.*ri)**2 +! + if(k >= kpblx(i)) then + prnum = 1.0 + 2.1*ri + prnum = min(prnum,prmax) + else + prnum = 1.0 + endif +! dku(i,k) = xkzmo(i,k) + tem1 * prnum +! dkt(i,k) = xkzo(i,k) + tem1 + dku(i,k) = tem1 * prnum + dkt(i,k) = tem1 + endif +! + dku(i,k) = min(dku(i,k),dkmax) + dku(i,k) = max(dku(i,k),xkzmo(i,k)) + dkt(i,k) = min(dkt(i,k),dkmax) + dkt(i,k) = max(dkt(i,k),xkzo(i,k)) +! + endif +! + enddo + enddo +! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! compute components for mass flux mixing by large thermals +!> ## If the PBL is convective, call the mass flux scheme to replace the countergradient terms. +!! If the PBL is convective, the updraft properties are initialized to be the same as the state variables and the subroutine mfpbl is called. + do k = 1, km + do i = 1, im + if(pcnvflg(i)) then + tcko(i,k) = t1(i,k) + ucko(i,k) = u1(i,k) + vcko(i,k) = v1(i,k) + xmf(i,k) = 0. + endif + enddo + enddo + do kk = 1, ntrac + do k = 1, km + do i = 1, im + if(pcnvflg(i)) then + qcko(i,k,kk) = q1(i,k,kk) + endif + enddo + enddo + enddo +!> For details of the mfpbl subroutine, step into its documentation ::mfpbl + call mfpbl(im,ix,km,ntrac,dt2,pcnvflg, + & zl,zi,thvx,q1,t1,u1,v1,hpbl,kpbl, + & sflux,ustar,wstar,xmf,tcko,qcko,ucko,vcko) +! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! compute diffusion coefficients for cloud-top driven diffusion +! if the condition for cloud-top instability is met, +! increase entrainment flux at cloud top +! +!> ## Compute enhanced diffusion coefficients related to stratocumulus-topped PBLs +!! If a stratocumulus layer has been identified in the PBL, the diffusion coefficients in the PBL are modified in the following way. +!! +!! -# First, the criteria for CTEI is checked, using the threshold from equation 13 of Macvean and Mason (1990) \cite macvean_and_mason_1990. If the criteria is met, the cloud top diffusion is increased: +!! \f[ +!! K_h^{Sc} = -c\frac{\Delta F_R}{\rho c_p}\frac{1}{\frac{\partial \theta_v}{\partial z}} +!! \f] +!! where the constant \f$c\f$ is set to 0.2 if the CTEI criterion is not met and 1.0 if it is. +!! +!! -# Calculate the diffusion coefficients due to stratocumulus mixing according to equation 5 in Lock et al. (2000) \cite lock_et_al_2000 for every level below the stratocumulus top using the characteristic stratocumulus velocity scale previously calculated. The diffusion coefficient for momentum is calculated assuming a constant inverse Prandtl number of 0.75. + do i = 1, im + if(scuflg(i)) then + k = krad(i) + tem = thetae(i,k) - thetae(i,k+1) + tem1 = qtx(i,k) - qtx(i,k+1) + if (tem > 0. .and. tem1 > 0.) then + cteit= cp*tem/(hvap*tem1) + if(cteit > actei) rent(i) = rentf2 + endif + endif + enddo + do i = 1, im + if(scuflg(i)) then + k = krad(i) + tem1 = max(bf(i,k),tdzmin) + ckt(i,k) = -rent(i)*radmin(i)/tem1 + cku(i,k) = ckt(i,k) + endif + enddo +! + do k = 1, kmpbl + do i=1,im + if(scuflg(i) .and. k < krad(i)) then + tem1=hrad(i)-zd(i) + tem2=zi(i,k+1)-tem1 + if(tem2 > 0.) then + ptem= tem2/zd(i) + if(ptem.ge.1.) ptem= 1. + ptem= tem2*ptem*sqrt(1.-ptem) + ckt(i,k) = radfac*vk*vrad(i)*ptem + cku(i,k) = 0.75*ckt(i,k) + ckt(i,k) = max(ckt(i,k),dkmin) + ckt(i,k) = min(ckt(i,k),dkmax) + cku(i,k) = max(cku(i,k),dkmin) + cku(i,k) = min(cku(i,k),dkmax) + endif + endif + enddo + enddo +! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! +!> After \f$K_h^{Sc}\f$ has been determined from the surface to the top of the stratocumulus layer, it is added to the value for the diffusion coefficient calculated previously using surface-based mixing [see equation 6 of Lock et al. (2000) \cite lock_et_al_2000 ]. + do k = 1, kmpbl + do i=1,im + if(scuflg(i)) then + ! dkt(i,k) = dkt(i,k)+ckt(i,k) + ! dku(i,k) = dku(i,k)+cku(i,k) + !! if K needs to be adjusted by alpha, then no need to add this term + if(alpha .ge. 0.0) dkt(i,k) = dkt(i,k)+ckt(i,k) + if(alpha .ge. 0.0) dku(i,k) = dku(i,k)+cku(i,k) + + dkt(i,k) = min(dkt(i,k),dkmax) + dku(i,k) = min(dku(i,k),dkmax) + endif + enddo + enddo +! +! compute tridiagonal matrix elements for heat and moisture +! +!> ## Solve for the temperature and moisture tendencies due to vertical mixing. +!! The tendencies of heat, moisture, and momentum due to vertical diffusion are calculated using a two-part process. First, a solution is obtained using an implicit time-stepping scheme, then the time tendency terms are "backed out". The tridiagonal matrix elements for the implicit solution for temperature and moisture are prepared in this section, with differing algorithms depending on whether the PBL was convective (substituting the mass flux term for counter-gradient term), unstable but not convective (using the computed counter-gradient terms), or stable (no counter-gradient terms). + do i=1,im + ad(i,1) = 1. + a1(i,1) = t1(i,1) + beta(i) * heat(i) + a2(i,1) = q1(i,1,1) + beta(i) * evap(i) + enddo + + if(ntrac >= 2) then + do k = 2, ntrac + is = (k-1) * km + do i = 1, im + a2(i,1+is) = q1(i,1,k) + enddo + enddo + endif +! + do k = 1,km1 + do i = 1,im + dtodsd = dt2/del(i,k) + dtodsu = dt2/del(i,k+1) + dsig = prsl(i,k)-prsl(i,k+1) + rdz = rdzt(i,k) + tem1 = dsig * dkt(i,k) * rdz + dsdz2 = tem1 * rdz + au(i,k) = -dtodsd*dsdz2 + al(i,k) = -dtodsu*dsdz2 +! + if(pcnvflg(i) .and. k < kpbl(i)) then + tem2 = dsig * rdz + ptem = 0.5 * tem2 * xmf(i,k) + ptem1 = dtodsd * ptem + ptem2 = dtodsu * ptem + ad(i,k) = ad(i,k)-au(i,k)-ptem1 + ad(i,k+1) = 1.-al(i,k)+ptem2 + au(i,k) = au(i,k)-ptem1 + al(i,k) = al(i,k)+ptem2 + ptem = tcko(i,k) + tcko(i,k+1) + dsdzt = tem1 * gocp + a1(i,k) = a1(i,k)+dtodsd*dsdzt-ptem1*ptem + a1(i,k+1) = t1(i,k+1)-dtodsu*dsdzt+ptem2*ptem + ptem = qcko(i,k,1) + qcko(i,k+1,1) + a2(i,k) = a2(i,k) - ptem1 * ptem + a2(i,k+1) = q1(i,k+1,1) + ptem2 * ptem + elseif(ublflg(i) .and. k < kpbl(i)) then + ptem1 = dsig * dktx(i,k) * rdz + tem = 1.0 / hpbl(i) + dsdzt = tem1 * gocp - ptem1 * hgamt(i) * tem + dsdzq = - ptem1 * hgamq(i) * tem + ad(i,k) = ad(i,k)-au(i,k) + ad(i,k+1) = 1.-al(i,k) + a1(i,k) = a1(i,k)+dtodsd*dsdzt + a1(i,k+1) = t1(i,k+1)-dtodsu*dsdzt + a2(i,k) = a2(i,k)+dtodsd*dsdzq + a2(i,k+1) = q1(i,k+1,1)-dtodsu*dsdzq + else + ad(i,k) = ad(i,k)-au(i,k) + ad(i,k+1) = 1.-al(i,k) + dsdzt = tem1 * gocp + a1(i,k) = a1(i,k)+dtodsd*dsdzt + a1(i,k+1) = t1(i,k+1)-dtodsu*dsdzt + a2(i,k+1) = q1(i,k+1,1) + endif +! + enddo + enddo +! + if(ntrac >= 2) then + do kk = 2, ntrac + is = (kk-1) * km + do k = 1, km1 + do i = 1, im + if(pcnvflg(i) .and. k < kpbl(i)) then + dtodsd = dt2/del(i,k) + dtodsu = dt2/del(i,k+1) + dsig = prsl(i,k)-prsl(i,k+1) + tem = dsig * rdzt(i,k) + ptem = 0.5 * tem * xmf(i,k) + ptem1 = dtodsd * ptem + ptem2 = dtodsu * ptem + tem1 = qcko(i,k,kk) + qcko(i,k+1,kk) + a2(i,k+is) = a2(i,k+is) - ptem1*tem1 + a2(i,k+1+is)= q1(i,k+1,kk) + ptem2*tem1 + else + a2(i,k+1+is) = q1(i,k+1,kk) + endif + enddo + enddo + enddo + endif +! +! solve tridiagonal problem for heat and moisture +! +!> The tridiagonal system is solved by calling the internal ::tridin subroutine. + call tridin99(im,km,ntrac,al,ad,au,a1,a2,au,a1,a2) + +! +! recover tendencies of heat and moisture +! +!> After returning with the solution, the tendencies for temperature and moisture are recovered. + do k = 1,km + do i = 1,im + ttend = (a1(i,k)-t1(i,k)) * rdt + qtend = (a2(i,k)-q1(i,k,1))*rdt + tau(i,k) = tau(i,k)+ttend + rtg(i,k,1) = rtg(i,k,1)+qtend + dtsfc(i) = dtsfc(i)+cont*del(i,k)*ttend + dqsfc(i) = dqsfc(i)+conq*del(i,k)*qtend + enddo + enddo + if(ntrac >= 2) then + do kk = 2, ntrac + is = (kk-1) * km + do k = 1, km + do i = 1, im + qtend = (a2(i,k+is)-q1(i,k,kk))*rdt + rtg(i,k,kk) = rtg(i,k,kk)+qtend + enddo + enddo + enddo + endif +! +! compute tke dissipation rate +! +!> ## Calculate heating due to TKE dissipation and add to the tendency for temperature +!! Following Han et al. (2015) \cite han_et_al_2015 , turbulence dissipation contributes to the tendency of temperature in the following way. First, turbulence dissipation is calculated by equation 17 of Han et al. (2015) \cite han_et_al_2015 for the PBL and equation 16 for the surface layer. + if(dspheat) then +! + do k = 1,km1 + do i = 1,im + diss(i,k) = dku(i,k)*shr2(i,k)-g*ti(i,k)*dkt(i,k)*bf(i,k) +! diss(i,k) = dku(i,k)*shr2(i,k) + enddo + enddo +! +! add dissipative heating at the first model layer +! +!> Next, the temperature tendency is updated following equation 14. + do i = 1,im + tem = govrth(i)*sflux(i) + tem1 = tem + stress(i)*spd1(i)/zl(i,1) + tem2 = 0.5 * (tem1+diss(i,1)) + tem2 = max(tem2, 0.) + ttend = tem2 / cp + if (alpha .gt. 0.0) then + tau(i,1) = tau(i,1)+0.5*ttend + else + tau(i,1) = tau(i,1)+0.7*ttend ! in HWRF/HMON, use 0.7 + endif + enddo +! +! add dissipative heating above the first model layer +! + do k = 2,km1 + do i = 1,im + tem = 0.5 * (diss(i,k-1)+diss(i,k)) + tem = max(tem, 0.) + ttend = tem / cp + tau(i,k) = tau(i,k) + 0.5*ttend + enddo + enddo +! + endif +! +! compute tridiagonal matrix elements for momentum +! +!> ## Solve for the horizontal momentum tendencies and add them to the output tendency terms +!! As with the temperature and moisture tendencies, the horizontal momentum tendencies are calculated by solving tridiagonal matrices after the matrices are prepared in this section. + do i=1,im + ad(i,1) = 1.0 + beta(i) * stress(i) / spd1(i) + a1(i,1) = u1(i,1) + a2(i,1) = v1(i,1) + enddo +! + do k = 1,km1 + do i=1,im + dtodsd = dt2/del(i,k) + dtodsu = dt2/del(i,k+1) + dsig = prsl(i,k)-prsl(i,k+1) + rdz = rdzt(i,k) + tem1 = dsig*dku(i,k)*rdz + dsdz2 = tem1 * rdz + au(i,k) = -dtodsd*dsdz2 + al(i,k) = -dtodsu*dsdz2 +! + if(pcnvflg(i) .and. k < kpbl(i)) then + tem2 = dsig * rdz + ptem = 0.5 * tem2 * xmf(i,k) + ptem1 = dtodsd * ptem + ptem2 = dtodsu * ptem + ad(i,k) = ad(i,k)-au(i,k)-ptem1 + ad(i,k+1) = 1.-al(i,k)+ptem2 + au(i,k) = au(i,k)-ptem1 + al(i,k) = al(i,k)+ptem2 + ptem = ucko(i,k) + ucko(i,k+1) + a1(i,k) = a1(i,k) - ptem1 * ptem + a1(i,k+1) = u1(i,k+1) + ptem2 * ptem + ptem = vcko(i,k) + vcko(i,k+1) + a2(i,k) = a2(i,k) - ptem1 * ptem + a2(i,k+1) = v1(i,k+1) + ptem2 * ptem + else + ad(i,k) = ad(i,k)-au(i,k) + ad(i,k+1) = 1.-al(i,k) + a1(i,k+1) = u1(i,k+1) + a2(i,k+1) = v1(i,k+1) + endif +! + enddo + enddo +! +! solve tridiagonal problem for momentum +! + call tridi299(im,km,al,ad,au,a1,a2,au,a1,a2) +! +! recover tendencies of momentum +! +!> Finally, the tendencies are recovered from the tridiagonal solutions. + do k = 1,km + do i = 1,im + utend = (a1(i,k)-u1(i,k))*rdt + vtend = (a2(i,k)-v1(i,k))*rdt + du(i,k) = du(i,k) + utend + dv(i,k) = dv(i,k) + vtend + dusfc(i) = dusfc(i) + conw*del(i,k)*utend + dvsfc(i) = dvsfc(i) + conw*del(i,k)*vtend +! +! for dissipative heating for ecmwf model +! +! tem1 = 0.5*(a1(i,k)+u1(i,k)) +! tem2 = 0.5*(a2(i,k)+v1(i,k)) +! diss(i,k) = -(tem1*utend+tem2*vtend) +! diss(i,k) = max(diss(i,k),0.) +! ttend = diss(i,k) / cp +! tau(i,k) = tau(i,k) + ttend +! + enddo + enddo +! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! + do i = 1, im + hpbl(i) = hpblx(i) + kpbl(i) = kpblx(i) + enddo +! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + return + end subroutine hedmf_hafs_run + +!> @} + +c----------------------------------------------------------------------- +!> \ingroup PBL +!! \brief Routine to solve the tridiagonal system to calculate temperature and moisture at \f$ t + \Delta t \f$; part of two-part process to calculate time tendencies due to vertical diffusion. +!! +!! Origin of subroutine unknown. + subroutine tridi299(l,n,cl,cm,cu,r1,r2,au,a1,a2) +cc + use machine , only : kind_phys + implicit none + integer k,n,l,i + real(kind=kind_phys) fk +cc + real(kind=kind_phys) cl(l,2:n),cm(l,n),cu(l,n-1),r1(l,n),r2(l,n), & + & au(l,n-1),a1(l,n),a2(l,n) +c----------------------------------------------------------------------- + do i=1,l + fk = 1./cm(i,1) + au(i,1) = fk*cu(i,1) + a1(i,1) = fk*r1(i,1) + a2(i,1) = fk*r2(i,1) + enddo + do k=2,n-1 + do i=1,l + fk = 1./(cm(i,k)-cl(i,k)*au(i,k-1)) + au(i,k) = fk*cu(i,k) + a1(i,k) = fk*(r1(i,k)-cl(i,k)*a1(i,k-1)) + a2(i,k) = fk*(r2(i,k)-cl(i,k)*a2(i,k-1)) + enddo + enddo + do i=1,l + fk = 1./(cm(i,n)-cl(i,n)*au(i,n-1)) + a1(i,n) = fk*(r1(i,n)-cl(i,n)*a1(i,n-1)) + a2(i,n) = fk*(r2(i,n)-cl(i,n)*a2(i,n-1)) + enddo + do k=n-1,1,-1 + do i=1,l + a1(i,k) = a1(i,k)-au(i,k)*a1(i,k+1) + a2(i,k) = a2(i,k)-au(i,k)*a2(i,k+1) + enddo + enddo +c----------------------------------------------------------------------- + return + end subroutine tridi299 +c----------------------------------------------------------------------- +!> \ingroup PBL +!! \brief Routine to solve the tridiagonal system to calculate u- and v-momentum at \f$ t + \Delta t \f$; part of two-part process to calculate time tendencies due to vertical diffusion. +!! +!! Origin of subroutine unknown. + subroutine tridin99(l,n,nt,cl,cm,cu,r1,r2,au,a1,a2) +cc + use machine , only : kind_phys + implicit none + integer is,k,kk,n,nt,l,i + real(kind=kind_phys) fk(l) +cc + real(kind=kind_phys) cl(l,2:n), cm(l,n), cu(l,n-1), & + & r1(l,n), r2(l,n*nt), & + & au(l,n-1), a1(l,n), a2(l,n*nt), & + & fkk(l,2:n-1) +c----------------------------------------------------------------------- + do i=1,l + fk(i) = 1./cm(i,1) + au(i,1) = fk(i)*cu(i,1) + a1(i,1) = fk(i)*r1(i,1) + enddo + do k = 1, nt + is = (k-1) * n + do i = 1, l + a2(i,1+is) = fk(i) * r2(i,1+is) + enddo + enddo + do k=2,n-1 + do i=1,l + fkk(i,k) = 1./(cm(i,k)-cl(i,k)*au(i,k-1)) + au(i,k) = fkk(i,k)*cu(i,k) + a1(i,k) = fkk(i,k)*(r1(i,k)-cl(i,k)*a1(i,k-1)) + enddo + enddo + do kk = 1, nt + is = (kk-1) * n + do k=2,n-1 + do i=1,l + a2(i,k+is) = fkk(i,k)*(r2(i,k+is)-cl(i,k)*a2(i,k+is-1)) + enddo + enddo + enddo + do i=1,l + fk(i) = 1./(cm(i,n)-cl(i,n)*au(i,n-1)) + a1(i,n) = fk(i)*(r1(i,n)-cl(i,n)*a1(i,n-1)) + enddo + do k = 1, nt + is = (k-1) * n + do i = 1, l + a2(i,n+is) = fk(i)*(r2(i,n+is)-cl(i,n)*a2(i,n+is-1)) + enddo + enddo + do k=n-1,1,-1 + do i=1,l + a1(i,k) = a1(i,k) - au(i,k)*a1(i,k+1) + enddo + enddo + do kk = 1, nt + is = (kk-1) * n + do k=n-1,1,-1 + do i=1,l + a2(i,k+is) = a2(i,k+is) - au(i,k)*a2(i,k+is+1) + enddo + enddo + enddo +c----------------------------------------------------------------------- + return + end subroutine tridin99 + +!> @} + + end module hedmf_hafs diff --git a/physics/moninedmf_hafs.meta b/physics/moninedmf_hafs.meta new file mode 100644 index 000000000..bc1461ada --- /dev/null +++ b/physics/moninedmf_hafs.meta @@ -0,0 +1,526 @@ +[ccpp-arg-table] + name = hedmf_hafs_init + type = scheme +[moninq_fac] + standard_name = atmosphere_diffusivity_coefficient_factor + long_name = multiplicative constant for atmospheric diffusivities + units = none + dimensions = () + type = real + kind = kind_phys + 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 = hedmf_hafs_run + type = scheme +[ix] + standard_name = horizontal_dimension + long_name = horizontal dimension + units = count + dimensions = () + type = integer + intent = in + optional = F +[im] + standard_name = horizontal_loop_extent + long_name = horizontal loop extent + units = count + dimensions = () + type = integer + intent = in + optional = F +[km] + standard_name = vertical_dimension + long_name = vertical layer dimension + units = count + dimensions = () + type = integer + intent = in + optional = F +[ntrac] + standard_name = number_of_vertical_diffusion_tracers + long_name = number of tracers to diffuse vertically + units = count + dimensions = () + type = integer + intent = in + optional = F +[ntcw] + standard_name = index_for_liquid_cloud_condensate + long_name = cloud condensate index in tracer array + units = index + dimensions = () + type = integer + intent = in + optional = F +[dv] + standard_name = tendency_of_y_wind_due_to_model_physics + long_name = updated tendency of the y wind + units = m s-2 + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[du] + standard_name = tendency_of_x_wind_due_to_model_physics + long_name = updated tendency of the x wind + units = m s-2 + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[tau] + standard_name = tendency_of_air_temperature_due_to_model_physics + long_name = updated tendency of the temperature + units = K s-1 + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[rtg] + standard_name = tendency_of_vertically_diffused_tracer_concentration + long_name = updated tendency of the tracers due to vertical diffusion in PBL scheme + units = kg kg-1 s-1 + dimensions = (horizontal_dimension,vertical_dimension,number_of_vertical_diffusion_tracers) + type = real + kind = kind_phys + intent = inout + optional = F +[u1] + standard_name = x_wind + long_name = x component of layer 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 = y component of layer 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 = layer mean air temperature + units = K + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[q1] + standard_name = vertically_diffused_tracer_concentration + long_name = tracer concentration diffused by PBL scheme + units = kg kg-1 + dimensions = (horizontal_dimension,vertical_dimension,number_of_vertical_diffusion_tracers) + type = real + kind = kind_phys + intent = in + optional = F +[swh] + standard_name = tendency_of_air_temperature_due_to_shortwave_heating_on_radiation_time_step + long_name = total sky shortwave heating rate + units = K s-1 + dimensions = (horizontal_dimension,adjusted_vertical_layer_dimension_for_radiation) + type = real + kind = kind_phys + intent = in + optional = F +[hlw] + standard_name = tendency_of_air_temperature_due_to_longwave_heating_on_radiation_time_step + long_name = total sky longwave heating rate + units = K s-1 + dimensions = (horizontal_dimension,adjusted_vertical_layer_dimension_for_radiation) + type = real + kind = kind_phys + intent = in + optional = F +[xmu] + standard_name = zenith_angle_temporal_adjustment_factor_for_shortwave_fluxes + long_name = zenith angle temporal adjustment factor for shortwave + units = none + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[psk] + standard_name = dimensionless_exner_function_at_lowest_model_interface + long_name = dimensionless Exner function at the surface interface + units = none + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[rbsoil] + 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 +[zorl] + standard_name = surface_roughness_length + long_name = surface roughness length in cm + units = cm + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[u10m] + standard_name = x_wind_at_10m + long_name = x component of wind at 10 m + units = m s-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[v10m] + standard_name = y_wind_at_10m + long_name = y component of wind at 10 m + units = m s-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[fm] + standard_name = Monin_Obukhov_similarity_function_for_momentum + long_name = Monin-Obukhov similarity function for momentum + units = none + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[fh] + standard_name = Monin_Obukhov_similarity_function_for_heat + long_name = Monin-Obukhov similarity function for heat + units = none + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[tsea] + standard_name = surface_skin_temperature + long_name = surface skin temperature + units = K + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[heat] + standard_name = kinematic_surface_upward_sensible_heat_flux + long_name = kinematic surface upward sensible heat flux + units = K m s-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[evap] + standard_name = kinematic_surface_upward_latent_heat_flux + long_name = kinematic surface upward latent heat flux + units = kg kg-1 m s-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[stress] + standard_name = surface_wind_stress + long_name = surface wind stress + units = m2 s-2 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[spd1] + standard_name = wind_speed_at_lowest_model_layer + long_name = wind speed at lowest model level + units = m s-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[kpbl] + standard_name = vertical_index_at_top_of_atmosphere_boundary_layer + long_name = PBL top model level index + units = index + dimensions = (horizontal_dimension) + type = integer + intent = out + optional = F +[prsi] + standard_name = air_pressure_at_interface + long_name = air pressure at model layer interfaces + 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 = pres(k) - pres(k+1) + units = Pa + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[prsl] + standard_name = air_pressure + long_name = mean 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 = Exner function at layers + units = none + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[phii] + standard_name = geopotential_at_interface + long_name = geopotential at model layer interfaces + 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 = geopotential at model layer centers + units = m2 s-2 + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[delt] + standard_name = time_step_for_physics + long_name = time step for physics + units = s + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F +[dspheat] + standard_name = flag_TKE_dissipation_heating + long_name = flag for using TKE dissipation heating + units = flag + dimensions = () + type = logical + intent = in + optional = F +[dusfc] + standard_name = instantaneous_surface_x_momentum_flux + long_name = x momentum flux + units = Pa + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = out + optional = F +[dvsfc] + standard_name = instantaneous_surface_y_momentum_flux + long_name = y momentum flux + units = Pa + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = out + optional = F +[dtsfc] + standard_name = instantaneous_surface_upward_sensible_heat_flux + long_name = surface upward sensible heat flux + units = W m-2 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = out + optional = F +[dqsfc] + standard_name = instantaneous_surface_upward_latent_heat_flux + long_name = surface upward latent heat flux + units = W m-2 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = out + optional = F +[hpbl] + standard_name = atmosphere_boundary_layer_thickness + long_name = PBL thickness + units = m + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = out + optional = F +[hgamt] + standard_name = countergradient_mixing_term_for_temperature + long_name = countergradient mixing term for temperature + units = K + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[hgamq] + standard_name = countergradient_mixing_term_for_water_vapor + long_name = countergradient mixing term for water vapor + units = kg kg-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[dkt] + standard_name = atmosphere_heat_diffusivity + long_name = diffusivity for heat + units = m2 s-1 + dimensions = (horizontal_dimension,vertical_dimension_minus_one) + type = real + kind = kind_phys + intent = out + optional = F +[kinver] + standard_name = index_of_highest_temperature_inversion + long_name = index of highest temperature inversion + units = index + dimensions = (horizontal_dimension) + type = integer + intent = in + optional = F +[xkzm_m] + standard_name = atmosphere_momentum_diffusivity_background + long_name = background value of momentum diffusivity + units = m2 s-1 + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F +[xkzm_h] + standard_name = atmosphere_heat_diffusivity_background + long_name = background value of heat diffusivity + units = m2 s-1 + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F +[xkzm_s] + standard_name = diffusivity_background_sigma_level + long_name = sigma level threshold for background diffusivity + units = none + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F +[lprnt] + standard_name = flag_print + long_name = flag for printing diagnostics to output + units = flag + dimensions = () + type = logical + intent = in + optional = F +[ipr] + standard_name = horizontal_index_of_printed_column + long_name = horizontal index of printed column + units = index + dimensions = () + type = integer + intent = in + optional = F +[xkzminv] + standard_name = atmosphere_heat_diffusivity_background_maximum + long_name = maximum background value of heat diffusivity + units = m2 s-1 + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F +[moninq_fac] + standard_name = atmosphere_diffusivity_coefficient_factor + long_name = multiplicative constant for atmospheric diffusivities + units = none + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F +[islimsk] + 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 +[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 From ed16475af3e9710368b388614cf0c26d4830d24f Mon Sep 17 00:00:00 2001 From: "Chunxi.Zhang-NOAA" Date: Wed, 22 Jan 2020 16:42:07 +0000 Subject: [PATCH 3/6] Fixed the bugs related to the closure of shallow convection. --- physics/cu_ntiedtke.F90 | 40 +++++++++++++++++----------------------- physics/cu_ntiedtke.meta | 22 ++++++++++++++++++++-- 2 files changed, 37 insertions(+), 25 deletions(-) diff --git a/physics/cu_ntiedtke.F90 b/physics/cu_ntiedtke.F90 index 8e42ebdd4..c06f3ecc7 100644 --- a/physics/cu_ntiedtke.F90 +++ b/physics/cu_ntiedtke.F90 @@ -148,7 +148,7 @@ end subroutine cu_ntiedtke_finalize !----------------------------------------------------------------------- ! level 1 subroutine 'tiecnvn' !----------------------------------------------------------------- - subroutine cu_ntiedtke_run(pu,pv,pt,pqv,pqvf,ptf,clw,poz,pzz,prsl,prsi,pomg, & + subroutine cu_ntiedtke_run(pu,pv,pt,pqv,tdi,qvdi,pqvf,ptf,clw,poz,pzz,prsl,prsi,pomg, & evap,hfx,zprecc,lmask,lq,ix,km,dt,dx,kbot,ktop,kcnv,& ktrac,ud_mf,dd_mf,dt_mf,cnvw,cnvc,errmsg,errflg) !----------------------------------------------------------------- @@ -162,13 +162,9 @@ subroutine cu_ntiedtke_run(pu,pv,pt,pqv,pqvf,ptf,clw,poz,pzz,prsl,prsi,pomg, & integer, dimension( lq ), intent(in) :: lmask real(kind=kind_phys), dimension( lq ), intent(in ) :: evap, hfx, dx real(kind=kind_phys), dimension( ix , km ), intent(inout) :: pu, pv, pt, pqv - real(kind=kind_phys), dimension( ix , km ), intent(in ) :: poz, prsl, pomg, pqvf, ptf + real(kind=kind_phys), dimension( ix , km ), intent(in ) :: tdi, qvdi, poz, prsl, pomg, pqvf, ptf real(kind=kind_phys), dimension( ix , km+1 ), intent(in ) :: pzz, prsi - ! DH* TODO - check dimensions of clw, ktrac+2 seems to be smaller - ! than the actual dimensions (ok as long as only indices 1 and 2 - ! are accessed here, and as long as these contain what is expected); - ! better to expand into the cloud-ice and cloud-water components *DH - real(kind=kind_phys), dimension( ix , km, ktrac+2 ), intent(inout ) :: clw + real(kind=kind_phys), dimension( ix , km, ktrac ), intent(inout ) :: clw integer, dimension( lq ), intent(out) :: kbot, ktop, kcnv real(kind=kind_phys), dimension( lq ), intent(out) :: zprecc @@ -188,7 +184,7 @@ subroutine cu_ntiedtke_run(pu,pv,pt,pqv,pqvf,ptf,clw,poz,pzz,prsl,prsi,pomg, & real(kind=kind_phys) ztp1(lq,km), zqp1(lq,km), ztu(lq,km), zqu(lq,km),& & zlu(lq,km), zlude(lq,km), zmfu(lq,km), zmfd(lq,km), zmfude_rate(lq,km),& & zqsat(lq,km), zrain(lq) - real(kind=kind_phys) pcen(lq,km,ktrac),ptenc(lq,km,ktrac) + real(kind=kind_phys) pcen(lq,km,ktrac-2),ptenc(lq,km,ktrac-2) integer icbot(lq), ictop(lq), ktype(lq), lndj(lq) logical locum(lq) @@ -246,9 +242,9 @@ subroutine cu_ntiedtke_run(pu,pv,pt,pqv,pqvf,ptf,clw,poz,pzz,prsl,prsi,pomg, & zqs = min(0.5,zqs) zcor = 1./(1.-vtmpc1*zqs) zqsat(j,k1)=zqs*zcor - pqte(j,k1)=pqvf(j,k) + pqte(j,k1)=pqvf(j,k)+(pqv(j,k)-qvdi(j,k))/ztmst zqq(j,k1) =pqte(j,k1) - ptte(j,k1)=ptf(j,k) + ptte(j,k1)=ptf(j,k)+(pt(j,k)-tdi(j,k))/ztmst ztt(j,k1) =ptte(j,k1) ud_mf(j,k1)=0. dd_mf(j,k1)=0. @@ -258,7 +254,7 @@ subroutine cu_ntiedtke_run(pu,pv,pt,pqv,pqvf,ptf,clw,poz,pzz,prsl,prsi,pomg, & end do end do - do n=1,ktrac + do n=1,ktrac-2 do k=1,km k1=km-k+1 do j=1,lq @@ -289,7 +285,7 @@ subroutine cu_ntiedtke_run(pu,pv,pt,pqv,pqvf,ptf,clw,poz,pzz,prsl,prsi,pomg, & & zqp1, pum1, pvm1, pverv, zqsat,& & pqhfl, ztmst, pap, paph, pgeo, & & ptte, pqte, pvom, pvol, prsfc,& - & pssfc, locum, ktrac, pcen, ptenc,& + & pssfc, locum, ktrac-2, pcen, ptenc,& & ktype, icbot, ictop, ztu, zqu, & & zlu, zlude, zmfu, zmfd, zrain,& & pcte, phhfl, lndj, pgeoh, zmfude_rate, dx) @@ -314,7 +310,7 @@ subroutine cu_ntiedtke_run(pu,pv,pt,pqv,pqvf,ptf,clw,poz,pzz,prsl,prsi,pomg, & pt(j,k) = ztp1(j,k1)+(ptte(j,k1)-ztt(j,k1))*ztmst pqv(j,k)= zqp1(j,k1)+(pqte(j,k1)-zqq(j,k1))*ztmst ud_mf(j,k)= zmfu(j,k1)*ztmst - dd_mf(j,k)= zmfd(j,k1)*ztmst + dd_mf(j,k)= -zmfd(j,k1)*ztmst dt_mf(j,k)= zmfude_rate(j,k1)*ztmst cnvw(j,k) = zlude(j,k1)*ztmst*g/(prsi(j,k)-prsi(j,k+1)) cnvc(j,k) = 0.04 * log(1. + 675. * ud_mf(j,k)) @@ -344,16 +340,14 @@ subroutine cu_ntiedtke_run(pu,pv,pt,pqv,pqvf,ptf,clw,poz,pzz,prsl,prsi,pomg, & end do endif ! - if (ktrac > 0) then - do n=1,ktrac - do k=1,km - k1=km-k+1 - do j=1,lq - clw(j,k,n+2)=pcen(j,k,n)+ptenc(j,k1,n)*ztmst - end do - end do - end do - end if +! do n=1,ktrac-2 +! do k=1,km +! k1=km-k+1 +! do j=1,lq +! clw(j,k,n+2)=pcen(j,k,n)+ptenc(j,k1,n)*ztmst +! end do +! end do +! end do ! return end subroutine cu_ntiedtke_run diff --git a/physics/cu_ntiedtke.meta b/physics/cu_ntiedtke.meta index da9219c10..4208b6e46 100644 --- a/physics/cu_ntiedtke.meta +++ b/physics/cu_ntiedtke.meta @@ -80,6 +80,24 @@ kind = kind_phys intent = inout optional = F +[tdi] + 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 +[qvdi] + standard_name = water_vapor_specific_humidity + long_name = water vapor specific humidity + units = kg kg-1 + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = in + optional = F [pqvf] standard_name = moisture_tendency_due_to_dynamics long_name = moisture tendency due to dynamics only @@ -254,8 +272,8 @@ intent = out optional = F [ktrac] - standard_name = number_of_total_tracers - long_name = number of total tracers + standard_name = number_of_tracers_for_convective_transport + long_name = number of tracers for convective transport units = count dimensions = () type = integer From 7db1d7ec8d2832cf372bc2d86d83a0d391c3cd0d Mon Sep 17 00:00:00 2001 From: "Chunxi.Zhang-NOAA" Date: Thu, 23 Jan 2020 15:40:02 +0000 Subject: [PATCH 4/6] Revised the code the new Tiedtke scheme --- physics/cu_ntiedtke.F90 | 57 +++++++++++++++++++++++++++++------------ 1 file changed, 40 insertions(+), 17 deletions(-) diff --git a/physics/cu_ntiedtke.F90 b/physics/cu_ntiedtke.F90 index c06f3ecc7..156e75c70 100644 --- a/physics/cu_ntiedtke.F90 +++ b/physics/cu_ntiedtke.F90 @@ -184,13 +184,13 @@ subroutine cu_ntiedtke_run(pu,pv,pt,pqv,tdi,qvdi,pqvf,ptf,clw,poz,pzz,prsl,prsi, real(kind=kind_phys) ztp1(lq,km), zqp1(lq,km), ztu(lq,km), zqu(lq,km),& & zlu(lq,km), zlude(lq,km), zmfu(lq,km), zmfd(lq,km), zmfude_rate(lq,km),& & zqsat(lq,km), zrain(lq) - real(kind=kind_phys) pcen(lq,km,ktrac-2),ptenc(lq,km,ktrac-2) + real(kind=kind_phys),allocatable :: pcen(:,:,:),ptenc(:,:,:) integer icbot(lq), ictop(lq), ktype(lq), lndj(lq) logical locum(lq) ! real(kind=kind_phys) ztmst,fliq,fice,ztc,zalf,tt - integer i,j,k,k1,n,km1 + integer i,j,k,k1,n,km1,ktracer real(kind=kind_phys) ztpp1 real(kind=kind_phys) zew,zqs,zcor ! @@ -254,16 +254,33 @@ subroutine cu_ntiedtke_run(pu,pv,pt,pqv,tdi,qvdi,pqvf,ptf,clw,poz,pzz,prsl,prsi, end do end do - do n=1,ktrac-2 - do k=1,km - k1=km-k+1 - do j=1,lq - pcen(j,k1,n) = clw(j,k,n+2) - ptenc(j,k1,n)= 0. + if(ktrac > 2) then + ktracer = ktrac - 2 + allocate(pcen(lq,km,ktracer)) + allocate(ptenc(lq,km,ktracer)) + do n=1,ktracer + do k=1,km + k1=km-k+1 + do j=1,lq + pcen(j,k1,n) = clw(j,k,n+2) + ptenc(j,k1,n)= 0. + end do end do end do - end do - + else + ktracer = 2 + allocate(pcen(lq,km,ktracer)) + allocate(ptenc(lq,km,ktracer)) + do n=1,ktracer + do k=1,km + do j=1,lq + pcen(j,k,n) = 0. + ptenc(j,k,n)= 0. + end do + end do + end do + end if + ! print *, "pgeo=",pgeo(1,:) ! print *, "pgeoh=",pgeoh(1,:) ! print *, "pap=",pap(1,:) @@ -285,7 +302,7 @@ subroutine cu_ntiedtke_run(pu,pv,pt,pqv,tdi,qvdi,pqvf,ptf,clw,poz,pzz,prsl,prsi, & zqp1, pum1, pvm1, pverv, zqsat,& & pqhfl, ztmst, pap, paph, pgeo, & & ptte, pqte, pvom, pvol, prsfc,& - & pssfc, locum, ktrac-2, pcen, ptenc,& + & pssfc, locum, ktracer, pcen, ptenc,& & ktype, icbot, ictop, ztu, zqu, & & zlu, zlude, zmfu, zmfd, zrain,& & pcte, phhfl, lndj, pgeoh, zmfude_rate, dx) @@ -339,15 +356,21 @@ subroutine cu_ntiedtke_run(pu,pv,pt,pqv,tdi,qvdi,pqvf,ptf,clw,poz,pzz,prsl,prsi, end do end do endif + ! -! do n=1,ktrac-2 -! do k=1,km -! k1=km-k+1 -! do j=1,lq -! clw(j,k,n+2)=pcen(j,k,n)+ptenc(j,k1,n)*ztmst +! Currently, vertical mixing of tracers are turned off +! if(ktrac > 2) then +! do n=1,ktrac-2 +! do k=1,km +! k1=km-k+1 +! do j=1,lq +! clw(j,k,n+2)=pcen(j,k,n)+ptenc(j,k1,n)*ztmst +! end do ! end do ! end do -! end do +! end if + deallocate(pcen) + deallocate(ptenc) ! return end subroutine cu_ntiedtke_run From 8223afed5e66b0e124702eab47feca896543a1ee Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Thu, 23 Jan 2020 11:50:13 -0700 Subject: [PATCH 5/6] physics/module_mp_thompson.F90: bugfix, remove threaded computation/read of lookup tables --- physics/module_mp_thompson.F90 | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index b1ca6ba07..dfaea5c2f 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -924,11 +924,6 @@ SUBROUTINE thompson_init(nwfa2d, nifa2d, nwfa, nifa, & call cpu_time(stime) -!$OMP parallel num_threads(threads) - -!$OMP sections - -!$OMP section !> - Call qr_acr_qg() to create rain collecting graupel & graupel collecting rain table if (mpirank==mpiroot) write(0,*) ' creating rain collecting graupel table' call cpu_time(stime) @@ -936,7 +931,6 @@ SUBROUTINE thompson_init(nwfa2d, nifa2d, nwfa, nifa, & call cpu_time(etime) if (mpirank==mpiroot) print '("Computing rain collecting graupel table took ",f10.3," seconds.")', etime-stime -!$OMP section !> - Call qr_acr_qs() to create rain collecting snow & snow collecting rain table if (mpirank==mpiroot) write (*,*) ' creating rain collecting snow table' call cpu_time(stime) @@ -944,10 +938,6 @@ SUBROUTINE thompson_init(nwfa2d, nifa2d, nwfa, nifa, & call cpu_time(etime) if (mpirank==mpiroot) print '("Computing rain collecting snow table took ",f10.3," seconds.")', etime-stime -!$OMP end sections - -!$OMP end parallel - !> - Call freezeh2o() to create cloud water and rain freezing (Bigg, 1953) table if (mpirank==mpiroot) write(0,*) ' creating freezing of water drops table' call cpu_time(stime) From 4c7dcaa8ae1e5465c5358647e75c239c6dafb30c Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Mon, 27 Jan 2020 10:08:39 -0700 Subject: [PATCH 6/6] Add missing updates from IPD physics commit 7ffe6471c20404091fbbf8f321fbb9ee84a4f36d --- physics/module_gfdl_cloud_microphys.F90 | 2 +- physics/module_sf_noahmp_glacier.f90 | 0 physics/module_sf_noahmplsm.f90 | 0 physics/noahmp_tables.f90 | 0 physics/sfc_noahmp_drv.f | 0 5 files changed, 1 insertion(+), 1 deletion(-) mode change 100755 => 100644 physics/module_sf_noahmp_glacier.f90 mode change 100755 => 100644 physics/module_sf_noahmplsm.f90 mode change 100755 => 100644 physics/noahmp_tables.f90 mode change 100755 => 100644 physics/sfc_noahmp_drv.f diff --git a/physics/module_gfdl_cloud_microphys.F90 b/physics/module_gfdl_cloud_microphys.F90 index 01ab4655c..5750d27fd 100644 --- a/physics/module_gfdl_cloud_microphys.F90 +++ b/physics/module_gfdl_cloud_microphys.F90 @@ -3320,7 +3320,7 @@ subroutine fall_speed (ktop, kbot, den, qs, qi, qg, ql, tk, vts, vti, vtg) else tc (k) = tk (k) - tice vti (k) = (3. + log10 (qi (k) * den (k))) * (tc (k) * (aa * tc (k) + bb) + cc) + dd * tc (k) + ee - vti (k) = vi0 * exp (log_10 * vti (k)) * 0.8 + vti (k) = vi0 * exp (log_10 * vti (k)) * 0.9 vti (k) = min (vi_max, max (vf_min, vti (k))) endif enddo diff --git a/physics/module_sf_noahmp_glacier.f90 b/physics/module_sf_noahmp_glacier.f90 old mode 100755 new mode 100644 diff --git a/physics/module_sf_noahmplsm.f90 b/physics/module_sf_noahmplsm.f90 old mode 100755 new mode 100644 diff --git a/physics/noahmp_tables.f90 b/physics/noahmp_tables.f90 old mode 100755 new mode 100644 diff --git a/physics/sfc_noahmp_drv.f b/physics/sfc_noahmp_drv.f old mode 100755 new mode 100644