diff --git a/physics/module_samf.F90 b/physics/module_samf.F90 new file mode 100644 index 000000000..ef1e8114a --- /dev/null +++ b/physics/module_samf.F90 @@ -0,0 +1,195 @@ +! ######################################################################################### +!> \file module_samf.F90 +!! This file contains default parameters for the Scale-Aware mass flux Convection scheme. +!! +!! These parameters can be overwritten by the host by calling ty_cfg_samf%setup(). +!! +! ######################################################################################### +module module_samf + use machine, only : kind_phys + implicit none + + public ty_cfg_samf_shal, ty_cfg_samf_deep + +! ######################################################################################### +!! \section arg_table_ty_cfg_samf_shal Argument Table +!! \htmlinclude ty_cfg_samf_shal.html +!! +! ######################################################################################### + type ty_cfg_samf_shal + ! Scheme parameters (default) + real(kind_phys) :: d0 = 0.001 ! From han_et_al_2017 equation 8 + real(kind_phys) :: cm = 1.0 ! Enhancement factor in entrainment rates for momentum + real(kind_phys) :: cq = 1.0 ! + real(kind_phys) :: clamd = 0.1 ! + real(kind_phys) :: tkemx = 0.65 ! Maximum value for turbulent kinetic energy (TKE). + real(kind_phys) :: tkemn = 0.05 ! Minimum value for turbulent kinetic energy (TKE). + real(kind_phys) :: cthk = 200. ! Maximum cloud depth for shallow convection. + real(kind_phys) :: cthkmn = 0. ! Minimum cloud depth for shallow convection. + real(kind_phys) :: dthk = 25. ! With entrainment, recalculate the LFC as the first level + ! where buoyancy is positive. The difference in pressure + ! levels between LFCs calculated with/without entrainment + ! must be less than a threshold (dthk), in hPa. + real(kind_phys) :: sfclfac = 0.2 ! + real(kind_phys) :: rhcrt = 0.75 ! Minimum mean relative humidity in the subcloud layers + ! for convection to be triggered. + real(kind_phys) :: cinpcrmx = 180. ! Maximum pressure depth between parcel source level + ! for convection to be triggered. + real(kind_phys) :: cinpcrmn = 120. ! Minimum pressure depth between parcel source level + ! for convection to be triggered. + real(kind_phys) :: cinacrmx = -120 ! Minimum CIN for convection to be triggered. + real(kind_phys) :: shevf = 2.0 ! Enhancement factor for evaporation. + real(kind_phys) :: dtmax = 10800 ! Maximum timestep + real(kind_phys) :: dtmin = 600 ! Minimum timestep + real(kind_phys) :: bb1 = 4.0 ! From han_et_al_2017 equation 7 + real(kind_phys) :: bb2 = 0.8 ! From han_et_al_2017 equation 7 + real(kind_phys) :: csmf = 0.2 ! Parameter for updraft velocity calculation. + real(kind_phys) :: tkcrt = 2. ! + real(kind_phys) :: cmxfac = 15. ! + real(kind_phys) :: betaw = .03 ! From han_et_al_2017 equation 6 + real(kind_phys) :: dxcrtc0 = 9.e3 ! Parameter for scale-aware rain conversion calculation. + real(kind_phys) :: escav = 0.8 ! + real(kind_phys) :: tf = 233.16! + real(kind_phys) :: tcr = 263.16! + real(kind_phys) :: dxcrtas = 30.e3 ! + real(kind_phys) :: sigmind = 0.01 ! progsigma: Minimum area-fraction for deep clouds + real(kind_phys) :: sigmins = 0.03 ! progsigma: Minimum area-fraction for shallow clouds + real(kind_phys) :: sigminm = 0.01 ! progsigma: Minimum area-fraction for mid-level clouds + ! Derived scheme parameters + real(kind_phys) :: tcrf ! + real(kind_phys) :: dtke ! Valid range for TKE + real(kind_phys) :: cinacrmn ! Set by host (NOTE. No default value is provided!!!!) + contains + procedure, public :: setup => setup_param_shal + end type ty_cfg_samf_shal + +! ######################################################################################### +!! \section arg_table_ty_cfg_samf_deep Argument Table +!! \htmlinclude ty_cfg_samf_deep.html +!! +! ######################################################################################### + type ty_cfg_samf_deep + ! Scheme parameters (default) + real(kind_phys) :: d0 = 0.001 ! From han_et_al_2017 equation 8 + real(kind_phys) :: cm = 1.0 ! Enhancement factor in entrainment rates for momentum + real(kind_phys) :: cq = 1.0 ! + real(kind_phys) :: clamd = 0.03 ! + real(kind_phys) :: tkemx = 0.65 ! Maximum value for turbulent kinetic energy (TKE). + real(kind_phys) :: tkemn = 0.05 ! Minimum value for turbulent kinetic energy (TKE). + real(kind_phys) :: clamca = 0.03 ! + real(kind_phys) :: cthk = 200. ! Maximum cloud depth for deep convection. + real(kind_phys) :: dthk = 25. ! With entrainment, recalculate the LFC as the first level + ! where buoyancy is positive. The difference in pressure + ! levels between LFCs calculated with/without entrainment + ! must be less than a threshold (dthk), in hPa. + real(kind_phys) :: sfclfac = 0.2 ! + real(kind_phys) :: rhcrt = 0.75 ! Minimum mean relative humidity in the subcloud layers + ! for convection to be triggered. + real(kind_phys) :: cinpcrmx = 180. ! Maximum pressure depth between parcel source level + ! for convection to be triggered. + real(kind_phys) :: cinpcrmn = 120. ! Minimum pressure depth between parcel source level + ! for convection to be triggered. + real(kind_phys) :: cinacrmx = -120 ! Minimum CIN for convection to be triggered. + real(kind_phys) :: cinacrmn = -80 ! Maximum CIN for convection to be triggered. + real(kind_phys) :: bb1 = 4.0 ! From han_et_al_2017 equation 7 + real(kind_phys) :: bb2 = 0.8 ! From han_et_al_2017 equation 7 + real(kind_phys) :: csmf = 0.2 ! Parameter for updraft velocity calculation. + real(kind_phys) :: tkcrt = 2. ! + real(kind_phys) :: cmxfac = 15. ! + real(kind_phys) :: betaw = .03 ! From han_et_al_2017 equation 6 + ! Derived scheme parameters + real(kind_phys) :: dtke ! Valid range for TKE + contains + procedure, public :: setup => setup_param_deep + end type ty_cfg_samf_deep +contains + ! ######################################################################################### + ! Procedure (type-bound) to setup parameters needed by scheme. + ! ######################################################################################### + subroutine setup_param_shal(this,d0,cm,cq,clamd,tkemx,tkemn,cthk,cthkmn,dthk,sfclfac, & + rhcrt,cinpcrmx,cinpcrmn,cinacrmx,cinacrmn,shevf,dtmax,dtmin,bb1,bb2,csmf,tkcrt, & + cmxfac,betaw,dxcrtc0,h1,dxcrtas,sigmind,sigmins,sigminm,escav,tf,tcr) + ! + class(ty_cfg_samf_shal), intent(inout) :: this + real(kind_phys), intent(in), optional :: d0,cm,cq,clamd,tkemx,tkemn,cthk,cthkmn,dthk, & + sfclfac,rhcrt,cinpcrmx,cinpcrmn,cinacrmx,cinacrmn,shevf,dtmax,dtmin,bb1,bb2,csmf, & + tkcrt,cmxfac,betaw,dxcrtc0,h1,dxcrtas,sigmind,sigmins,sigminm,escav,tf,tcr + + ! If provided, override default values for parameters + if (present(d0 )) this%d0 = d0 + if (present(cm )) this%cm = cm + if (present(cq )) this%cq = cq + if (present(clamd )) this%clamd = clamd + if (present(tkemx )) this%tkemx = tkemx + if (present(tkemn )) this%tkemn = tkemn + if (present(cthk )) this%cthk = cthk + if (present(cthkmn )) this%cthkmn = cthkmn + if (present(dthk )) this%dthk = dthk + if (present(sfclfac )) this%sfclfac = sfclfac + if (present(rhcrt )) this%rhcrt = rhcrt + if (present(cinpcrmx)) this%cinpcrmx = cinpcrmx + if (present(cinpcrmn)) this%cinpcrmn = cinpcrmn + if (present(cinacrmx)) this%cinacrmx = cinacrmx + if (present(cinacrmn)) this%cinacrmn = cinacrmn + if (present(shevf )) this%shevf = shevf + if (present(dtmax )) this%dtmax = dtmax + if (present(dtmin )) this%dtmin = dtmin + if (present(bb1 )) this%bb1 = bb1 + if (present(bb2 )) this%bb2 = bb2 + if (present(csmf )) this%csmf = csmf + if (present(tkcrt )) this%tkcrt = tkcrt + if (present(cmxfac )) this%cmxfac = cmxfac + if (present(betaw )) this%betaw = betaw + if (present(dxcrtc0 )) this%dxcrtc0 = dxcrtc0 + if (present(dxcrtas )) this%dxcrtas = dxcrtas + if (present(sigmind )) this%sigmind = sigmind + if (present(sigmins )) this%sigmins = sigmins + if (present(sigminm )) this%sigminm = sigminm + if (present(escav )) this%escav = escav + if (present(tf )) this%tf = tf + if (present(tcr )) this%tcr = tcr + + ! Compute derived parameters + this%tcrf = 1._kind_phys/(this%tcr-this%tf) + this%dtke = this%tkemx-this%tkemn + + end subroutine setup_param_shal + + ! ######################################################################################### + ! Procedure (type-bound) to setup parameters needed by scheme. + ! ######################################################################################### + subroutine setup_param_deep(this,d0,cm,cq,clamd,tkemx,tkemn,clamca,cthk,dthk,sfclfac, & + rhcrt,cinpcrmx,cinpcrmn,cinacrmx,cinacrmn,bb1,bb2,csmf,tkcrt,cmxfac,betaw) + ! + class(ty_cfg_samf_deep), intent(inout) :: this + real(kind_phys), intent(in), optional :: d0,cm,cq,clamd,tkemx,tkemn,clamca,cthk,dthk, & + sfclfac,rhcrt,cinpcrmx,cinpcrmn,cinacrmx,cinacrmn,bb1,bb2,csmf,tkcrt,cmxfac,betaw + + ! If provided, override default values for parameters + if (present(d0 )) this%d0 = d0 + if (present(cm )) this%cm = cm + if (present(cq )) this%cq = cq + if (present(clamd )) this%clamd = clamd + if (present(tkemx )) this%tkemx = tkemx + if (present(tkemn )) this%tkemn = tkemn + if (present(cthk )) this%cthk = cthk + if (present(dthk )) this%dthk = dthk + if (present(sfclfac )) this%sfclfac = sfclfac + if (present(rhcrt )) this%rhcrt = rhcrt + if (present(cinpcrmx)) this%cinpcrmx = cinpcrmx + if (present(cinpcrmn)) this%cinpcrmn = cinpcrmn + if (present(cinacrmx)) this%cinacrmx = cinacrmx + if (present(cinacrmn)) this%cinacrmn = cinacrmn + if (present(bb1 )) this%bb1 = bb1 + if (present(bb2 )) this%bb2 = bb2 + if (present(csmf )) this%csmf = csmf + if (present(tkcrt )) this%tkcrt = tkcrt + if (present(cmxfac )) this%cmxfac = cmxfac + if (present(betaw )) this%betaw = betaw + + ! Compute derived parameters + this%dtke = this%tkemx-this%tkemn + + end subroutine setup_param_deep + ! +end module module_samf diff --git a/physics/module_samf.meta b/physics/module_samf.meta new file mode 100644 index 000000000..ee86484fb --- /dev/null +++ b/physics/module_samf.meta @@ -0,0 +1,40 @@ +[ccpp-table-properties] + name = ty_cfg_samf_shal + type = ddt + dependencies = + +[ccpp-arg-table] + name = ty_cfg_samf_shal + type = ddt + +######################################################################## +[ccpp-table-properties] + name = ty_cfg_samf_deep + type = ddt + dependencies = + +[ccpp-arg-table] + name = ty_cfg_samf_deep + type = ddt + +######################################################################## +[ccpp-table-properties] + name = module_samf + type = module + dependencies = machine.F + +[ccpp-arg-table] + name = module_samf + type = module +[ty_cfg_samf_shal] + standard_name = ty_cfg_samf_shal + long_name = definition of type ty_cfg_samf_shal + units = DDT + dimensions = () + type = ty_cfg_samf_shal +[ty_cfg_samf_deep] + standard_name = ty_cfg_samf_deep + long_name = definition of type ty_cfg_samf_deep + units = DDT + dimensions = () + type = ty_cfg_samf_deep \ No newline at end of file diff --git a/physics/module_satmedmf.F90 b/physics/module_satmedmf.F90 new file mode 100644 index 000000000..64d0d943b --- /dev/null +++ b/physics/module_satmedmf.F90 @@ -0,0 +1,285 @@ +! ######################################################################################### +!> \file module_satmedmf.F90 +!! This file contains default real(kind_phys)s for the Scale-Aware mass flux Convection scheme. +!! +!! These real(kind_phys)s can be overwritten by the host by calling ty_cfg_satmedmf%setup(). +!! +! ######################################################################################### +module module_satmedmf + use machine, only : kind_phys + implicit none + + public ty_cfg_satmedmf,ty_cfg_satmedmfq + +! ######################################################################################### +!! \section arg_table_ty_cfg_satmedmf Argument Table +!! \htmlinclude ty_cfg_satmedmf.html +!! +! ######################################################################################### + type ty_cfg_satmedmf + ! Scheme parameters (default) + real(kind_phys) :: wfac = 7.0 + real(kind_phys) :: cfac = 4.5 + real(kind_phys) :: gamcrt = 3. + real(kind_phys) :: gamcrq = 0. + real(kind_phys) :: sfcfrac = 0.1 + real(kind_phys) :: vk = 0.4 + real(kind_phys) :: rimin = -100. + real(kind_phys) :: rbcr = 0.25 + real(kind_phys) :: zolcru = -0.02 + real(kind_phys) :: tdzmin = 1.e-3 + real(kind_phys) :: rlmn = 30. + real(kind_phys) :: rlmx = 500. + real(kind_phys) :: elmx = 500. + real(kind_phys) :: prmin = 0.25 + real(kind_phys) :: prmax = 4.0 + real(kind_phys) :: pr0 = 1.0 + real(kind_phys) :: prtke = 1.0 + real(kind_phys) :: prscu = 0.67 + real(kind_phys) :: f0 = 1.e-4 + real(kind_phys) :: crbmin = 0.15 + real(kind_phys) :: crbmax = 0.35 + real(kind_phys) :: tkmin = 1.e-9 + real(kind_phys) :: dspmax = 10.0 + real(kind_phys) :: dspfac = 0.5 + real(kind_phys) :: qmin = 1.e-8 + real(kind_phys) :: qlmin = 1.e-12 + real(kind_phys) :: zfmin = 1.e-8 + real(kind_phys) :: aphi5 = 5. + real(kind_phys) :: aphi16 = 16. + real(kind_phys) :: elmfac = 1.0 + real(kind_phys) :: elefac = 1.0 + real(kind_phys) :: cql = 100. + real(kind_phys) :: dw2min = 1.e-4 + real(kind_phys) :: dkmax = 1000. + real(kind_phys) :: xkgdx = 25000. + real(kind_phys) :: qlcr = 3.5e-5 + real(kind_phys) :: zstblmax = 2500. + real(kind_phys) :: xkzinv = 0.15 + real(kind_phys) :: h1 = 0.33333333 + real(kind_phys) :: ck0 = 0.4 + real(kind_phys) :: ch0 = 0.4 + real(kind_phys) :: ck1 = 0.15 + real(kind_phys) :: ch1 = 0.15 + real(kind_phys) :: ce0 = 0.4 + real(kind_phys) :: rchck = 1.5 + real(kind_phys) :: cdtn = 25 + contains + procedure, public :: setup => setup_satmedmf + end type ty_cfg_satmedmf + +! ######################################################################################### +!! \section arg_table_ty_cfg_satmedmfq Argument Table +!! \htmlinclude ty_cfg_satmedmfq.html +!! +! ######################################################################################### + type ty_cfg_satmedmfq + ! Scheme parameters (default) + real(kind_phys) :: bfac = 100. + real(kind_phys) :: wfac = 7.0 + real(kind_phys) :: cfac = 4.5 + real(kind_phys) :: gamcrt = 3. + real(kind_phys) :: gamcrq = 0. + real(kind_phys) :: sfcfrac = 0.1 + real(kind_phys) :: vk = 0.4 + real(kind_phys) :: rimin = -100. + real(kind_phys) :: slfac = 0.1 + real(kind_phys) :: rbcr = 0.25 + real(kind_phys) :: zolcru = -0.02 + real(kind_phys) :: tdzmin = 1.e-3 + real(kind_phys) :: rlmn = 30. + real(kind_phys) :: rlmn0 = 5. + real(kind_phys) :: rlmn1 = 5. + real(kind_phys) :: rlmn2 = 10. + real(kind_phys) :: prmin = 0.25 + real(kind_phys) :: prmax = 4.0 + real(kind_phys) :: pr0 = 1.0 + real(kind_phys) :: prtke = 1.0 + real(kind_phys) :: prscu = 0.67 + real(kind_phys) :: f0 = 1.e-4 + real(kind_phys) :: crbmin = 0.15 + real(kind_phys) :: crbmax = 0.35 + real(kind_phys) :: tkmin = 1.e-9 + real(kind_phys) :: tkbmx = 0.2 + real(kind_phys) :: dspmax = 10.0 + real(kind_phys) :: qmin = 1.e-8 + real(kind_phys) :: qlmin = 1.e-12 + real(kind_phys) :: zfmin = 1.e-8 + real(kind_phys) :: aphi5 = 5. + real(kind_phys) :: aphi16 = 16. + real(kind_phys) :: elmfac = 1.0 + real(kind_phys) :: elefac = 1.0 + real(kind_phys) :: cql = 100. + real(kind_phys) :: dw2min = 1.e-4 + real(kind_phys) :: dkmax = 1000. + real(kind_phys) :: xkgdx = 1000. + real(kind_phys) :: qlcr = 3.5e-5 + real(kind_phys) :: zstblmax = 2500. + real(kind_phys) :: xkinv1 = 0.15 + real(kind_phys) :: xkinv2 = 0.3 + real(kind_phys) :: h1 = 0.33333333 + real(kind_phys) :: hcrinv = 250. + real(kind_phys) :: vegflo = 0.1 + real(kind_phys) :: vegfup = 1.0 + real(kind_phys) :: z0lo = 0.1 + real(kind_phys) :: z0up = 1.0 + real(kind_phys) :: vc0 = 1.0 + real(kind_phys) :: zc0 = 1.0 + real(kind_phys) :: ck1 = 0.15 + real(kind_phys) :: ch1 = 0.15 + real(kind_phys) :: cs0 = 0.4 + real(kind_phys) :: csmf = 0.5 + real(kind_phys) :: rchck = 1.5 + real(kind_phys) :: ndt = 20 + real(kind_phys) :: ck0 = 0.4 + real(kind_phys) :: ch0 = 0.4 + real(kind_phys) :: ce0 = 0.4 + contains + procedure, public :: setup => setup_satmedmfq + end type ty_cfg_satmedmfq + +contains + ! ######################################################################################### + ! Procedure (type-bound) to setup parameters needed by scheme. + ! ######################################################################################### + subroutine setup_satmedmf(this,wfac,cfac,gamcrt,gamcrq,sfcfrac,vk,rimin,rbcr,zolcru, & + tdzmin,rlmn,rlmx,elmx,prmin,prmax,pr0,prtke,prscu,f0,crbmin,crbmax,tkmin,dspmax, & + dspfac,qmin,qlmin,zfmin,aphi5,aphi16,elmfac,elefac,cql,dw2min,dkmax,xkgdx,qlcr, & + zstblmax,xkzinv,h1,ck0,ch0,ck1,ch1,ce0,rchck,cdtn) + ! + class(ty_cfg_satmedmf), intent(inout) :: this + real(kind_phys),intent(in),optional :: wfac,cfac,gamcrt,gamcrq,sfcfrac,vk,rimin,rbcr, & + zolcru,tdzmin,rlmn,rlmx,elmx,prmin,prmax,pr0,prtke,prscu,f0,crbmin,crbmax,tkmin, & + dspmax,dspfac,qmin,qlmin,zfmin,aphi5,aphi16,elmfac,elefac,cql,dw2min,dkmax,xkgdx, & + qlcr,zstblmax,xkzinv,h1,rchck,ck0,ch0,ck1,ch1,ce0,cdtn + + ! If provided, override default values for parameters + if (present(wfac )) this%wfac = wfac + if (present(cfac )) this%cfac = cfac + if (present(gamcrt )) this%gamcrt = gamcrt + if (present(gamcrq )) this%gamcrq = gamcrq + if (present(sfcfrac )) this%sfcfrac = sfcfrac + if (present(vk )) this%vk = vk + if (present(rimin )) this%rimin = rimin + if (present(rbcr )) this%rbcr = rbcr + if (present(zolcru )) this%zolcru = zolcru + if (present(tdzmin )) this%tdzmin = tdzmin + if (present(rlmn )) this%rlmn = rlmn + if (present(rlmx )) this%rlmx = rlmx + if (present(elmx )) this%elmx = elmx + if (present(prmin )) this%prmin = prmin + if (present(prmax )) this%prmax = prmax + if (present(pr0 )) this%pr0 = pr0 + if (present(prtke )) this%prtke = prtke + if (present(prscu )) this%prscu = prscu + if (present(f0 )) this%f0 = f0 + if (present(crbmin )) this%crbmin = crbmin + if (present(crbmax )) this%crbmax = crbmax + if (present(tkmin )) this%tkmin = tkmin + if (present(dspmax )) this%dspmax = dspmax + if (present(dspfac )) this%dspfac = dspfac + if (present(qmin )) this%qmin = qmin + if (present(qlmin )) this%qlmin = qlmin + if (present(zfmin )) this%zfmin = zfmin + if (present(aphi5 )) this%aphi5 = aphi5 + if (present(aphi16 )) this%aphi16 = aphi16 + if (present(elmfac )) this%elmfac = elmfac + if (present(elefac )) this%elefac = elefac + if (present(cql )) this%cql = cql + if (present(dw2min )) this%dw2min = dw2min + if (present(dkmax )) this%dkmax = dkmax + if (present(xkgdx )) this%xkgdx = xkgdx + if (present(qlcr )) this%qlcr = qlcr + if (present(zstblmax)) this%zstblmax = zstblmax + if (present(xkzinv )) this%xkzinv = xkzinv + if (present(h1 )) this%h1 = h1 + if (present(ck0 )) this%ck0 = ck0 + if (present(ch0 )) this%ch0 = ch0 + if (present(ck1 )) this%ck1 = ck1 + if (present(ch1 )) this%ch1 = ch1 + if (present(ce0 )) this%ce0 = ce0 + if (present(rchck )) this%rchck = rchck + if (present(cdtn )) this%cdtn = cdtn + ! + end subroutine setup_satmedmf + + ! ######################################################################################### + ! Procedure (type-bound) to setup parameters needed by scheme. + ! ######################################################################################### + subroutine setup_satmedmfq(this,bfac,wfac,cfac,gamcrt,gamcrq,sfcfrac,vk,rimin,slfac,rbcr, & + zolcru,tdzmin,rlmn,rlmn0,rlmn1,rlmn2,prmin,prmax,pr0,prtke,prscu,f0,crbmin,crbmax, & + tkmin,tkbmx,dspmax,qmin,qlmin,zfmin,aphi5,aphi16,elmfac,elefac,cql,dw2min,dkmax, & + xkgdx,qlcr,zstblmax,xkinv1,xkinv2,h1,hcrinv,vegflo,vegfup,z0lo,z0up,vc0,zc0,ck1,ch1, & + cs0,ck0,ch0,ce0,csmf,rchck,ndt) + ! + class(ty_cfg_satmedmfq), intent(inout) :: this + real(kind_phys),intent(in),optional :: bfac,wfac,cfac,gamcrt,gamcrq,sfcfrac,vk,rimin, & + slfac,rbcr,zolcru,tdzmin,rlmn,rlmn0,rlmn1,rlmn2,prmin,prmax,pr0,prtke,prscu,f0, & + crbmin,crbmax,tkmin,tkbmx,dspmax,qmin,qlmin,zfmin,aphi5,aphi16,elmfac,elefac,cql, & + dw2min,dkmax,xkgdx,qlcr,zstblmax,xkinv1,xkinv2,h1,hcrinv,vegflo,vegfup,z0lo,z0up, & + vc0,zc0,ck1,ch1,cs0,ck0,ch0,ce0,csmf,rchck,ndt + + ! If provided, override default values for parameters + if (present(bfac)) this%bfac = bfac + if (present(wfac)) this%wfac = wfac + if (present(cfac)) this%cfac = cfac + if (present(gamcrt)) this%gamcrt = gamcrt + if (present(gamcrq)) this%gamcrq = gamcrq + if (present(sfcfrac)) this%sfcfrac = sfcfrac + if (present(vk)) this%vk = vk + if (present(rimin)) this%rimin = rimin + if (present(slfac)) this%slfac = slfac + if (present(rbcr)) this%rbcr = rbcr + if (present(zolcru)) this%zolcru = zolcru + if (present(tdzmin)) this%tdzmin = tdzmin + if (present(rlmn)) this%rlmn = rlmn + if (present(rlmn0)) this%rlmn0 = rlmn0 + if (present(rlmn1)) this%rlmn1 = rlmn1 + if (present(rlmn2)) this%rlmn2 = rlmn2 + if (present(prmin)) this%prmin = prmin + if (present(prmax)) this%prmax = prmax + if (present(pr0)) this%pr0 = pr0 + if (present(prtke)) this%prtke = prtke + if (present(prscu)) this%prscu = prscu + if (present(f0)) this%f0 = f0 + if (present(crbmin)) this%crbmin = crbmin + if (present(crbmax)) this%crbmax = crbmax + if (present(tkmin)) this%tkmin = tkmin + if (present(tkbmx)) this%tkbmx = tkbmx + if (present(dspmax)) this%dspmax = dspmax + if (present(qmin)) this%qmin = qmin + if (present(qlmin)) this%qlmin = qlmin + if (present(zfmin)) this%zfmin = zfmin + if (present(aphi5)) this%aphi5 = aphi5 + if (present(aphi16)) this%aphi16 = aphi16 + if (present(elmfac)) this%elmfac = elmfac + if (present(elefac)) this%elefac = elefac + if (present(cql)) this%cql = cql + if (present(dw2min)) this%dw2min = dw2min + if (present(dkmax)) this%dkmax = dkmax + if (present(xkgdx)) this%xkgdx = xkgdx + if (present(qlcr)) this%qlcr = qlcr + if (present(zstblmax)) this%zstblmax = zstblmax + if (present(xkinv1)) this%xkinv1 = xkinv1 + if (present(xkinv2)) this%xkinv2 = xkinv2 + if (present(h1)) this%h1 = h1 + if (present(hcrinv)) this%hcrinv = hcrinv + if (present(vegflo)) this%vegflo = vegflo + if (present(vegfup)) this%vegfup = vegfup + if (present(z0lo)) this%z0lo = z0lo + if (present(z0up)) this%z0up = z0up + if (present(vc0)) this%vc0 = vc0 + if (present(zc0)) this%zc0 = zc0 + if (present(ck1)) this%ck1 = ck1 + if (present(ch1)) this%ch1 = ch1 + if (present(cs0)) this%cs0 = cs0 + if (present(ck0)) this%ck0 = ck0 + if (present(ch0)) this%ch0 = ch0 + if (present(ce0)) this%ce0 = ce0 + if (present(csmf)) this%csmf = csmf + if (present(rchck)) this%rchck = rchck + if (present(ndt)) this%ndt = ndt + ! + end subroutine setup_satmedmfq + ! +end module module_satmedmf diff --git a/physics/module_satmedmf.meta b/physics/module_satmedmf.meta new file mode 100644 index 000000000..53f52f199 --- /dev/null +++ b/physics/module_satmedmf.meta @@ -0,0 +1,40 @@ +[ccpp-table-properties] + name = ty_cfg_satmedmf + type = ddt + dependencies = + +[ccpp-arg-table] + name = ty_cfg_satmedmf + type = ddt + +######################################################################## +[ccpp-table-properties] + name = ty_cfg_satmedmfq + type = ddt + dependencies = + +[ccpp-arg-table] + name = ty_cfg_satmedmfq + type = ddt + +######################################################################## +[ccpp-table-properties] + name = module_satmedmf + type = module + dependencies = machine.F + +[ccpp-arg-table] + name = module_satmedmf + type = module +[ty_cfg_satmedmf] + standard_name = ty_cfg_satmedmf + long_name = definition of type ty_cfg_satmedmf + units = DDT + dimensions = () + type = ty_cfg_satmedmf +[ty_cfg_satmedmfq] + standard_name = ty_cfg_satmedmfq + long_name = definition of type ty_cfg_satmedmfq + units = DDT + dimensions = () + type = ty_cfg_satmedmfq diff --git a/physics/samfshalcnv.f b/physics/samfshalcnv.f index d0bab05dd..f4fb41e3b 100644 --- a/physics/samfshalcnv.f +++ b/physics/samfshalcnv.f @@ -49,7 +49,7 @@ end subroutine samfshalcnv_init !! -# Calculate the tendencies of the state variables (per unit cloud base mass flux) and the cloud base mass flux. !! -# For the "feedback control", calculate updated values of the state variables by multiplying the cloud base mass flux and the tendencies calculated per unit cloud base mass flux from the static control. !! \section det_samfshalcnv GFS samfshalcnv Detailed Algorithm - subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & + subroutine samfshalcnv_run(cfg,im,km,itc,ntc,cliq,cp,cvap, & & eps,epsm1,fv,grav,hvap,rd,rv, & & t0c,delt,ntk,ntr,delp,first_time_step,restart, & & tmf,qmicro,progsigma, & @@ -61,9 +61,11 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & ! use machine , only : kind_phys use funcphys , only : fpvs + use module_samf, only: ty_cfg_samf_shal implicit none ! + type(ty_cfg_samf_shal), intent(in) :: cfg integer, intent(in) :: im, km, itc, ntc, ntk, ntr, ncloud integer, intent(in) :: islimsk(:) real(kind=kind_phys), intent(in) :: cliq, cp, cvap, & @@ -96,24 +98,20 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & ! local variables integer i,j,indx, k, kk, km1, n integer kpbl(im) -! - real(kind=kind_phys) clamd, tkemx, tkemn, dtke ! real(kind=kind_phys) dellat, - & c0l, d0, + & c0l, & desdt, dp, & dq, dqsdp, dqsdt, dt, - & dt2, dtmax, dtmin, - & dxcrt, dxcrtc0, + & dt2, + & dxcrt, & dv1h, dv2h, dv3h, & dz, dz1, e1, & el2orc, elocp, aafac, - & cm, cq, - & es, etah, h1, shevf, + & es, etah, ! & evfact, evfactl, & fact1, fact2, factor, - & cthk, cthkmn, dthk, - & gamma, pprime, betaw, tauadv, + & gamma, pprime, tauadv, & qlk, qrch, qs, & rfact, shear, tfac, & val, val1, val2, @@ -147,56 +145,20 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & c real(kind=kind_phys) crtlame, crtlamd ! - real(kind=kind_phys) cinpcr, cinpcrmx, cinpcrmn, - & cinacr, cinacrmx, cinacrmn, - & sfclfac, rhcrt, - & tkcrt, cmxfac + real(kind=kind_phys) cinpcr, cinacr + ! ! parameters for updraft velocity calculation - real(kind=kind_phys) bb1, bb2, csmf, wucb + real(kind=kind_phys) wucb cc ! parameters for prognostic sigma closure real(kind=kind_phys) omega_u(im,km),zdqca(im,km),tmfq(im,km), & omegac(im),zeta(im,km),dbyo1(im,km), & sigmab(im),qadv(im,km) - real(kind=kind_phys) gravinv,dxcrtas,invdelt,sigmind,sigmins, - & sigminm + real(kind=kind_phys) gravinv,invdelt logical flag_shallow,flag_mid -c physical parameters -! parameter(g=grav,asolfac=0.89) -! parameter(g=grav) -! parameter(elocp=hvap/cp, -! & el2orc=hvap*hvap/(rv*cp)) -! parameter(c0s=0.002,c1=5.e-4,d0=.01) -! parameter(d0=.01) - parameter(d0=.001) -! parameter(c0l=c0s*asolfac) -! -! asolfac: aerosol-aware parameter based on Lim & Hong (2012) -! asolfac= cx / c0s(=.002) -! cx = min([-0.7 ln(Nccn) + 24]*1.e-4, c0s) -! Nccn: CCN number concentration in cm^(-3) -! Until a realistic Nccn is provided, Nccns are assumed -! as Nccn=100 for sea and Nccn=1000 for land -! - parameter(cm=1.0,cq=1.0) -! parameter(fact1=(cvap-cliq)/rv,fact2=hvap/rv-fact1*t0c) - parameter(clamd=0.1,tkemx=0.65,tkemn=0.05) - parameter(dtke=tkemx-tkemn) - parameter(cthk=200.,cthkmn=0.,dthk=25.) - parameter(sfclfac=0.2,rhcrt=0.75) - parameter(cinpcrmx=180.,cinpcrmn=120.) -! shevf is an enhancing evaporation factor for shallow convection - parameter(cinacrmx=-120.,shevf=2.0) - parameter(dtmax=10800.,dtmin=600.) - parameter(bb1=4.0,bb2=0.8,csmf=0.2) - parameter(tkcrt=2.,cmxfac=10.) -! parameter(bet1=1.875,cd1=.506,f1=2.0,gam1=.5) - parameter(betaw=.03,dxcrtc0=9.e3) - parameter(h1=0.33333333) -! progsigma - parameter(dxcrtas=30.e3,sigmind=0.01,sigmins=0.03,sigminm=0.01) +! c local variables and arrays real(kind=kind_phys) pfld(im,km), to(im,km), qo(im,km), & uo(im,km), vo(im,km), qeso(im,km), @@ -206,7 +168,6 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & c variables for tracer wet deposition, real(kind=kind_phys), dimension(im,km,ntc) :: chem_c, chem_pw, & wet_dep - real(kind=kind_phys), parameter :: escav = 0.8 ! wet scavenging efficiency ! ! for updraft velocity calculation real(kind=kind_phys) wu2(im,km), buo(im,km), drag(im,km), @@ -240,10 +201,6 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & real(kind=kind_phys) tsumn(im), tsump(im), rtnp(im) ! logical do_aerosols, totflg, cnvflg(im), flg(im) -! - real(kind=kind_phys) tf, tcr, tcrf - parameter (tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf)) - c----------------------------------------------------------------------- ! @@ -260,10 +217,6 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & fact1 = (cvap-cliq)/rv fact2 = hvap/rv-fact1*t0c - if (.not.hwrf_samfshal) then - cinacrmn=-80. - endif - if (progsigma) then dxcrt=10.e3 else @@ -304,7 +257,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & kbot(i)=km+1 ktop(i)=0 endif - sfcpbl(i) = sfclfac * hpbl(i) + sfcpbl(i) = cfg%sfclfac * hpbl(i) rn(i)=0. kbcon(i)=km ktcon(i)=1 @@ -330,7 +283,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & kbot(i)=km+1 ktop(i)=0 endif - sfcpbl(i) = sfclfac * hpbl(i) + sfcpbl(i) = cfg%sfclfac * hpbl(i) rn(i)=0. kbcon(i)=km ktcon(i)=1 @@ -365,8 +318,8 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & ! !> - determine scale-aware rain conversion parameter decreasing with decreasing grid size do i=1,im - if(gdx(i) < dxcrtc0) then - tem = gdx(i) / dxcrtc0 + if(gdx(i) < cfg%dxcrtc0) then + tem = gdx(i) / cfg%dxcrtc0 tem1 = tem**3 c0(i) = c0(i) * tem1 endif @@ -378,7 +331,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & if(t1(i,k) > 273.16) then c0t(i,k) = c0(i) else - tem = d0 * (t1(i,k) - 273.16) + tem = cfg%d0 * (t1(i,k) - 273.16) tem1 = exp(tem) c0t(i,k) = c0(i) * tem1 endif @@ -745,8 +698,8 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & val2 = 1. tem = min(tem,val2) ptem = 1. - tem - ptem1= .5*(cinpcrmx-cinpcrmn) - cinpcr = cinpcrmx - ptem * ptem1 + ptem1= .5*(cfg%cinpcrmx-cfg%cinpcrmn) + cinpcr = cfg%cinpcrmx - ptem * ptem1 tem1 = pfld(i,kb(i)) - pfld(i,kbcon(i)) if(tem1 > cinpcr .and. @@ -815,7 +768,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & endif enddo ! -!> - if the mean relative humidity in the subcloud layers is less than a threshold value (rhcrt), convection is not triggered. +!> - if the mean relative humidity in the subcloud layers is less than a threshold value (cfg%rhcrt), convection is not triggered. ! do i = 1, im rhbar(i) = 0. @@ -835,7 +788,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & do i= 1, im if(cnvflg(i)) then rhbar(i) = rhbar(i) / sumx(i) - if(rhbar(i) < rhcrt) then + if(rhbar(i) < cfg%rhcrt) then cnvflg(i) = .false. endif endif @@ -886,23 +839,23 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & do i= 1, im if(cnvflg(i)) then tkemean(i) = tkemean(i) / sumx(i) - if(tkemean(i) > tkemx) then - clamt(i) = clam + clamd - else if(tkemean(i) < tkemn) then - clamt(i) = clam - clamd + if(tkemean(i) > cfg%tkemx) then + clamt(i) = clam + cfg%clamd + else if(tkemean(i) < cfg%tkemn) then + clamt(i) = clam - cfg%clamd else - tem = tkemx - tkemean(i) - tem1 = 1. - 2. * tem / dtke - clamt(i) = clam + clamd * tem1 + tem = cfg%tkemx - tkemean(i) + tem1 = 1. - 2. * tem / cfg%dtke + clamt(i) = clam + cfg%clamd * tem1 endif endif enddo ! do i=1,im if(cnvflg(i)) then - if(tkemean(i) > tkcrt) then - tem = 1. + tkemean(i)/tkcrt - tem1 = min(tem, cmxfac) + if(tkemean(i) > cfg%tkcrt) then + tem = 1. + tkemean(i)/cfg%tkcrt + tem1 = min(tem, cfg%cmxfac) clamt(i) = tem1 * clam endif endif @@ -1031,7 +984,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & (heo(i,k)+heo(i,k-1)))/factor dbyo(i,k) = hcko(i,k) - heso(i,k) ! - tem = 0.5 * cm * tem + tem = 0.5 * cfg%cm * tem factor = 1. + tem ptem = tem + pgcon ptem1= tem - pgcon @@ -1057,7 +1010,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & if(k > kb(i) .and. k < kmax(i)) then dz = zi(i,k) - zi(i,k-1) tem = 0.25 * (xlamue(i,k)+xlamue(i,k-1)) * dz - tem = cq * tem + tem = cfg%cq * tem factor = 1. + tem ecko(i,k,n) = ((1.-tem)*ecko(i,k-1,n)+tem* & (ctro(i,k,n)+ctro(i,k-1,n)))/factor @@ -1076,12 +1029,12 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & if(k > kb(i) .and. k < kmax(i)) then dz = zi(i,k) - zi(i,k-1) tem = 0.25 * (xlamue(i,k)+xlamue(i,k-1)) * dz - tem = cq * tem + tem = cfg%cq * tem factor = 1. + tem ecko(i,k,kk) = ((1. - tem) * ecko(i,k-1,kk) + tem * & (ctro(i,k,kk) + ctro(i,k-1,kk))) / factor ercko(i,k,kk) = ecko(i,k,kk) - chem_c(i,k,n) = escav * fscav(n) * ecko(i,k,kk) + chem_c(i,k,n) = cfg%escav * fscav(n) * ecko(i,k,kk) tem = chem_c(i,k,n) / (1. + c0t(i,k) * dz) chem_pw(i,k,n) = c0t(i,k) * dz * tem * eta(i,k-1) ecko(i,k,kk) = tem + ecko(i,k,kk) - chem_c(i,k,n) @@ -1119,7 +1072,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & do i=1,im if(cnvflg(i)) then tem = pfld(i,kbcon(i)) - pfld(i,kbcon1(i)) - if(tem > dthk) then + if(tem > cfg%dthk) then cnvflg(i) = .false. endif endif @@ -1162,7 +1115,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & if (hwrf_samfshal) then do i = 1, im if(cnvflg(i)) then - cinacr = cinacrmx + cinacr = cfg%cinacrmx if(cina(i) < cinacr) cnvflg(i) = .false. endif enddo @@ -1193,8 +1146,8 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & val2 = 1. tem = min(tem,val2) tem = 1. - tem - tem1= .5*(cinacrmx-cinacrmn) - cinacr = cinacrmx - tem * tem1 + tem1= .5*(cfg%cinacrmx-cfg%cinacrmn) + cinacr = cfg%cinacrmx - tem * tem1 if(cina(i) < cinacr) cnvflg(i) = .false. endif enddo @@ -1226,12 +1179,12 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & enddo enddo c -c turn off shallow convection if cloud depth is larger than cthk or less than cthkmn +c turn off shallow convection if cloud depth is larger than cfg%cthk or less than cfg%cthkmn c do i = 1, im if(cnvflg(i)) then tem = pfld(i,kbcon(i))-pfld(i,ktcon(i)) - if(tem > cthk .or. tem < cthkmn) then + if(tem > cfg%cthk .or. tem < cfg%cthkmn) then cnvflg(i) = .false. endif endif @@ -1270,7 +1223,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & + gamma * dbyo(i,k) / (hvap * (1. + gamma)) cj tem = 0.25 * (xlamue(i,k)+xlamue(i,k-1)) * dz - tem = cq * tem + tem = cfg%cq * tem factor = 1. + tem qcko(i,k) = ((1.-tem)*qcko(i,k-1)+tem* & (qo(i,k)+qo(i,k-1)))/factor @@ -1314,7 +1267,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & ! tem = ((uo(i,k)-uo(i,k-1))/dz)**2 tem = tem+((vo(i,k)-vo(i,k-1))/dz)**2 - wush(i,k) = csmf * sqrt(tem) + wush(i,k) = cfg%csmf * sqrt(tem) ! endif ! @@ -1443,7 +1396,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & + gamma * dbyo(i,k) / (hvap * (1. + gamma)) cj tem = 0.25 * (xlamue(i,k)+xlamue(i,k-1)) * dz - tem = cq * tem + tem = cfg%cq * tem factor = 1. + tem qcko(i,k) = ((1.-tem)*qcko(i,k-1)+tem* & (qo(i,k)+qo(i,k-1)))/factor @@ -1495,8 +1448,8 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & if (cnvflg(i)) then if(k > kbcon1(i) .and. k < ktcon(i)) then dz = zi(i,k) - zi(i,k-1) - tem = 0.25 * bb1 * (drag(i,k-1)+drag(i,k)) * dz - tem1 = 0.5 * bb2 * (buo(i,k-1)+buo(i,k)) + tem = 0.25 * cfg%bb1 * (drag(i,k-1)+drag(i,k)) * dz + tem1 = 0.5 * cfg%bb2 * (buo(i,k-1)+buo(i,k)) tem2 = wush(i,k) * sqrt(wu2(i,k-1)) tem2 = (tem1 - tem2) * dz ptem = (1. - tem) * wu2(i,k-1) @@ -1914,9 +1867,9 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & tfac = 1. + gdx(i) / 75000. dtconv(i) = tfac * dtconv(i) endif - dtconv(i) = max(dtconv(i),dtmin) + dtconv(i) = max(dtconv(i),cfg%dtmin) dtconv(i) = max(dtconv(i),dt2) - dtconv(i) = min(dtconv(i),dtmax) + dtconv(i) = min(dtconv(i),cfg%dtmax) endif enddo ! @@ -1980,7 +1933,8 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & call progsigma_calc(im,km,first_time_step,restart,flag_shallow, & flag_mid,del,tmfq,qmicro,dbyo1,zdqca,omega_u,zeta,hvap, & delt,qadv,kbcon1,ktcon,cnvflg,betascu,betamcu,betadcu, - & sigmind,sigminm,sigmins,sigmain,sigmaout,sigmab) + & cfg%sigmind,cfg%sigminm,cfg%sigmins, + & sigmain,sigmaout,sigmab) endif !> - From Han et al.'s (2017) \cite han_et_al_2017 equation 6, calculate cloud base mass flux as a function of the mean updraft velcoity. @@ -1990,10 +1944,10 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & if(cnvflg(i)) then k = kbcon(i) rho = po(i,k)*100. / (rd*to(i,k)) - if(progsigma .and. gdx(i) < dxcrtas)then + if(progsigma .and. gdx(i) < cfg%dxcrtas)then xmb(i) = advfac(i)*sigmab(i)*((-1.0*omegac(i))*gravinv) else - xmb(i) = advfac(i)*betaw*rho*wc(i) + xmb(i) = advfac(i)*cfg%betaw*rho*wc(i) endif endif enddo @@ -2330,7 +2284,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & ! evef = edt(i) * evfact ! if(islimsk(i) == 1) evef=edt(i) * evfactl ! if(islimsk(i) == 1) evef=.07 - qcond(i) = shevf * evef * (q1(i,k) - qeso(i,k)) + qcond(i) = cfg%shevf * evef * (q1(i,k) - qeso(i,k)) & / (1. + el2orc * qeso(i,k) / t1(i,k)**2) dp = 1000. * del(i,k) factor = dp / grav @@ -2434,7 +2388,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & ! if (k > kb(i) .and. k <= ktcon(i)) then if (k >= kbcon(i) .and. k <= ktcon(i)) then tem = dellal(i,k) * xmb(i) * dt2 - tem1 = max(0.0, min(1.0, (tcr-t1(i,k))*tcrf)) + tem1 = max(0.0, min(1.0, (cfg%tcr-t1(i,k))*cfg%tcrf)) if (qtr(i,k,2) > -999.0) then qtr(i,k,1) = qtr(i,k,1) + tem * tem1 ! ice qtr(i,k,2) = qtr(i,k,2) + tem *(1.0-tem1) ! water @@ -2499,7 +2453,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & if(progsigma)then tem2 = sigmab(i) else - tem2 = max(sigmagfm(i), betaw) + tem2 = max(sigmagfm(i), cfg%betaw) endif ptem = tem / (tem2 * tem1) qtr(i,k,ntk)=qtr(i,k,ntk)+0.5*tem2*ptem*ptem diff --git a/physics/samfshalcnv.meta b/physics/samfshalcnv.meta index 200e33707..bbd326808 100644 --- a/physics/samfshalcnv.meta +++ b/physics/samfshalcnv.meta @@ -2,6 +2,7 @@ name = samfshalcnv type = scheme dependencies = funcphys.f90,machine.F,samfaerosols.F,progsigma_calc.f90 + dependencies = module_samf.F90 ######################################################################## [ccpp-arg-table] @@ -41,6 +42,13 @@ [ccpp-arg-table] name = samfshalcnv_run type = scheme +[cfg] + standard_name = parameters_for_scale_aware_mass_flux_shallow_convection + long_name = parameters for SAMF shallow convection scheme + units = mixed + dimensions = () + type = ty_cfg_samf_shal + intent = in [im] standard_name = horizontal_loop_extent long_name = horizontal loop extent diff --git a/physics/satmedmfvdif.F b/physics/satmedmfvdif.F index 79f7bbea1..8327fd16a 100644 --- a/physics/satmedmfvdif.F +++ b/physics/satmedmfvdif.F @@ -61,7 +61,7 @@ end subroutine satmedmfvdif_init !! (mfscu.f). !! \section detail_satmedmfvidf GFS satmedmfvdif Detailed Algorithm !> @{ - subroutine satmedmfvdif_run(im,km,ntrac,ntcw,ntiw,ntke, & + subroutine satmedmfvdif_run(cfg,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, & & psk,rbsoil,zorl,u10m,v10m,fm,fh, & @@ -75,10 +75,12 @@ subroutine satmedmfvdif_run(im,km,ntrac,ntcw,ntiw,ntke, & ! use machine , only : kind_phys use funcphys , only : fpvs + use module_satmedmf, only: ty_cfg_satmedmf ! implicit none ! !---------------------------------------------------------------------- + type(ty_cfg_satmedmf), intent(in) :: cfg integer, intent(in) :: im, km, ntrac, ntcw, ntiw, ntke integer, intent(in) :: kinver(:) integer, intent(out) :: kpbl(:) @@ -183,58 +185,31 @@ subroutine satmedmfvdif_run(im,km,ntrac,ntcw,ntiw,ntke, & ! ! pcnvflg: true for unstable pbl ! - real(kind=kind_phys) aphi16, aphi5, - & wfac, cfac, - & gamcrt, gamcrq, sfcfrac, + real(kind=kind_phys) & conq, cont, conw, - & dsdz2, dsdzt, dkmax, + & dsdz2, dsdzt, & dsig, dt2, dtodsd, & dtodsu, g, factor, dz, - & gocp, gravi, zol1, zolcru, + & gocp, gravi, zol1, & buop, shrp, dtn, cdtn, - & prnum, prmax, prmin, prtke, - & prscu, dw2, dw2min, zk, - & elmfac, elefac, dspmax, - & alp, clwt, cql, - & f0, robn, crbmin, crbmax, + & prnum, + & dw2, zk, + & alp, clwt, + & robn, & es, qs, value, onemrh, & cfh, gamma, elocp, el2orc, & epsi, beta, chx, cqx, - & rdt, rdz, qmin, qlmin, - & ri, rimin, - & rbcr, rbint, tdzmin, - & rlmn, rlmx, elmx, + & rdt, rdz, + & ri, + & rbint, & ttend, utend, vtend, qtend, - & zfac, zfmin, vk, spdk2, - & tkmin, xkzinv, dspfac, xkgdx, + & zfac, spdk2, & zlup, zldn, bsum, & tem, tem1, tem2, & ptem, ptem0, ptem1, ptem2 ! - real(kind=kind_phys) ck0, ck1, ch0, ch1, ce0, rchck -! - real(kind=kind_phys) qlcr, zstblmax -! - real(kind=kind_phys) h1 integer :: idtend !! - parameter(wfac=7.0,cfac=4.5) - parameter(gamcrt=3.,gamcrq=0.,sfcfrac=0.1) - parameter(vk=0.4,rimin=-100.) - parameter(rbcr=0.25,zolcru=-0.02,tdzmin=1.e-3) - parameter(rlmn=30.,rlmx=500.,elmx=500.) - parameter(prmin=0.25,prmax=4.0,prtke=1.0,prscu=0.67) - parameter(f0=1.e-4,crbmin=0.15,crbmax=0.35) - parameter(tkmin=1.e-9,dspfac=0.5,dspmax=10.0) - parameter(qmin=1.e-8,qlmin=1.e-12,zfmin=1.e-8) - parameter(aphi5=5.,aphi16=16.) - parameter(elmfac=1.0,elefac=1.0,cql=100.) - parameter(dw2min=1.e-4,dkmax=1000.,xkgdx=25000.) - parameter(qlcr=3.5e-5,zstblmax=2500.,xkzinv=0.15) - parameter(h1=0.33333333) - parameter(ck0=0.4,ck1=0.15,ch0=0.4,ch1=0.15,ce0=0.4) - parameter(rchck=1.5,cdtn=25.) - gravi=1.0/grav g=grav gocp=g/cp @@ -272,8 +247,8 @@ subroutine satmedmfvdif_run(im,km,ntrac,ntcw,ntiw,ntke, & xmfd(i,k) = 0. buou(i,k) = 0. buod(i,k) = 0. - ckz(i,k) = ck1 - chz(i,k) = ch1 + ckz(i,k) = cfg%ck1 + chz(i,k) = cfg%ch1 enddo enddo do i=1,im @@ -292,7 +267,7 @@ subroutine satmedmfvdif_run(im,km,ntrac,ntcw,ntiw,ntke, & !! from tracer (\p tke and \p tkeh) do k=1,km do i=1,im - tke(i,k) = max(q1(i,k,ntke), tkmin) + tke(i,k) = max(q1(i,k,ntke), cfg%tkmin) enddo enddo do k=1,km1 @@ -313,18 +288,18 @@ subroutine satmedmfvdif_run(im,km,ntrac,ntcw,ntiw,ntke, & !! - 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.) +!!\n xkzm_hx = 0.01 + (xkzm_h - 0.01)/(cfg%xkgdx-5.) * (gdx-5.) +!!\n xkzm_mx = 0.01 + (xkzm_h - 0.01)/(cfg%xkgdx-5.) * (gdx-5.) do i=1,im kx1(i) = 1 tx1(i) = 1.0 / prsi(i,1) tx2(i) = tx1(i) - if(gdx(i) >= xkgdx) then + if(gdx(i) >= cfg%xkgdx) then xkzm_hx(i) = xkzm_h xkzm_mx(i) = xkzm_m else - tem = 1. / (xkgdx - 5.) + tem = 1. / (cfg%xkgdx - 5.) tem1 = (xkzm_h - 0.01) * tem tem2 = (xkzm_m - 0.01) * tem ptem = gdx(i) - 5. @@ -388,37 +363,37 @@ subroutine satmedmfvdif_run(im,km,ntrac,ntcw,ntiw,ntke, & pix(i,k) = psk(i) / prslk(i,k) theta(i,k) = t1(i,k) * pix(i,k) if(ntiw > 0) then - tem = max(q1(i,k,ntcw),qlmin) - tem1 = max(q1(i,k,ntiw),qlmin) + tem = max(q1(i,k,ntcw),cfg%qlmin) + tem1 = max(q1(i,k,ntiw),cfg%qlmin) qlx(i,k) = tem + tem1 ptem = hvap*tem + (hvap+hfus)*tem1 slx(i,k) = cp * t1(i,k) + phil(i,k) - ptem else - qlx(i,k) = max(q1(i,k,ntcw),qlmin) + qlx(i,k) = max(q1(i,k,ntcw),cfg%qlmin) slx(i,k) = cp * t1(i,k) + phil(i,k) - hvap*qlx(i,k) endif - tem2 = 1.+fv*max(q1(i,k,1),qmin)-qlx(i,k) + tem2 = 1.+fv*max(q1(i,k,1),cfg%qmin)-qlx(i,k) thvx(i,k) = theta(i,k) * tem2 tvx(i,k) = t1(i,k) * tem2 - qtx(i,k) = max(q1(i,k,1),qmin)+qlx(i,k) + qtx(i,k) = max(q1(i,k,1),cfg%qmin)+qlx(i,k) thlx(i,k) = theta(i,k) - pix(i,k)*elocp*qlx(i,k) thlvx(i,k) = thlx(i,k) * (1. + fv * qtx(i,k)) svx(i,k) = cp * tvx(i,k) - ptem1 = elocp * pix(i,k) * max(q1(i,k,1),qmin) + ptem1 = elocp * pix(i,k) * max(q1(i,k,1),cfg%qmin) thetae(i,k)= theta(i,k) + ptem1 gotvx(i,k) = g / tvx(i,k) enddo enddo ! !> - The background vertical diffusivities in the inversion layers are limited -!! to be less than or equal to xkzinv +!! to be less than or equal to cfg%xkzinv ! do k = 1,km1 do i=1,im tem1 = (tvx(i,k+1)-tvx(i,k)) * rdzt(i,k) if(tem1 > 1.e-5) then - xkzo(i,k) = min(xkzo(i,k),xkzinv) - xkzmo(i,k) = min(xkzmo(i,k),xkzinv) + xkzo(i,k) = min(xkzo(i,k),cfg%xkzinv) + xkzmo(i,k) = min(xkzmo(i,k),cfg%xkzinv) endif enddo enddo @@ -430,8 +405,8 @@ subroutine satmedmfvdif_run(im,km,ntrac,ntcw,ntiw,ntke, & plyr(i,k) = 0.01 * prsl(i,k) ! pa to mb (hpa) ! --- ... compute relative humidity es = 0.01 * fpvs(t1(i,k)) ! fpvs in pa - qs = max(qmin, eps * es / (plyr(i,k) + epsm1*es)) - rhly(i,k) = max(0.0, min(1.0, max(qmin, q1(i,k,1))/qs)) + qs = max(cfg%qmin, eps * es / (plyr(i,k) + epsm1*es)) + rhly(i,k) = max(0.0, min(1.0, max(cfg%qmin, q1(i,k,1))/qs)) qstl(i,k) = qs enddo enddo @@ -443,7 +418,7 @@ subroutine satmedmfvdif_run(im,km,ntrac,ntcw,ntiw,ntke, & if (qlx(i,k) > clwt) then onemrh= max(1.e-10, 1.0-rhly(i,k)) tem1 = min(max((onemrh*qstl(i,k))**0.49,0.0001),1.0) - tem1 = cql / tem1 + tem1 = cfg%cql / tem1 value = max(min( tem1*qlx(i,k), 50.0), 0.0) tem2 = sqrt(sqrt(rhly(i,k))) cfly(i,k) = min(max(tem2*(1.0-exp(-value)), 0.0), 1.0) @@ -514,15 +489,15 @@ subroutine satmedmfvdif_run(im,km,ntrac,ntcw,ntiw,ntke, & if(pblflg(i)) then ! thermal(i) = thvx(i,1) thermal(i) = thlvx(i,1) - crb(i) = rbcr + crb(i) = cfg%rbcr else - thermal(i) = tsea(i)*(1.+fv*max(q1(i,1,1),qmin)) + thermal(i) = tsea(i)*(1.+fv*max(q1(i,1,1),cfg%qmin)) tem = sqrt(u10m(i)**2+v10m(i)**2) tem = max(tem, 1.) - robn = tem / (f0 * z0(i)) + robn = tem / (cfg%f0 * z0(i)) tem1 = 1.e-7 * robn crb(i) = 0.16 * (tem1 ** (-0.18)) - crb(i) = max(min(crb(i), crbmax), crbmin) + crb(i) = max(min(crb(i), cfg%crbmax), cfg%crbmin) endif enddo ! @@ -544,7 +519,7 @@ subroutine satmedmfvdif_run(im,km,ntrac,ntcw,ntiw,ntke, & ! bf(i,k) = gotvx(i,k)*(thvx(i,k+1)-thvx(i,k))*rdz 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 + shr2(i,k) = max(dw2,cfg%dw2min)*rdz*rdz enddo enddo ! @@ -609,11 +584,11 @@ subroutine satmedmfvdif_run(im,km,ntrac,ntcw,ntiw,ntke, & !! 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) + zol(i) = max(rbsoil(i)*fm(i)*fm(i)/fh(i),cfg%rimin) if(sfcflg(i)) then - zol(i) = min(zol(i),-zfmin) + zol(i) = min(zol(i),-cfg%zfmin) else - zol(i) = max(zol(i),zfmin) + zol(i) = max(zol(i),cfg%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: @@ -627,13 +602,13 @@ subroutine satmedmfvdif_run(im,km,ntrac,ntcw,ntiw,ntke, & !!\f[ !! \phi_m=\phi_t=(1+5\frac{0.1h}{L}) !!\f] - zol1 = zol(i)*sfcfrac*hpbl(i)/zl(i,1) + zol1 = zol(i)*cfg%sfcfrac*hpbl(i)/zl(i,1) if(sfcflg(i)) then - tem = 1.0 / (1. - aphi16*zol1) + tem = 1.0 / (1. - cfg%aphi16*zol1) phih(i) = sqrt(tem) phim(i) = sqrt(phih(i)) else - phim(i) = 1. + aphi5*zol1 + phim(i) = 1. + cfg%aphi5*zol1 phih(i) = phim(i) endif enddo @@ -647,7 +622,7 @@ subroutine satmedmfvdif_run(im,km,ntrac,ntcw,ntiw,ntke, & !! 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), +!! of the surface layer height to the PBL height (specified as cfg%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[ @@ -655,14 +630,15 @@ subroutine satmedmfvdif_run(im,km,ntrac,ntcw,ntiw,ntke, & !!\f] do i=1,im if(pblflg(i)) then - if(zol(i) < zolcru) then + if(zol(i) < cfg%zolcru) then pcnvflg(i) = .true. endif wst3(i) = gotvx(i,1)*sflux(i)*hpbl(i) - wstar(i)= wst3(i)**h1 + wstar(i)= wst3(i)**cfg%h1 ust3(i) = ustar(i)**3. - wscale(i)=(ust3(i)+wfac*vk*wst3(i)*sfcfrac)**h1 - ptem = ustar(i)/aphi5 + wscale(i)=(ust3(i)+cfg%wfac*cfg%vk*wst3(i)* + & cfg%sfcfrac)**cfg%h1 + ptem = ustar(i)/cfg%aphi5 wscale(i) = max(wscale(i),ptem) endif enddo @@ -675,7 +651,7 @@ subroutine satmedmfvdif_run(im,km,ntrac,ntcw,ntiw,ntke, & hgamq(i) = evap(i)/wscale(i) vpert(i) = hgamt(i) + hgamq(i)*fv*theta(i,1) vpert(i) = max(vpert(i),0.) - tem = min(cfac*vpert(i),gamcrt) + tem = min(cfg%cfac*vpert(i),cfg%gamcrt) thermal(i)= thermal(i) + tem endif enddo @@ -736,7 +712,7 @@ subroutine satmedmfvdif_run(im,km,ntrac,ntcw,ntiw,ntke, & enddo do k = 1, km1 do i=1,im - if(flg(i).and.zl(i,k) >= zstblmax) then + if(flg(i).and.zl(i,k) >= cfg%zstblmax) then lcld(i)=k flg(i)=.false. endif @@ -748,7 +724,7 @@ subroutine satmedmfvdif_run(im,km,ntrac,ntcw,ntiw,ntke, & do k = kmscu,1,-1 do i = 1, im if(flg(i) .and. k <= lcld(i)) then - if(qlx(i,k) >= qlcr) then + if(qlx(i,k) >= cfg%qlcr) then kcld(i)=k flg(i)=.false. endif @@ -769,7 +745,7 @@ subroutine satmedmfvdif_run(im,km,ntrac,ntcw,ntiw,ntke, & do k = kmscu,1,-1 do i = 1, im if(flg(i) .and. k <= kcld(i)) then - if(qlx(i,k) >= qlcr) then + if(qlx(i,k) >= cfg%qlcr) then if(radx(i,k) < radmin(i)) then radmin(i)=radx(i,k) krad(i)=k @@ -837,29 +813,29 @@ subroutine satmedmfvdif_run(im,km,ntrac,ntcw,ntiw,ntke, & !! is the PBL height) non-dimensional gradient functions for heat and momentum, !! respectively. !! - We assume that \f$P_r=0.67\f$ in the stratocumulus and the unstable -!! layers above PBL (prscu). +!! layers above PBL (cfg%prscu). !! - For the stable layers above PBL, \f$P_r\f$ is given as \f$P_r=1+2.1R_i\f$ !! (prnum) (Kennedy and Shapiro (1980)\cite kennedy_and_shapiro_1980 ) do k = 1, kmpbl do i = 1, im if(k < kpbl(i)) then tem = phih(i)/phim(i) - ptem = -3.*(max(zi(i,k+1)-sfcfrac*hpbl(i),0.))**2. + ptem = -3.*(max(zi(i,k+1)-cfg%sfcfrac*hpbl(i),0.))**2. & /hpbl(i)**2. if(pcnvflg(i)) then prn(i,k) = 1. + (tem-1.)*exp(ptem) else prn(i,k) = tem endif - prn(i,k) = min(prn(i,k),prmax) - prn(i,k) = max(prn(i,k),prmin) -! - ckz(i,k) = ck1 + (ck0-ck1)*exp(ptem) - ckz(i,k) = min(ckz(i,k),ck0) - ckz(i,k) = max(ckz(i,k),ck1) - chz(i,k) = ch1 + (ch0-ch1)*exp(ptem) - chz(i,k) = min(chz(i,k),ch0) - chz(i,k) = max(chz(i,k),ch1) + prn(i,k) = min(prn(i,k),cfg%prmax) + prn(i,k) = max(prn(i,k),cfg%prmin) +! + ckz(i,k) = cfg%ck1 + (cfg%ck0-cfg%ck1)*exp(ptem) + ckz(i,k) = min(ckz(i,k),cfg%ck0) + ckz(i,k) = max(ckz(i,k),cfg%ck1) + chz(i,k) = cfg%ch1 + (cfg%ch0-cfg%ch1)*exp(ptem) + chz(i,k) = min(chz(i,k),cfg%ch0) + chz(i,k) = max(chz(i,k),cfg%ch1) endif enddo enddo @@ -882,9 +858,9 @@ subroutine satmedmfvdif_run(im,km,ntrac,ntcw,ntiw,ntke, & zlup = zlup + dz if(bsum >= tke(i,k)) then if(ptem >= 0.) then - tem2 = max(ptem, zfmin) + tem2 = max(ptem, cfg%zfmin) else - tem2 = min(ptem, -zfmin) + tem2 = min(ptem, -cfg%zfmin) endif ptem1 = (bsum - tke(i,k)) / tem2 zlup = zlup - ptem1 * dz @@ -900,7 +876,7 @@ subroutine satmedmfvdif_run(im,km,ntrac,ntcw,ntiw,ntke, & if(mlenflg) then if(n == 1) then dz = zl(i,1) - tem1 = tsea(i)*(1.+fv*max(q1(i,1,1),qmin)) + tem1 = tsea(i)*(1.+fv*max(q1(i,1,1),cfg%qmin)) else dz = zl(i,n) - zl(i,n-1) tem1 = thvx(i,n-1) @@ -912,9 +888,9 @@ subroutine satmedmfvdif_run(im,km,ntrac,ntcw,ntiw,ntke, & zldn = zldn + dz if(bsum >= tke(i,k)) then if(ptem >= 0.) then - tem2 = max(ptem, zfmin) + tem2 = max(ptem, cfg%zfmin) else - tem2 = min(ptem, -zfmin) + tem2 = min(ptem, -cfg%zfmin) endif ptem1 = (bsum - tke(i,k)) / tem2 zldn = zldn - ptem1 * dz @@ -925,7 +901,7 @@ subroutine satmedmfvdif_run(im,km,ntrac,ntcw,ntiw,ntke, & enddo ! tem = 0.5 * (zi(i,k+1)-zi(i,k)) - tem1 = min(tem, rlmn) + tem1 = min(tem, cfg%rlmn) !> - 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[ @@ -939,14 +915,14 @@ subroutine satmedmfvdif_run(im,km,ntrac,ntcw,ntiw,ntke, & !! 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) = cfg%elmfac * ptem2 rlam(i,k) = max(rlam(i,k), tem1) - rlam(i,k) = min(rlam(i,k), rlmx) + rlam(i,k) = min(rlam(i,k), cfg%rlmx) ! ptem2 = sqrt(zlup*zldn) - ele(i,k) = elefac * ptem2 + ele(i,k) = cfg%elefac * ptem2 ele(i,k) = max(ele(i,k), tem1) - ele(i,k) = min(ele(i,k), elmx) + ele(i,k) = min(ele(i,k), cfg%elmx) ! enddo enddo @@ -954,7 +930,7 @@ subroutine satmedmfvdif_run(im,km,ntrac,ntcw,ntiw,ntke, & !! 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) + tem = cfg%vk * zl(i,k) if (zol(i) < 0.) then ptem = 1. - 100. * zol(i) ptem1 = ptem**0.2 @@ -994,14 +970,14 @@ subroutine satmedmfvdif_run(im,km,ntrac,ntcw,ntiw,ntke, & dku(i,k) = dkt(i,k) * prn(i,k) endif else - ri = max(bf(i,k)/shr2(i,k),rimin) + ri = max(bf(i,k)/shr2(i,k),cfg%rimin) if(ri < 0.) then ! unstable regime - dku(i,k) = ck1 * tem + dku(i,k) = cfg%ck1 * tem dkt(i,k) = rchck * dku(i,k) else ! stable regime - dkt(i,k) = ch1 * tem + dkt(i,k) = cfg%ch1 * tem prnum = 1.0 + 2.1*ri - prnum = min(prnum,prmax) + prnum = min(prnum,cfg%prmax) dku(i,k) = dkt(i,k) * prnum endif endif @@ -1009,19 +985,19 @@ subroutine satmedmfvdif_run(im,km,ntrac,ntcw,ntiw,ntke, & if(scuflg(i)) then if(k >= mrad(i) .and. k < krad(i)) then tem1 = ckz(i,k) * tem - ptem1 = tem1 / prscu + ptem1 = tem1 / cfg%prscu dku(i,k) = max(dku(i,k), tem1) dkt(i,k) = max(dkt(i,k), ptem1) endif endif ! - dkq(i,k) = prtke * dkt(i,k) + dkq(i,k) = cfg%prtke * dkt(i,k) ! - dkt(i,k) = min(dkt(i,k),dkmax) + dkt(i,k) = min(dkt(i,k),cfg%dkmax) dkt(i,k) = max(dkt(i,k),xkzo(i,k)) - dkq(i,k) = min(dkq(i,k),dkmax) + dkq(i,k) = min(dkq(i,k),cfg%dkmax) dkq(i,k) = max(dkq(i,k),xkzo(i,k)) - dku(i,k) = min(dku(i,k),dkmax) + dku(i,k) = min(dku(i,k),cfg%dkmax) dku(i,k) = max(dku(i,k),xkzmo(i,k)) ! enddo @@ -1031,7 +1007,7 @@ subroutine satmedmfvdif_run(im,km,ntrac,ntcw,ntiw,ntke, & if(scuflg(i)) then k = krad(i) tem = bf(i,k) / gotvx(i,k) - tem1 = max(tem, tdzmin) + tem1 = max(tem, cfg%tdzmin) ptem = radj(i) / tem1 dkt(i,k) = dkt(i,k) + ptem dku(i,k) = dku(i,k) + ptem @@ -1092,7 +1068,7 @@ subroutine satmedmfvdif_run(im,km,ntrac,ntcw,ntiw,ntke, & ptem2 = ptem2 + ptem ! ! tem2 = stress(i)*spd1(i)/zl(i,1) - tem2 = stress(i)*ustar(i)*phim(i)/(vk*zl(i,1)) + tem2 = stress(i)*ustar(i)*phim(i)/(cfg%vk*zl(i,1)) shrp = 0.5 * (tem1 + ptem1 + ptem2 + tem2) else tem1 = -dkt(i,k-1) * bf(i,k-1) @@ -1167,7 +1143,7 @@ subroutine satmedmfvdif_run(im,km,ntrac,ntcw,ntiw,ntke, & ! do k = 1,km1 do i=1,im - rle(i,k) = ce0 / ele(i,k) + rle(i,k) = cfg%ce0 / ele(i,k) enddo enddo kk = max(nint(dt2/cdtn), 1) @@ -1180,7 +1156,7 @@ subroutine satmedmfvdif_run(im,km,ntrac,ntcw,ntiw,ntke, & tem1 = prod(i,k) + tke(i,k) / dtn diss(i,k)=max(min(diss(i,k), tem1), 0.) tke(i,k) = tke(i,k) + dtn * (prod(i,k)-diss(i,k)) - tke(i,k) = max(tke(i,k), tkmin) + tke(i,k) = max(tke(i,k), cfg%tkmin) enddo enddo enddo @@ -1279,7 +1255,7 @@ subroutine satmedmfvdif_run(im,km,ntrac,ntcw,ntiw,ntke, & ! do k = 1,km do i = 1,im -! f1(i,k) = max(f1(i,k), tkmin) +! f1(i,k) = max(f1(i,k), cfg%tkmin) qtend = (f1(i,k)-q1(i,k,ntke))*rdt rtg(i,k,ntke) = rtg(i,k,ntke)+qtend enddo @@ -1438,10 +1414,10 @@ subroutine satmedmfvdif_run(im,km,ntrac,ntcw,ntiw,ntke, & if(dspheat) then do k = 1,km1 do i = 1,im -! tem = min(diss(i,k), dspmax) +! tem = min(diss(i,k), cfg%dspmax) ! ttend = tem / cp ttend = diss(i,k) / cp - tdt(i,k) = tdt(i,k) + dspfac * ttend + tdt(i,k) = tdt(i,k) + cfg%dspfac * ttend enddo enddo endif diff --git a/physics/satmedmfvdif.meta b/physics/satmedmfvdif.meta index 3609ed50f..b5014cdce 100644 --- a/physics/satmedmfvdif.meta +++ b/physics/satmedmfvdif.meta @@ -2,6 +2,7 @@ name = satmedmfvdif type = scheme dependencies = funcphys.f90,machine.F,mfpblt.f,mfscu.f,tridi.f + dependencies = module_satmedmf.F90 ######################################################################## [ccpp-arg-table] @@ -49,6 +50,13 @@ [ccpp-arg-table] name = satmedmfvdif_run type = scheme +[cfg] + standard_name = parameters_for_scale_aware_TKE_moist_EDMF_PBL_2018 + long_name = parameters for scale-aware TKE moist EDMF PBL scheme + units = mixed + dimensions = () + type = ty_cfg_satmedmf + intent = in [im] standard_name = horizontal_loop_extent long_name = horizontal loop extent diff --git a/physics/satmedmfvdifq.F b/physics/satmedmfvdifq.F index 7b54b6d12..3020efa23 100644 --- a/physics/satmedmfvdifq.F +++ b/physics/satmedmfvdifq.F @@ -73,7 +73,7 @@ end subroutine satmedmfvdifq_init !! -# A mass-flux approach is also used to represent the stratocumulus-top-induced turbulence !! (mfscuq.f). !! \section detail_satmedmfvidfq GFS satmedmfvdifq Detailed Algorithm - subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & + subroutine satmedmfvdifq_run(cfg,im,km,ntrac,ntcw,ntrw, & & ntiw,ntke,grav,rd,cp,rv,hvap,hfus,fv,eps,epsm1, & & dv,du,tdt,rtg,u1,v1,t1,q1,swh,hlw,xmu, & & garea,zvfun,sigmaf, & @@ -89,10 +89,12 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & ! use machine , only : kind_phys use funcphys , only : fpvs + use module_satmedmf, only: ty_cfg_satmedmfq ! implicit none ! !---------------------------------------------------------------------- + type(ty_cfg_satmedmfq), intent(in) :: cfg integer, intent(in) :: im, km, ntrac, ntcw, ntrw, ntiw, & & ntke, ntqv integer, intent(in) :: sfc_rlm @@ -143,7 +145,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & !*** !*** local variables !*** - integer i,is,k,n,ndt,km1,kmpbl,kmscu,ntrac1,idtend + integer i,is,k,n,km1,kmpbl,kmscu,ntrac1,idtend integer kps,kbx,kmx integer lcld(im),kcld(im),krad(im),mrad(im) integer kx1(im), kb1(im), kpblx(im) @@ -215,79 +217,31 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & ! ! pcnvflg: true for unstable pbl ! - real(kind=kind_phys) aphi16, aphi5, - & wfac, cfac, - & gamcrt, gamcrq, sfcfrac, + real(kind=kind_phys) ! & conq, cont, conw, - & dsdz2, dsdzt, dkmax, + & dsdz2, dsdzt, & dsig, dt2, dtodsd, & dtodsu, g, factor, dz, - & gocp, gravi, zol1, zolcru, + & gocp, gravi, zol1, & buop, shrp, dtn, - & prnum, prmax, prmin, prtke, - & prscu, pr0, ri, - & dw2, dw2min, zk, - & elmfac, elefac, dspmax, - & alp, clwt, cql, - & f0, robn, crbmin, crbmax, + & prnum, + & ri, + & dw2, zk, + & alp, clwt, + & robn, & es, qs, value, onemrh, & cfh, gamma, elocp, el2orc, & epsi, beta, chx, cqx, - & rdt, rdz, qmin, qlmin, - & rimin, rbcr, rbint, tdzmin, - & rlmn, rlmn0, rlmn1, rlmn2, + & rdt, rdz, + & rbint, & ttend, utend, vtend, qtend, - & zfac, zfmin, vk, spdk2, - & tkmin, tkbmx, xkgdx, - & xkinv1, xkinv2, - & zlup, zldn, cs0, csmf, + & zfac, spdk2, + & zlup, zldn, & tem, tem1, tem2, tem3, & ptem, ptem0, ptem1, ptem2 ! - real(kind=kind_phys) slfac -! - real(kind=kind_phys) vegflo, vegfup, z0lo, z0up, vc0, zc0 -! - real(kind=kind_phys) ck0, ck1, ch0, ch1, ce0, rchck -! - real(kind=kind_phys) qlcr, zstblmax, hcrinv -! - real(kind=kind_phys) h1 - - real(kind=kind_phys) bfac, mffac + real(kind=kind_phys) mffac !! - parameter(bfac=100.) - parameter(wfac=7.0,cfac=4.5) - parameter(gamcrt=3.,gamcrq=0.,sfcfrac=0.1) - parameter(vk=0.4,rimin=-100.,slfac=0.1) - parameter(rbcr=0.25,zolcru=-0.02,tdzmin=1.e-3) - parameter(rlmn=30.,rlmn0=5.,rlmn1=5.,rlmn2=10.) - parameter(prmin=0.25,prmax=4.0) - parameter(pr0=1.0,prtke=1.0,prscu=0.67) - parameter(f0=1.e-4,crbmin=0.15,crbmax=0.35) - parameter(tkmin=1.e-9,tkbmx=0.2,dspmax=10.0) - parameter(qmin=1.e-8,qlmin=1.e-12,zfmin=1.e-8) - parameter(aphi5=5.,aphi16=16.) - parameter(elmfac=1.0,elefac=1.0,cql=100.) - parameter(dw2min=1.e-4,dkmax=1000.,xkgdx=1000.) - parameter(qlcr=3.5e-5,zstblmax=2500.) - parameter(xkinv1=0.4,xkinv2=0.3) - parameter(h1=0.33333333,hcrinv=250.) - parameter(vegflo=0.1,vegfup=1.0,z0lo=0.1,z0up=1.0) - parameter(vc0=1.0,zc0=1.0) - parameter(ck1=0.15,ch1=0.15) - parameter(cs0=0.4,csmf=0.5) - parameter(rchck=1.5,ndt=20) - - if (tc_pbl == 0) then - ck0 = 0.4 - ch0 = 0.4 - ce0 = 0.4 - else if (tc_pbl == 1) then - ck0 = 0.55 - ch0 = 0.55 - ce0 = 0.12 - endif gravi = 1.0 / grav g = grav gocp = g / cp @@ -325,9 +279,9 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & buou(i,k) = 0. buod(i,k) = 0. wush(i,k) = 0. - ckz(i,k) = ck1 - chz(i,k) = ch1 - rlmnz(i,k) = rlmn0 + ckz(i,k) = cfg%ck1 + chz(i,k) = cfg%ch1 + rlmnz(i,k) = cfg%rlmn0 enddo enddo do i=1,im @@ -346,7 +300,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & !! from tracer (\p tke and \p tkeh) do k=1,km do i=1,im - tke(i,k) = max(q1(i,k,ntke), tkmin) + tke(i,k) = max(q1(i,k,ntke), cfg%tkmin) enddo enddo do k=1,km1 @@ -358,7 +312,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & do k = 1,km1 do i=1,im rdzt(i,k) = 1.0 / (zl(i,k+1) - zl(i,k)) - prn(i,k) = pr0 + prn(i,k) = cfg%pr0 enddo enddo ! @@ -368,20 +322,20 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & !> - Compute background vertical diffusivities for scalars and momentum (xkzo and xkzmo) -!> - set background diffusivities with xkzm_h & xkzm_m for gdx >= xkgdx and -!! as a function of horizontal grid size for gdx < xkgdx -!! \n xkzm_hx = xkzm_h * (gdx / xkgdx) -!! \n xkzm_mx = xkzm_m * (gdx / xkgdx) +!> - set background diffusivities with xkzm_h & xkzm_m for gdx >= cfg%xkgdx and +!! as a function of horizontal grid size for gdx < cfg%xkgdx +!! \n xkzm_hx = xkzm_h * (gdx / cfg%xkgdx) +!! \n xkzm_mx = xkzm_m * (gdx / cfg%xkgdx) ! do i=1,im kx1(i) = 1 tx1(i) = 1.0 / prsi(i,1) tx2(i) = tx1(i) - if(gdx(i) >= xkgdx) then + if(gdx(i) >= cfg%xkgdx) then xkzm_hx(i) = xkzm_h xkzm_mx(i) = xkzm_m else - tem = gdx(i) / xkgdx + tem = gdx(i) / cfg%xkgdx xkzm_hx(i) = xkzm_h * tem xkzm_mx(i) = xkzm_m * tem endif @@ -396,8 +350,8 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & tem1 = 1.0 - ptem tem2 = tem1 * tem1 * 2.5 tem2 = min(1.0, exp(-tem2)) - rlmnz(i,k)= rlmn * tem2 - rlmnz(i,k)= max(rlmnz(i,k), rlmn0) + rlmnz(i,k)= cfg%rlmn * tem2 + rlmnz(i,k)= max(rlmnz(i,k), cfg%rlmn0) ! vertical background diffusivity tem2 = tem1 * tem1 * 10.0 tem2 = min(1.0, exp(-tem2)) @@ -420,7 +374,8 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & !> - Some output variables and logical flags are initialized do i = 1,im z0(i) = 0.01 * zorl(i) - rho_a(i) = prsl(i,1)/(rd*t1(i,1)*(1.+fv*max(q1(i,1,1),qmin))) + rho_a(i) = prsl(i,1)/(rd*t1(i,1)* + & (1.+fv*max(q1(i,1,1),cfg%qmin))) dusfc(i) = 0. dvsfc(i) = 0. dtsfc(i) = 0. @@ -448,12 +403,12 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & !! roughness length ! do i = 1,im - tem = (sigmaf(i) - vegflo) / (vegfup - vegflo) + tem = (sigmaf(i) - cfg%vegflo) / (cfg%vegfup - cfg%vegflo) tem = min(max(tem, 0.), 1.) tem1 = sqrt(tem) - ptem = (z0(i) - z0lo) / (z0up - z0lo) + ptem = (z0(i) - cfg%z0lo) / (cfg%z0up - cfg%z0lo) ptem = min(max(ptem, 0.), 1.) - vez0fun(i) = (1. + vc0 * tem1) * (1. + zc0 * ptem) + vez0fun(i) = (1. + cfg%vc0 * tem1) * (1. + cfg%zc0 * ptem) enddo ! !> - Compute \f$\theta\f$(theta), and \f$q_l\f$(qlx), \f$\theta_e\f$(thetae), @@ -463,23 +418,23 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & pix(i,k) = psk(i) / prslk(i,k) theta(i,k) = t1(i,k) * pix(i,k) if(ntiw > 0) then - tem = max(q1(i,k,ntcw),qlmin) - tem1 = max(q1(i,k,ntiw),qlmin) + tem = max(q1(i,k,ntcw),cfg%qlmin) + tem1 = max(q1(i,k,ntiw),cfg%qlmin) qlx(i,k) = tem + tem1 ptem = hvap*tem + (hvap+hfus)*tem1 slx(i,k) = cp * t1(i,k) + phil(i,k) - ptem else - qlx(i,k) = max(q1(i,k,ntcw),qlmin) + qlx(i,k) = max(q1(i,k,ntcw),cfg%qlmin) slx(i,k) = cp * t1(i,k) + phil(i,k) - hvap*qlx(i,k) endif - tem2 = 1.+fv*max(q1(i,k,1),qmin)-qlx(i,k) + tem2 = 1.+fv*max(q1(i,k,1),cfg%qmin)-qlx(i,k) thvx(i,k) = theta(i,k) * tem2 tvx(i,k) = t1(i,k) * tem2 - qtx(i,k) = max(q1(i,k,1),qmin)+qlx(i,k) + qtx(i,k) = max(q1(i,k,1),cfg%qmin)+qlx(i,k) thlx(i,k) = theta(i,k) - pix(i,k)*elocp*qlx(i,k) thlvx(i,k) = thlx(i,k) * (1. + fv * qtx(i,k)) svx(i,k) = cp * tvx(i,k) - ptem1 = elocp * pix(i,k) * max(q1(i,k,1),qmin) + ptem1 = elocp * pix(i,k) * max(q1(i,k,1),cfg%qmin) thetae(i,k)= theta(i,k) + ptem1 ! gotvx(i,k) = g / tvx(i,k) gotvx(i,k) = g / thvx(i,k) @@ -493,8 +448,8 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & plyr(i,k) = 0.01 * prsl(i,k) ! pa to mb (hpa) ! --- ... compute relative humidity es = 0.01 * fpvs(t1(i,k)) ! fpvs in pa - qs = max(qmin, eps * es / (plyr(i,k) + epsm1*es)) - rhly(i,k) = max(0.0, min(1.0, max(qmin, q1(i,k,1))/qs)) + qs = max(cfg%qmin, eps * es / (plyr(i,k) + epsm1*es)) + rhly(i,k) = max(0.0, min(1.0, max(cfg%qmin, q1(i,k,1))/qs)) qstl(i,k) = qs enddo enddo @@ -506,7 +461,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & if (qlx(i,k) > clwt) then onemrh = max(1.e-10, 1.0-rhly(i,k)) tem1 = min(max((onemrh*qstl(i,k))**0.49,0.0001),1.0) - tem1 = cql / tem1 + tem1 = cfg%cql / tem1 value = max(min( tem1*qlx(i,k), 50.0), 0.0) tem2 = sqrt(sqrt(rhly(i,k))) cfly(i,k) = min(max(tem2*(1.0-exp(-value)), 0.0), 1.0) @@ -583,15 +538,15 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & if(pblflg(i)) then ! thermal(i) = thvx(i,1) thermal(i) = thlvx(i,1) - crb(i) = rbcr + crb(i) = cfg%rbcr else - thermal(i) = tsea(i)*(1.+fv*max(q1(i,1,1),qmin)) + thermal(i) = tsea(i)*(1.+fv*max(q1(i,1,1),cfg%qmin)) tem = sqrt(u10m(i)**2+v10m(i)**2) tem = max(tem, 1.) - robn = tem / (f0 * z0(i)) + robn = tem / (cfg%f0 * z0(i)) tem1 = 1.e-7 * robn crb(i) = 0.16 * (tem1 ** (-0.18)) - crb(i) = max(min(crb(i), crbmax), crbmin) + crb(i) = max(min(crb(i), cfg%crbmax), cfg%crbmin) endif enddo !> - Compute \f$\frac{\Delta t}{\Delta z}\f$ , \f$u_*\f$ @@ -612,7 +567,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & ! bf(i,k) = gotvx(i,k)*(thvx(i,k+1)-thvx(i,k))*rdz 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 + shr2(i,k) = max(dw2,cfg%dw2min)*rdz*rdz enddo enddo ! @@ -640,7 +595,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & & (g*zl(i,k)/thlvx(i,1))/spdk2 else if (tc_pbl == 1) then spdk2 = max((u1(i,k)-u1(i,1))**2+(v1(i,k)-v1(i,1))**2,1.) - & + bfac*ustar(i)**2 + & + cfg%bfac*ustar(i)**2 rbup(i) = (thlvx(i,k)-thermal(i))* & (g*(zl(i,k)-zl(i,1))/thlvx(i,1))/spdk2 endif @@ -673,10 +628,10 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & if(kpbl(i) <= 1) pblflg(i)=.false. enddo ! -! update thermal at a level of slfac*hpbl for unstable pbl +! update thermal at a level of cfg%slfac*hpbl for unstable pbl ! do i=1,im - sfcpbl(i) = slfac * hpbl(i) + sfcpbl(i) = cfg%slfac * hpbl(i) kb1(i) = 1 flg(i) = .false. if(pblflg(i)) then @@ -721,7 +676,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & & (g*zl(i,k)/thlvx(i,1))/spdk2 else if (tc_pbl == 1) then spdk2 = max((u1(i,k)-u1(i,1))**2+(v1(i,k)-v1(i,1))**2,1.) - & + bfac*ustar(i)**2 + & + cfg%bfac*ustar(i)**2 rbup(i) = (thlvx(i,k)-thermal(i))* & (g*(zl(i,k)-zl(i,1))/thlvx(i,1))/spdk2 endif @@ -777,7 +732,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & dz = zi(i,k+1) - zi(i,k) tem = (0.5*(u1(i,k-1)-u1(i,k+1))/dz)**2 tem1 = tem+(0.5*(v1(i,k-1)-v1(i,k+1))/dz)**2 - wush(i,k) = csmf * sqrt(tem1) + wush(i,k) = cfg%csmf * sqrt(tem1) enddo enddo ! @@ -791,11 +746,11 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & !! 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) + zol(i) = max(rbsoil(i)*fm(i)*fm(i)/fh(i),cfg%rimin) if(sfcflg(i)) then - zol(i) = min(zol(i),-zfmin) + zol(i) = min(zol(i),-cfg%zfmin) else - zol(i) = max(zol(i),zfmin) + zol(i) = max(zol(i),cfg%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: @@ -808,17 +763,17 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & !! \f[ !! \phi_m=\phi_t=(1+5\frac{0.1h}{L}) !! \f] - zol1 = zol(i)*sfcfrac*hpbl(i)/zl(i,1) + zol1 = zol(i)*cfg%sfcfrac*hpbl(i)/zl(i,1) if(sfcflg(i)) then - tem = 1.0 / (1. - aphi16*zol1) + tem = 1.0 / (1. - cfg%aphi16*zol1) phih(i) = sqrt(tem) phim(i) = sqrt(phih(i)) - tem1 = 1.0 / (1. - aphi16*zol(i)) + tem1 = 1.0 / (1. - cfg%aphi16*zol(i)) phims(i) = sqrt(sqrt(tem1)) else - phim(i) = 1. + aphi5*zol1 + phim(i) = 1. + cfg%aphi5*zol1 phih(i) = phim(i) - phims(i) = 1. + aphi5*zol(i) + phims(i) = 1. + cfg%aphi5*zol(i) endif enddo ! @@ -831,7 +786,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & !! 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), +!! of the surface layer height to the PBL height (specified as cfg%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[ @@ -839,14 +794,15 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & !! \f] do i=1,im if(pblflg(i)) then - if(zol(i) < zolcru) then + if(zol(i) < cfg%zolcru) then pcnvflg(i) = .true. endif wst3(i) = gotvx(i,1)*sflux(i)*hpbl(i) - wstar(i)= wst3(i)**h1 + wstar(i)= wst3(i)**cfg%h1 ust3(i) = ustar(i)**3. - wscale(i)=(ust3(i)+wfac*vk*wst3(i)*sfcfrac)**h1 - ptem = ustar(i)/aphi5 + wscale(i)=(ust3(i)+cfg%wfac*cfg%vk*wst3(i)* + & cfg%sfcfrac)**cfg%h1 + ptem = ustar(i)/cfg%aphi5 wscale(i) = max(wscale(i),ptem) endif enddo @@ -860,7 +816,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & hgamq(i) = evap(i)/wscale(i) vpert(i) = hgamt(i) + hgamq(i)*fv*theta(i,1) vpert(i) = max(vpert(i),0.) - tem = min(cfac*vpert(i),gamcrt) + tem = min(cfg%cfac*vpert(i),cfg%gamcrt) thermal(i)= thermal(i) + tem endif enddo @@ -885,7 +841,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & & (g*zl(i,k)/thlvx(i,1))/spdk2 else if (tc_pbl == 1) then spdk2 = max((u1(i,k)-u1(i,1))**2+(v1(i,k)-v1(i,1))**2,1.) - & + bfac*ustar(i)**2 + & + cfg%bfac*ustar(i)**2 rbup(i) = (thlvx(i,k)-thermal(i))* & (g*(zl(i,k)-zl(i,1))/thlvx(i,1))/spdk2 endif @@ -926,7 +882,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & enddo do k = 1, km1 do i=1,im - if(flg(i).and.zl(i,k) >= zstblmax) then + if(flg(i).and.zl(i,k) >= cfg%zstblmax) then lcld(i)=k flg(i)=.false. endif @@ -938,7 +894,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & do k = kmscu,1,-1 do i = 1, im if(flg(i) .and. k <= lcld(i)) then - if(qlx(i,k) >= qlcr) then + if(qlx(i,k) >= cfg%qlcr) then kcld(i)=k flg(i)=.false. endif @@ -959,7 +915,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & do k = kmscu,1,-1 do i = 1, im if(flg(i) .and. k <= kcld(i)) then - if(qlx(i,k) >= qlcr) then + if(qlx(i,k) >= cfg%qlcr) then if(radx(i,k) < radmin(i)) then radmin(i)=radx(i,k) krad(i)=k @@ -1047,56 +1003,56 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & do i = 1, im if(k < kpbl(i)) then tem = phih(i)/phim(i) - ptem = sfcfrac*hpbl(i) + ptem = cfg%sfcfrac*hpbl(i) tem1 = max(zi(i,k+1)-ptem, 0.) tem2 = tem1 / (hpbl(i) - ptem) if(pcnvflg(i)) then - tem = min(tem, pr0) - prn(i,k) = tem + (pr0 - tem) * tem2 + tem = min(tem, cfg%pr0) + prn(i,k) = tem + (cfg%pr0 - tem) * tem2 else - tem = max(tem, pr0) + tem = max(tem, cfg%pr0) prn(i,k) = tem endif - prn(i,k) = min(prn(i,k),prmax) - prn(i,k) = max(prn(i,k),prmin) + prn(i,k) = min(prn(i,k),cfg%prmax) + prn(i,k) = max(prn(i,k),cfg%prmin) ! - ckz(i,k) = ck0 + (ck1 - ck0) * tem2 - ckz(i,k) = max(min(ckz(i,k), ck0), ck1) - chz(i,k) = ch0 + (ch1 - ch0) * tem2 - chz(i,k) = max(min(chz(i,k), ch0), ch1) + ckz(i,k) = cfg%ck0 + (cfg%ck1 - cfg%ck0) * tem2 + ckz(i,k) = max(min(ckz(i,k), cfg%ck0), cfg%ck1) + chz(i,k) = cfg%ch0 + (cfg%ch1 - cfg%ch0) * tem2 + chz(i,k) = max(min(chz(i,k), cfg%ch0), cfg%ch1) ! endif enddo enddo ! -! Above a threshold height (hcrinv), the background vertical +! Above a threshold height (cfg%hcrinv), the background vertical ! diffusivities & mixing length -! in the inversion layers are set to much smaller values (xkinv1 & -! rlmn1) +! in the inversion layers are set to much smaller values (cfg%xkinv1 & +! cfg%rlmn1) ! -! Below the threshold height (hcrinv), the background vertical +! Below the threshold height (cfg%hcrinv), the background vertical ! diffusivities & mixing length ! in the inversion layers are increased with increasing roughness ! length & vegetation fraction ! do k = 1,km1 do i=1,im - if(zi(i,k+1) > hcrinv) then + if(zi(i,k+1) > cfg%hcrinv) then tem1 = tvx(i,k+1)-tvx(i,k) if(tem1 >= 0.) then - xkzo(i,k) = min(xkzo(i,k), xkinv1) - xkzmo(i,k) = min(xkzmo(i,k), xkinv1) - rlmnz(i,k) = min(rlmnz(i,k), rlmn1) + xkzo(i,k) = min(xkzo(i,k), cfg%xkinv1) + xkzmo(i,k) = min(xkzmo(i,k), cfg%xkinv1) + rlmnz(i,k) = min(rlmnz(i,k), cfg%rlmn1) endif else tem1 = tvx(i,k+1)-tvx(i,k) if(tem1 > 0.) then ptem = xkzo(i,k) * zvfun(i) - xkzo(i,k) = min(max(ptem, xkinv2), xkzo(i,k)) + xkzo(i,k) = min(max(ptem, cfg%xkinv2), xkzo(i,k)) ptem = xkzmo(i,k) * zvfun(i) - xkzmo(i,k) = min(max(ptem, xkinv2), xkzmo(i,k)) + xkzmo(i,k) = min(max(ptem, cfg%xkinv2), xkzmo(i,k)) ptem = rlmnz(i,k) * zvfun(i) - rlmnz(i,k) = min(max(ptem, rlmn2), rlmnz(i,k)) + rlmnz(i,k) = min(max(ptem, cfg%rlmn2), rlmnz(i,k)) endif endif enddo @@ -1119,7 +1075,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & if(mlenflg) then dz = zl(i,n+1) - zl(i,n) tem1 = 2.*gotvx(i,n+1)*(thvx(i,k)-thvx(i,n+1)) - tem2 = cs0*sqrt(e2(i,n))*sqrt(shr2(i,n)) + tem2 = cfg%cs0*sqrt(e2(i,n))*sqrt(shr2(i,n)) e2(i,n+1) = e2(i,n) + (tem1 - tem2) * dz zlup = zlup + dz if(e2(i,n+1) < 0.) then @@ -1136,15 +1092,15 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & if(mlenflg) then if(n == 1) then dz = zl(i,1) - tem = tsea(i)*(1.+fv*max(q1(i,1,1),qmin)) + tem = tsea(i)*(1.+fv*max(q1(i,1,1),cfg%qmin)) tem1 = 2.*gotvx(i,n)*(tem-thvx(i,k)) - tem2 = ustar(i)*phims(i)/(vk*dz) - tem2 = cs0*sqrt(e2(i,n))*tem2 + tem2 = ustar(i)*phims(i)/(cfg%vk*dz) + tem2 = cfg%cs0*sqrt(e2(i,n))*tem2 e2(i,n-1) = e2(i,n) + (tem1 - tem2) * dz else dz = zl(i,n) - zl(i,n-1) tem1 = 2.*gotvx(i,n-1)*(thvx(i,n-1)-thvx(i,k)) - tem2 = cs0*sqrt(e2(i,n))*sqrt(shr2(i,n-1)) + tem2 = cfg%cs0*sqrt(e2(i,n))*sqrt(shr2(i,n-1)) e2(i,n-1) = e2(i,n) + (tem1 - tem2) * dz endif zldn = zldn + dz @@ -1176,12 +1132,12 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & ! mixing length was included. ! ptem2 = min(zlup,zldn) - rlam(i,k) = elmfac * ptem2 + rlam(i,k) = cfg%elmfac * ptem2 rlam(i,k) = max(rlam(i,k), tem1) rlam(i,k) = min(rlam(i,k), rlmx) ! ptem2 = sqrt(zlup*zldn) - ele(i,k) = elefac * ptem2 + ele(i,k) = cfg%elefac * ptem2 ele(i,k) = max(ele(i,k), tem1) ele(i,k) = min(ele(i,k), elmx) ! @@ -1191,7 +1147,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & !! 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) + tem = cfg%vk * zl(i,k) if (zol(i) < 0.) then ptem = 1. - 100. * zol(i) ptem1 = ptem**0.2 @@ -1243,7 +1199,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & do i = 1, im tem = 0.5 * (elm(i,k) + elm(i,k+1)) tem = tem * sqrt(tkeh(i,k)) - ri = max(bf(i,k)/shr2(i,k),rimin) + ri = max(bf(i,k)/shr2(i,k),cfg%rimin) if(k < kpbl(i)) then if(pcnvflg(i)) then dku(i,k) = ckz(i,k) * tem @@ -1259,12 +1215,12 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & endif else if(ri < 0.) then ! unstable regime - dku(i,k) = ck1 * tem - dkt(i,k) = rchck * dku(i,k) + dku(i,k) = cfg%ck1 * tem + dkt(i,k) = cfg%rchck * dku(i,k) else ! stable regime - dkt(i,k) = ch1 * tem + dkt(i,k) = cfg%ch1 * tem prnum = 1.0 + 2.1 * ri - prnum = min(prnum,prmax) + prnum = min(prnum,cfg%prmax) dku(i,k) = dkt(i,k) * prnum endif endif @@ -1272,19 +1228,19 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & if(scuflg(i)) then if(k >= mrad(i) .and. k < krad(i)) then tem1 = ckz(i,k) * tem - ptem1 = tem1 / prscu + ptem1 = tem1 / cfg%prscu dku(i,k) = max(dku(i,k), tem1) dkt(i,k) = max(dkt(i,k), ptem1) endif endif ! - dkq(i,k) = prtke * dkt(i,k) + dkq(i,k) = cfg%prtke * dkt(i,k) ! - dkt(i,k) = min(dkt(i,k),dkmax) + dkt(i,k) = min(dkt(i,k),cfg%dkmax) dkt(i,k) = max(dkt(i,k),xkzo(i,k)) - dkq(i,k) = min(dkq(i,k),dkmax) + dkq(i,k) = min(dkq(i,k),cfg%dkmax) dkq(i,k) = max(dkq(i,k),xkzo(i,k)) - dku(i,k) = min(dku(i,k),dkmax) + dku(i,k) = min(dku(i,k),cfg%dkmax) dku(i,k) = max(dku(i,k),xkzmo(i,k)) ! enddo @@ -1303,8 +1259,8 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & endif ptem = tem1 / (tem * elm(i,k)) tkmnz(i,k) = ptem * ptem - tkmnz(i,k) = min(tkmnz(i,k), tkbmx) - tkmnz(i,k) = max(tkmnz(i,k), tkmin) + tkmnz(i,k) = min(tkmnz(i,k), cfg%tkbmx) + tkmnz(i,k) = max(tkmnz(i,k), cfg%tkmin) enddo enddo ! @@ -1361,7 +1317,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & endif ptem2 = ptem2 + ptem ! - tem2 = stress(i)*ustar(i)*phims(i)/(vk*zl(i,1)) + tem2 = stress(i)*ustar(i)*phims(i)/(cfg%vk*zl(i,1)) shrp = 0.5 * (tem1 + ptem1 + ptem2 + tem2) else tem1 = -dkt(i,k-1) * bf(i,k-1) @@ -1434,17 +1390,17 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & !---------------------------------------------------------------------- !> - First predict tke due to tke production & dissipation(diss) ! - dtn = dt2 / float(ndt) - do n = 1, ndt + dtn = dt2 / cfg%ndt + do n = 1, cfg%ndt do k = 1,km1 do i=1,im tem = sqrt(tke(i,k)) - ptem = ce0 / ele(i,k) + ptem = cfg%ce0 / ele(i,k) diss(i,k) = ptem * tke(i,k) * tem tem1 = prod(i,k) + tke(i,k) / dtn diss(i,k)=max(min(diss(i,k), tem1), 0.) tke(i,k) = tke(i,k) + dtn * (prod(i,k)-diss(i,k)) -! tke(i,k) = max(tke(i,k), tkmin) +! tke(i,k) = max(tke(i,k), cfg%tkmin) tke(i,k) = max(tke(i,k), tkmnz(i,k)) enddo enddo @@ -1776,7 +1732,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & c do k = 1,km do i = 1,im -! f1(i,k) = max(f1(i,k), tkmin) +! f1(i,k) = max(f1(i,k), cfg%tkmin) qtend = (f1(i,k)-q1(i,k,ntke))*rdt rtg(i,k,ntke) = rtg(i,k,ntke)+qtend enddo @@ -2282,7 +2238,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, & if(dspheat) then do k = 1,km1 do i = 1,im -! tem = min(diss(i,k), dspmax) +! tem = min(diss(i,k), cfg%dspmax) ! ttend = tem / cp ttend = diss(i,k) / cp tdt(i,k) = tdt(i,k) + dspfac * ttend diff --git a/physics/satmedmfvdifq.meta b/physics/satmedmfvdifq.meta index b6680dccb..aef1d5ab7 100644 --- a/physics/satmedmfvdifq.meta +++ b/physics/satmedmfvdifq.meta @@ -2,6 +2,7 @@ name = satmedmfvdifq type = scheme dependencies = funcphys.f90,machine.F,mfpbltq.f,mfscuq.f,tridi.f + dependencies = module_satmedmf.F90 ######################################################################## [ccpp-arg-table] @@ -48,6 +49,13 @@ [ccpp-arg-table] name = satmedmfvdifq_run type = scheme +[cfg] + standard_name = parameters_for_scale_aware_TKE_moist_EDMF_PBL_2019 + long_name = parameters for scale-aware TKE moist EDMF PBL scheme + units = mixed + dimensions = () + type = ty_cfg_satmedmfq + intent = in [im] standard_name = horizontal_loop_extent long_name = horizontal loop extent