From 805c62c1a89b867b676a50555dc43f323fe1bf56 Mon Sep 17 00:00:00 2001 From: Grant Firl Date: Tue, 7 Dec 2021 10:39:05 -0700 Subject: [PATCH 01/20] add single precision code changes from michalakes fork, jm-nrl-32bitfp-24cc09e branch --- physics/calpreciptype.f90 | 77 +++++++++-------- physics/funcphys.f90 | 138 ++++++++++++++++++++++++++----- physics/machine.F | 14 ++-- physics/module_bl_mynn.F90 | 14 ++-- physics/radlw_main.F90 | 10 ++- physics/radsw_main.F90 | 9 +- physics/sfc_diag_post.F90 | 9 +- physics/surface_perturbation.F90 | 2 +- 8 files changed, 197 insertions(+), 76 deletions(-) diff --git a/physics/calpreciptype.f90 b/physics/calpreciptype.f90 index dcc8ed49b..d3fbb253b 100644 --- a/physics/calpreciptype.f90 +++ b/physics/calpreciptype.f90 @@ -26,17 +26,18 @@ subroutine calpreciptype(kdt,nrcm,im,ix,lm,lp1,randomno, & ! -------------------------------------------------------------------- use funcphys, only : fpvs,ftdp,fpkap,ftlcl,stma,fthe use physcons + use machine , only : kind_phys !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none ! - real, parameter :: pthresh = 0.0, oneog = 1.0/con_g + real(kind=kind_phys), parameter :: pthresh = 0.0, oneog = 1.0/con_g integer,parameter :: nalg = 5 ! ! declare variables. ! integer,intent(in) :: kdt,nrcm,im,ix,lm,lp1 - real,intent(in) :: xlat(im),xlon(im) - real,intent(in) :: randomno(ix,nrcm) + real(kind=kind_phys),intent(in) :: xlat(im),xlon(im) + real(kind=kind_phys),intent(in) :: randomno(ix,nrcm) real(kind=kind_phys),dimension(im), intent(in) :: prec,tskin real(kind=kind_phys),dimension(ix,lm), intent(in) :: gt0,gq0,prsl real(kind=kind_phys),dimension(ix,lp1),intent(in) :: prsi,phii @@ -220,8 +221,9 @@ subroutine calpreciptype(kdt,nrcm,im,ix,lm,lp1,randomno, & !! This subroutine computes precipitation type using a decision tree approach that uses !! variables such as integrated wet bulb temperatue below freezing and lowest layer !! temperature (Baldwin et al. 1994 \cite baldwin_et_al_1994) - subroutine calwxt(lm,lp1,t,q,pmid,pint, & - d608,rog,epsq,zint,iwx,twet) + subroutine calwxt(lm,lp1,t,q,pmid,pint, & + d608,rog,epsq,zint,iwx,twet) + use machine , only : kind_phys ! ! file: calwxt.f ! written: 11 november 1993, michael baldwin @@ -247,10 +249,10 @@ subroutine calwxt(lm,lp1,t,q,pmid,pint, & ! t,q,pmid,htm,lmh,zint ! integer,intent(in) :: lm,lp1 - real,dimension(lm),intent(in) :: t,q,pmid,twet - real,dimension(lp1),intent(in) :: zint,pint + real(kind=kind_phys),dimension(lm),intent(in) :: t,q,pmid,twet + real(kind=kind_phys),dimension(lp1),intent(in) :: zint,pint integer,intent(out) :: iwx - real,intent(in) :: d608,rog,epsq + real(kind=kind_phys),intent(in) :: d608,rog,epsq ! output: @@ -264,10 +266,10 @@ subroutine calwxt(lm,lp1,t,q,pmid,pint, & ! ! internal: ! -! real, allocatable :: twet(:) - real, parameter :: d00=0.0 +! real(kind=kind_phys), allocatable :: twet(:) + real(kind=kind_phys), parameter :: d00=0.0 integer karr,licee - real tcold,twarm + real(kind=kind_phys) tcold,twarm ! subroutines called: ! wetbulb @@ -282,7 +284,7 @@ subroutine calwxt(lm,lp1,t,q,pmid,pint, & ! integer l,lice,iwrml,ifrzl - real psfck,tdchk,a,tdkl,tdpre,tlmhk,twrmk,areas8,areap4, & + real(kind=kind_phys) psfck,tdchk,a,tdkl,tdpre,tlmhk,twrmk,areas8,areap4, & surfw,surfc,dzkl,area1,pintk1,pintk2,pm150,pkl,tkl,qkl ! allocate ( twet(lm) ) @@ -486,27 +488,28 @@ subroutine calwxt_ramer(lm,lp1,t,q,pmid,rh,td,pint,ptyp) ! use params_mod ! use ctlblk_mod !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + use machine , only : kind_phys implicit none ! - real,parameter :: twice=266.55,rhprcp=0.80,deltag=1.02, & + real(kind=kind_phys),parameter :: twice=266.55,rhprcp=0.80,deltag=1.02, & & emelt=0.045,rlim=0.04,slim=0.85 - real,parameter :: twmelt=273.15,tz=273.15,efac=1.0 ! specify in params now + real(kind=kind_phys),parameter :: twmelt=273.15,tz=273.15,efac=1.0 ! specify in params now ! integer*4 i, k1, lll, k2, toodry ! - real xxx ,mye, icefrac + real(kind=kind_phys) xxx ,mye, icefrac integer, intent(in) :: lm,lp1 - real,dimension(lm), intent(in) :: t,q,pmid,rh,td - real,dimension(lp1),intent(in) :: pint + real(kind=kind_phys),dimension(lm), intent(in) :: t,q,pmid,rh,td + real(kind=kind_phys),dimension(lp1),intent(in) :: pint integer, intent(out) :: ptyp ! - real,dimension(lm) :: tq,pq,rhq,twq + real(kind=kind_phys),dimension(lm) :: tq,pq,rhq,twq ! integer j,l,lev,ii - real rhmax,twmax,ptop,dpdrh,twtop,rhtop,wgt1,wgt2, & + real(kind=kind_phys) rhmax,twmax,ptop,dpdrh,twtop,rhtop,wgt1,wgt2, & rhavg,dtavg,dpk,ptw,pbot -! real b,qtmp,rate,qc - real,external :: xmytw +! real(kind=kind_phys) b,qtmp,rate,qc +! real(kind=kind_phys),external :: xmytw (now inside the module) ! ! initialize. icefrac = -9999. @@ -521,7 +524,7 @@ subroutine calwxt_ramer(lm,lp1,t,q,pmid,rh,td,pint,ptyp) ! causing problems later in this subroutine ! qtmp=max(h1m12,q(l)) ! rhqtmp(lev)=qtmp/qc - rhq(lev) = rh(l) + rhq(lev) = rh(l) pq(lev) = pmid(l) * 0.01 tq(lev) = t(l) enddo @@ -753,10 +756,11 @@ subroutine calwxt_ramer(lm,lp1,t,q,pmid,rh,td,pint,ptyp) !-------------------------------------------------------------------------- function xmytw(t,td,p) ! + use machine , only : kind_phys implicit none ! integer*4 cflag, l - real f, c0, c1, c2, k, kd, kw, ew, t, td, p, ed, fp, s, & + real(kind=kind_phys) f, c0, c1, c2, k, kd, kw, ew, t, td, p, ed, fp, s, & & de, xmytw data f, c0, c1, c2 /0.0006355, 26.66082, 0.0091379024, 6106.3960/ ! @@ -877,19 +881,20 @@ function xmytw(t,td,p) !! \cite bourgouin_2000. !of aes (canada) 1992 subroutine calwxt_bourg(lm,lp1,rn,g,t,q,pmid,pint,zint,ptype) + use machine , only : kind_phys implicit none ! ! input: integer,intent(in) :: lm,lp1 - real,intent(in) :: g,rn(2) - real,intent(in), dimension(lm) :: t, q, pmid - real,intent(in), dimension(lp1) :: pint, zint + real(kind=kind_phys),intent(in) :: g,rn(2) + real(kind=kind_phys),intent(in), dimension(lm) :: t, q, pmid + real(kind=kind_phys),intent(in), dimension(lp1) :: pint, zint ! ! output: integer, intent(out) :: ptype ! integer ifrzl,iwrml,l,lhiwrm - real pintk1,areane,tlmhk,areape,pintk2,surfw,area1,dzkl,psfck,r1,r2 + real(kind=kind_phys) pintk1,areane,tlmhk,areape,pintk2,surfw,area1,dzkl,psfck,r1,r2 ! ! initialize weather type array to zero (ie, off). ! we do this since we want ptype to represent the @@ -1076,6 +1081,7 @@ subroutine calwxt_revised(lm,lp1,t,q,pmid,pint, & ! use params_mod ! use ctlblk_mod !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + use machine , only : kind_phys implicit none ! ! list of variables needed @@ -1087,9 +1093,9 @@ subroutine calwxt_revised(lm,lp1,t,q,pmid,pint, & ! t,q,pmid,htm,lmh,zint integer,intent(in) :: lm,lp1 - real,dimension(lm),intent(in) :: t,q,pmid,twet - real,dimension(lp1),intent(in) :: pint,zint - real,intent(in) :: d608,rog,epsq + real(kind=kind_phys),dimension(lm),intent(in) :: t,q,pmid,twet + real(kind=kind_phys),dimension(lp1),intent(in) :: pint,zint + real(kind=kind_phys),intent(in) :: d608,rog,epsq ! output: ! iwx - instantaneous weather type. ! acts like a 4 bit binary @@ -1101,12 +1107,12 @@ subroutine calwxt_revised(lm,lp1,t,q,pmid,pint, & integer, intent(out) :: iwx ! internal: ! - real, parameter :: d00=0.0 + real(kind=kind_phys), parameter :: d00=0.0 integer karr,licee - real tcold,twarm + real(kind=kind_phys) tcold,twarm ! integer l,lmhk,lice,iwrml,ifrzl - real psfck,tdchk,a,tdkl,tdpre,tlmhk,twrmk,areas8,areap4,area1, & + real(kind=kind_phys) psfck,tdchk,a,tdkl,tdpre,tlmhk,twrmk,areas8,areap4,area1, & surfw,surfc,dzkl,pintk1,pintk2,pm150,qkl,tkl,pkl,area0,areap0 ! subroutines called: @@ -1316,14 +1322,15 @@ subroutine calwxt_dominant(nalg,rain,freezr,sleet,snow, & ! algorithms and sums them up to give a dominant type ! !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + use machine , only : kind_phys implicit none ! ! input: integer,intent(in) :: nalg - real,intent(out) :: doms,domr,domzr,domip + real(kind=kind_phys),intent(out) :: doms,domr,domzr,domip integer,dimension(nalg),intent(in) :: rain,snow,sleet,freezr integer l - real totsn,totip,totr,totzr + real(kind=kind_phys) totsn,totip,totr,totzr !-------------------------------------------------------------------------- ! print* , 'into dominant' domr = 0. diff --git a/physics/funcphys.f90 b/physics/funcphys.f90 index 8cb4b1b15..3e81a0d5a 100644 --- a/physics/funcphys.f90 +++ b/physics/funcphys.f90 @@ -260,7 +260,7 @@ module funcphys ! Language: Fortran 90 ! !$$$ - use machine,only:kind_phys + use machine,only:kind_phys,r8=>kind_dbl_prec,r4=>kind_sngl_prec use physcons implicit none private @@ -308,6 +308,13 @@ module funcphys public grkap,frkap,frkapq,frkapx public gtlcl,ftlcl,ftlclq,ftlclo,ftlclx public gfuncphys + + interface fpvsl + module procedure fpvsl_r4, fpvsl_r8 + end interface fpvsl + interface fpvsi + module procedure fpvsi_r4, fpvsi_r8 + end interface fpvsi contains !------------------------------------------------------------------------------- !> This subroutine computes saturation vapor pressure table as a function of @@ -364,7 +371,8 @@ subroutine gpvsl !! in gpvsl(). See documentation for fpvslx() for details. Input values !! outside table range are reset to table extrema. !>\author N phillips - elemental function fpvsl(t) + + elemental function fpvsl_r4(t) !$$$ Subprogram Documentation Block ! ! Subprogram: fpvsl Compute saturation vapor pressure over liquid @@ -396,16 +404,62 @@ elemental function fpvsl(t) ! !$$$ implicit none - real(krealfp) fpvsl - real(krealfp),intent(in):: t + real(r4) fpvsl_r4 + real(r4),intent(in):: t integer jx - real(krealfp) xj + real(r4) xj ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - xj=min(max(c1xpvsl+c2xpvsl*t,1._krealfp),real(nxpvsl,krealfp)) - jx=min(xj,nxpvsl-1._krealfp) - fpvsl=tbpvsl(jx)+(xj-jx)*(tbpvsl(jx+1)-tbpvsl(jx)) + xj=min(max(c1xpvsl+c2xpvsl*t,1._r4),real(nxpvsl,r4)) + jx=min(xj,nxpvsl-1._r4) + fpvsl_r4=tbpvsl(jx)+(xj-jx)*(tbpvsl(jx+1)-tbpvsl(jx)) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function + end function fpvsl_r4 + + elemental function fpvsl_r8(t) +!$$$ Subprogram Documentation Block +! +! Subprogram: fpvsl Compute saturation vapor pressure over liquid +! Author: N Phillips w/NMC2X2 Date: 30 dec 82 +! +! Abstract: Compute saturation vapor pressure from the temperature. +! A linear interpolation is done between values in a lookup table +! computed in gpvsl. See documentation for fpvslx for details. +! Input values outside table range are reset to table extrema. +! The interpolation accuracy is almost 6 decimal places. +! On the Cray, fpvsl is about 4 times faster than exact calculation. +! This function should be expanded inline in the calling routine. +! +! Program History Log: +! 91-05-07 Iredell made into inlinable function +! 94-12-30 Iredell expand table +! 1999-03-01 Iredell f90 module +! +! Usage: pvsl=fpvsl(t) +! +! Input argument list: +! t Real(krealfp) temperature in Kelvin +! +! Output argument list: +! fpvsl Real(krealfp) saturation vapor pressure in Pascals +! +! Attributes: +! Language: Fortran 90. +! +!$$$ + implicit none + real(r8) fpvsl_r8 + real(r8),intent(in):: t + integer jx + real(r8) xj +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + xj=min(max(c1xpvsl+c2xpvsl*t,1._r8),real(nxpvsl,r8)) + jx=min(xj,nxpvsl-1._r8) + fpvsl_r8=tbpvsl(jx)+(xj-jx)*(tbpvsl(jx+1)-tbpvsl(jx)) +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + end function fpvsl_r8 + + + !------------------------------------------------------------------------------- !> This function computes saturation vapor pressure from the temperature. !! A quadratic interpolation is done between values in a lookup table @@ -576,7 +630,8 @@ subroutine gpvsi !! computed in gpvsi(). See documentation for fpvsix() for details. !! Input values outside table range are reset to table extrema. !>\author N Phillips - elemental function fpvsi(t) + + elemental function fpvsi_r4(t) !$$$ Subprogram Documentation Block ! ! Subprogram: fpvsi Compute saturation vapor pressure over ice @@ -609,16 +664,61 @@ elemental function fpvsi(t) ! !$$$ implicit none - real(krealfp) fpvsi - real(krealfp),intent(in):: t + real(r4) fpvsi_r4 + real(r4),intent(in):: t integer jx - real(krealfp) xj + real(r4) xj ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - xj=min(max(c1xpvsi+c2xpvsi*t,1._krealfp),real(nxpvsi,krealfp)) - jx=min(xj,nxpvsi-1._krealfp) - fpvsi=tbpvsi(jx)+(xj-jx)*(tbpvsi(jx+1)-tbpvsi(jx)) + xj=min(max(c1xpvsi+c2xpvsi*t,1._r4),real(nxpvsi,r4)) + jx=min(xj,nxpvsi-1._r4) + fpvsi_r4=tbpvsi(jx)+(xj-jx)*(tbpvsi(jx+1)-tbpvsi(jx)) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function + end function fpvsi_r4 + + elemental function fpvsi_r8(t) +!$$$ Subprogram Documentation Block +! +! Subprogram: fpvsi Compute saturation vapor pressure over ice +! Author: N Phillips w/NMC2X2 Date: 30 dec 82 +! +! Abstract: Compute saturation vapor pressure from the temperature. +! A linear interpolation is done between values in a lookup table +! computed in gpvsi. See documentation for fpvsix for details. +! Input values outside table range are reset to table extrema. +! The interpolation accuracy is almost 6 decimal places. +! On the Cray, fpvsi is about 4 times faster than exact calculation. +! This function should be expanded inline in the calling routine. +! +! Program History Log: +! 91-05-07 Iredell made into inlinable function +! 94-12-30 Iredell expand table +! 1999-03-01 Iredell f90 module +! 2001-02-26 Iredell ice phase +! +! Usage: pvsi=fpvsi(t) +! +! Input argument list: +! t Real(krealfp) temperature in Kelvin +! +! Output argument list: +! fpvsi Real(krealfp) saturation vapor pressure in Pascals +! +! Attributes: +! Language: Fortran 90. +! +!$$$ + implicit none + real(r8) fpvsi_r8 + real(r8),intent(in):: t + integer jx + real(r8) xj +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + xj=min(max(c1xpvsi+c2xpvsi*t,1._r8),real(nxpvsi,r8)) + jx=min(xj,nxpvsi-1._r8) + fpvsi_r8=tbpvsi(jx)+(xj-jx)*(tbpvsi(jx+1)-tbpvsi(jx)) +! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + end function fpvsi_r8 + !------------------------------------------------------------------------------- !> This function computes saturation vapor pressure from the temperature. !! A quadratic interpolation is done between values in a lookup table @@ -2375,7 +2475,7 @@ elemental subroutine stmaq(the,pk,tma,qma) !>\param[in] pk real, pressure over 1e5 Pa to the kappa power !>\param[out] tma real, parcel temperature in Kelvin !>\param[out] qma real, parcel specific humidity in kg/kg - elemental subroutine stmax(the,pk,tma,qma) + subroutine stmax(the,pk,tma,qma) !$$$ Subprogram Documentation Block ! ! Subprogram: stmax Compute moist adiabat temperature @@ -2443,7 +2543,7 @@ elemental subroutine stmax(the,pk,tma,qma) !>\param[in] pk real, pressure over 1e5 Pa to the kappa power !>\param[out] tma real, parcel temperature in Kelvin !>\param[out] qma real, parcel specific humidity in kg/kg - elemental subroutine stmaxg(tg,the,pk,tma,qma) + subroutine stmaxg(tg,the,pk,tma,qma) !$$$ Subprogram Documentation Block ! ! Subprogram: stmaxg Compute moist adiabat temperature diff --git a/physics/machine.F b/physics/machine.F index 896b665da..2ee7fb865 100644 --- a/physics/machine.F +++ b/physics/machine.F @@ -9,11 +9,12 @@ module machine #ifndef SINGLE_PREC integer, parameter :: kind_io4 = 4, kind_io8 = 8 , kind_ior = 8 & &, kind_evod = 8, kind_dbl_prec = 8 & -#ifdef __PGI + &, kind_sngl_prec = 4 +# ifdef __PGI &, kind_qdt_prec = 8 & -#else +# else &, kind_qdt_prec = 16 & -#endif +# endif &, kind_rad = 8 & &, kind_phys = 8 ,kind_taum=8 & &, kind_grid = 8 & @@ -24,11 +25,12 @@ module machine #else integer, parameter :: kind_io4 = 4, kind_io8 = 8 , kind_ior = 8 & &, kind_evod = 4, kind_dbl_prec = 8 & -#ifdef __PGI + &, kind_sngl_prec = 4 +# ifdef __PGI &, kind_qdt_prec = 8 & -#else +# else &, kind_qdt_prec = 16 & -#endif +# endif &, kind_rad = 4 & &, kind_phys = 4 ,kind_taum=4 & &, kind_grid = 4 & diff --git a/physics/module_bl_mynn.F90 b/physics/module_bl_mynn.F90 index d691de909..ff9574a27 100644 --- a/physics/module_bl_mynn.F90 +++ b/physics/module_bl_mynn.F90 @@ -142,6 +142,7 @@ MODULE module_bl_mynn & XLF => con_hfus, & & EP_1 => con_fvirt, & & EP_2 => con_eps + use machine, only : kind_phys IMPLICIT NONE @@ -1470,8 +1471,11 @@ SUBROUTINE boulac_length(kts,kte,zw,dz,qtke,theta,lb1,lb2) dld(iz) = min(dld(iz),zw(iz+1))!not used in PBL anyway, only free atmos lb1(iz) = min(dlu(iz),dld(iz)) !minimum !JOE-fight floating point errors +#ifdef SINGLE_PREC + !JM: keep up the fight, JOE dlu(iz)=MAX(0.1,MIN(dlu(iz),1000.)) dld(iz)=MAX(0.1,MIN(dld(iz),1000.)) +#endif lb2(iz) = sqrt(dlu(iz)*dld(iz)) !average - biased towards smallest !lb2(iz) = 0.5*(dlu(iz)+dld(iz)) !average @@ -2692,7 +2696,7 @@ SUBROUTINE mym_condensation (kts,kte, & !CLOUD WATER AND ICE IF (q1k < 0.) THEN !unstaurated - ql_water = sgm(k)*EXP(1.2*q1k-1) + ql_water = sgm(k)*EXP(1.2*q1k-1.) ql_ice = sgm(k)*EXP(1.2*q1k-1.) !Reduce ice mixing ratios in the upper troposphere ! low_weight = MIN(MAX(p(k)-40000.0, 0.0),40000.0)/40000.0 @@ -6723,15 +6727,15 @@ FUNCTION qsat_blend(t, P, waterice) IF ((t .GE. 273.16) .OR. (wrt .EQ. 'w')) THEN ESL = J0+XC*(J1+XC*(J2+XC*(J3+XC*(J4+XC*(J5+XC*(J6+XC*(J7+XC*J8))))))) - qsat_blend = 0.622*ESL/(P-ESL) + qsat_blend = 0.622*ESL/max((P-ESL),1.0E-7_kind_phys) ELSE IF (t .LE. 253.) THEN ESI = K0+XC*(K1+XC*(K2+XC*(K3+XC*(K4+XC*(K5+XC*(K6+XC*(K7+XC*K8))))))) - qsat_blend = 0.622*ESI/(P-ESI) + qsat_blend = 0.622*ESI/max((P-ESI),1.0E-7_kind_phys) ELSE ESL = J0+XC*(J1+XC*(J2+XC*(J3+XC*(J4+XC*(J5+XC*(J6+XC*(J7+XC*J8))))))) ESI = K0+XC*(K1+XC*(K2+XC*(K3+XC*(K4+XC*(K5+XC*(K6+XC*(K7+XC*K8))))))) - RSLF = 0.622*ESL/(P-ESL) - RSIF = 0.622*ESI/(P-ESI) + RSLF = 0.622*ESL/max((P-ESL),1.0E-7_kind_phys) + RSIF = 0.622*ESI/max((P-ESI),1.0E-7_kind_phys) chi = (273.16-t)/20.16 qsat_blend = (1.-chi)*RSLF + chi*RSIF END IF diff --git a/physics/radlw_main.F90 b/physics/radlw_main.F90 index 89609c283..b6e41b094 100644 --- a/physics/radlw_main.F90 +++ b/physics/radlw_main.F90 @@ -286,7 +286,8 @@ module rrtmg_lw & random_stat !mz use machine, only : kind_phys, & - & im => kind_io4, rb => kind_phys + & im => kind_io4, rb => kind_phys, & + & kind_dbl_prec use module_radlw_parameters ! @@ -2071,9 +2072,10 @@ subroutine mcica_subcol & logical, dimension(ngptlw,nlay), intent(out) :: lcloudy ! --- locals: - real (kind=kind_phys) :: cdfunc(ngptlw,nlay), rand1d(ngptlw), & - & rand2d(nlay*ngptlw), tem1, fac_lcf(nlay), & + real (kind=kind_phys) :: cdfunc(ngptlw,nlay), & + & tem1, fac_lcf(nlay), & & cdfun2(ngptlw,nlay) + real (kind=kind_dbl_prec) rand2d(nlay*ngptlw), rand1d(ngptlw) type (random_stat) :: stat ! for thread safe random generator @@ -8968,4 +8970,4 @@ end subroutine cldprmc !........................................!$ end module rrtmg_lw !$ -!========================================!$ \ No newline at end of file +!========================================!$ diff --git a/physics/radsw_main.F90 b/physics/radsw_main.F90 index 0f5a8b110..32097d868 100644 --- a/physics/radsw_main.F90 +++ b/physics/radsw_main.F90 @@ -310,7 +310,7 @@ module rrtmg_sw use physcons, only : con_g, con_cp, con_avgd, con_amd, & & con_amw, con_amo3 use machine, only : rb => kind_phys, im => kind_io4, & - & kind_phys + & kind_phys, kind_dbl_prec use module_radsw_parameters use mersenne_twister, only : random_setseed, random_number, & @@ -1733,6 +1733,10 @@ subroutine rswinit & tfn = float(i) / float(NTBMX-i) tau = bpade * tfn exp_tbl(i) = exp( -tau ) +#ifdef SINGLE_PREC + ! from WRF version, prevents zero at single prec + if (exp_tbl(i) .le. expeps) exp_tbl(i) = expeps +#endif enddo return @@ -2213,8 +2217,9 @@ subroutine mcica_subcol & ! --- locals: real (kind=kind_phys) :: cdfunc(nlay,ngptsw), tem1, & - & rand2d(nlay*ngptsw), rand1d(ngptsw), fac_lcf(nlay), & + & fac_lcf(nlay), & & cdfun2(nlay,ngptsw) + real (kind=kind_dbl_prec) :: rand2d(nlay*ngptsw), rand1d(ngptsw) type (random_stat) :: stat ! for thread safe random generator diff --git a/physics/sfc_diag_post.F90 b/physics/sfc_diag_post.F90 index 6f14fe93d..26f4f1ba8 100644 --- a/physics/sfc_diag_post.F90 +++ b/physics/sfc_diag_post.F90 @@ -19,7 +19,7 @@ subroutine sfc_diag_post_run (im, lsm, lsm_noahmp, dry, lssav, dtf, con_eps, con t2m, q2m, u10m, v10m, tmpmin, tmpmax, spfhmin, spfhmax, & wind10mmax, u10mmax, v10mmax, dpt2m, errmsg, errflg) - use machine, only: kind_phys + use machine, only: kind_phys, kind_dbl_prec implicit none @@ -35,7 +35,7 @@ subroutine sfc_diag_post_run (im, lsm, lsm_noahmp, dry, lssav, dtf, con_eps, con integer, intent(out) :: errflg integer :: i - real(kind=kind_phys) :: tem + real(kind=kind_dbl_prec) :: tem ! made dbl prec always, JM 20211104 ! Initialize CCPP error handling variables errmsg = '' @@ -57,8 +57,9 @@ subroutine sfc_diag_post_run (im, lsm, lsm_noahmp, dry, lssav, dtf, con_eps, con v10mmax(i) = v10m(i) endif ! Compute dew point, first using vapor pressure - tem = max(pgr(i) * q2m(i) / ( con_eps - con_epsm1 *q2m(i)), 1.e-8) - dpt2m(i) = 243.5 / ( ( 17.67 / log(tem/611.2) ) - 1.) + 273.14 + tem = max(pgr(i) * q2m(i) / ( con_eps - con_epsm1 *q2m(i)), 1.d-8) + dpt2m(i) = 243.5_kind_dbl_prec / & + ( ( 17.67_kind_dbl_prec / log(tem/611.2_kind_dbl_prec) ) - 1.) + 273.14 enddo endif diff --git a/physics/surface_perturbation.F90 b/physics/surface_perturbation.F90 index e0429a5fc..7ddbe5279 100644 --- a/physics/surface_perturbation.F90 +++ b/physics/surface_perturbation.F90 @@ -48,7 +48,7 @@ subroutine cdfnor(z,cdfz) cdfz = 0.5 else x = 0.5*z*z - call cdfgam(x,0.5,del,iflag, cdfx) + call cdfgam(x,0.5_kind_phys,del,iflag, cdfx) if (iflag.ne.0) return if (z.gt.0.0) then cdfz = 0.5+0.5*cdfx From 235ec3825715dfb646afe00dc71b0c010ff9b661 Mon Sep 17 00:00:00 2001 From: Lisa Bengtsson Date: Wed, 13 Apr 2022 17:44:27 +0000 Subject: [PATCH 02/20] add progsigma_calc --- physics/progsigma_calc.f90 | 260 +++++++++++++++++++++++++++++++++++++ 1 file changed, 260 insertions(+) create mode 100644 physics/progsigma_calc.f90 diff --git a/physics/progsigma_calc.f90 b/physics/progsigma_calc.f90 new file mode 100644 index 000000000..378d43ef4 --- /dev/null +++ b/physics/progsigma_calc.f90 @@ -0,0 +1,260 @@ +!>\file progsigma +!! This file contains the subroutine that calculates the prognostic +!! updraft area fraction that is used for closure computations in +!! saSAS deep and shallow convection. + +!>\ingroup samfdeepcnv +!! This subroutine computes a prognostic updraft area fraction +!! used in the closure computations in the samfdeepcnv.f scheme +!>\ingroup samfshalcnv +!! This subroutine computes a prognostic updraft area fracftion +!! used in the closure computations in the samfshalcnv. scheme +!!\section progsigma General Algorithm +!> @{ + + subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_deep, & + del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt, & + qgrs_dsave,q,kbcon1,ktcon,cnvflg,gdx, & + do_ca, ca_closure, ca_entr, ca_trigger, nthresh, ca_deep, & + ca_turb,ca_micro,ca_shal,ca_rad,convcount,ca1,ca2,ca3,ca4, & + sigmain,sigmaout,sigmab,errmsg,errflg) +! +! + use machine, only : kind_phys + use funcphys, only : fpvs + + implicit none + +! intent in + integer, intent(in) :: im,km,kbcon1(im),ktcon(im) + real, intent(in) :: hvap,delt + real, intent(in) :: qgrs_dsave(im,km), q(im,km),del(im,km), & + qmicro(im,km),tmf(im,km),dbyo1(im,km),zdqca(im,km), & + omega_u(im,km),zeta(im,km),gdx(im) + logical, intent(in) :: flag_init,flag_restart,flag_deep,cnvflg(im) + real(kind=kind_phys), intent(in) :: nthresh + real(kind=kind_phys), intent(in) :: ca_deep(im) + real(kind=kind_phys), intent(out):: ca_turb(im), & + ca_micro(im),ca_rad(im),ca_shal(im),convcount(im),ca1(im), & + ca2(im),ca3(im),ca4(im) + logical, intent(in) :: do_ca,ca_closure,ca_entr,ca_trigger + + real(kind=kind_phys), intent(in) :: sigmain(im,km) + +! intent out + real(kind=kind_phys), intent(out) :: sigmaout(im,km) + real(kind=kind_phys), intent(out) :: sigmab(im) + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + +! Local variables + integer :: i,k,km1 + real(kind=kind_phys) :: termA(im),termB(im),termC(im),termD(im), & + mcons(im),zfdqa(im),zform(im,km), & + qadv(im,km),sigmamax(im) + + + real(kind=kind_phys) :: gcvalmx,ZEPS7,ZZ,ZCVG,mcon,buy2, & + zfdqb,dtdyn,dxlim,rmulacvg,dp,tem, & + alpha,DEN + integer :: inbu(im,km) + + !Parameters + gcvalmx = 0.1 + rmulacvg=10. + ZEPS7=1.E-11 + km1=km-1 + alpha=7000. + + !Initialization 2D + do k = 1,km + do i = 1,im + sigmaout(i,k)=0. + inbu(i,k)=0 + zform(i,k)=0. + enddo + enddo + + !Initialization 1D + do i=1,im + sigmab(i)=0. + sigmamax(i)=0.95 + termA(i)=0. + termB(i)=0. + termC(i)=0. + termD(i)=0. + zfdqa(i)=0. + mcons(i)=0. + enddo + + !Temporary Initialization output: + do i = 1,im + if(flag_deep)then + !ca_turb(i)=0. + ca_shal(i)=0. + endif + if(.not. flag_deep)then + ca_rad(i)=0. + convcount(i)=0. + ca1(i)=0. + endif + enddo + + !Initial computations, place maximum sigmain in sigmab + + do k=2,km + do i=1,im + if(flag_init .and. .not. flag_restart)then + if(cnvflg(i))then + sigmab(i)=0.03 + endif + else + if(cnvflg(i))then + !if(sigmain(i,k)<1.E-5)then + ! sigmain(i,k)=0. + !endif + if(sigmain(i,k)>sigmab(i))then + sigmab(i)=sigmain(i,k) + endif + endif + endif + enddo + enddo + + do i=1,im + if(sigmab(i) < 1.E-5)then !after advection + sigmab(i)=0. + endif + enddo + + !Initial computations, sigmamax + do i=1,im + sigmamax(i)=alpha/gdx(i) + sigmamax(i)=MIN(0.95,sigmamax(i)) + enddo + + !Initial computations, dynamic q-tendency + do k = 1,km + do i = 1,im + if(flag_init .and. .not.flag_restart)then + qadv(i,k)=0. + else + qadv(i,k)=(q(i,k) - qgrs_dsave(i,k))/delt + endif + enddo + enddo + + !compute termD "The vertical integral of the latent heat convergence is limited to the + !buoyant layers with positive moisture convergence (accumulated from the surface). + !Lowest level: + do i = 1,im + dp = 1000. * del(i,1) + mcons(i)=(hvap*(qadv(i,1)+tmf(i,1)+qmicro(i,1))*dp) + enddo + !Levels above: + do k = 2,km1 + do i = 1,im + dp = 1000. * del(i,k) + if(cnvflg(i))then + mcon = (hvap*(qadv(i,k)+tmf(i,k)+qmicro(i,k))*dp) + buy2 = termD(i)+mcon+mcons(i) +! Do the integral over buoyant layers with positive mcon acc from surface + if(k > kbcon1(i) .and. k < ktcon(i) .and. buy2 > 0.)then + inbu(i,k)=1 + endif + inbu(i,k-1)=MAX(inbu(i,k-1),inbu(i,k)) + termD(i) = termD(i) + float(inbu(i,k-1))*mcons(i) + mcons(i)=mcon + endif + enddo + enddo + + !termA + do k = 2,km1 + do i = 1,im + dp = 1000. * del(i,k) + if(cnvflg(i))then + tem=(sigmab(i)*zeta(i,k)*float(inbu(i,k))*dbyo1(i,k))*dp + termA(i)=termA(i)+tem + endif + enddo + enddo + + !termB + do k = 2,km1 + do i = 1,im + dp = 1000. * del(i,k) + if(cnvflg(i))then + tem=(dbyo1(i,k)*float(inbu(i,k)))*dp + termB(i)=termB(i)+tem + endif + enddo + enddo + + !termC + do k = 2,km1 + do i = 1,im + if(cnvflg(i))then + dp = 1000. * del(i,k) + zform(i,k)=-1.0*float(inbu(i,k))*(omega_u(i,k)*delt) + zfdqb=0.5*((zform(i,k)*zdqca(i,k))) + termC(i)=termC(i)+(float(inbu(i,k))* & + (zfdqb+zfdqa(i))*hvap*zeta(i,k)) + zfdqa(i)=zfdqb + endif + enddo + enddo + + !sigmab + do i = 1,im + if(cnvflg(i))then + + DEN=MIN(termC(i)+termB(i),1.E8) !1.E8 + !DEN=MAX(termC(i)+termB(i),1.E7) !1.E7 + + ZCVG=termD(i)*delt + + ZZ=MAX(0.0,SIGN(1.0,termA(i))) & + *MAX(0.0,SIGN(1.0,termB(i))) & + *MAX(0.0,SIGN(1.0,termC(i)-ZEPS7)) + + + ZCVG=MAX(0.0,ZCVG) + + if(flag_init)then + sigmab(i)=0.03 + else + sigmab(i)=(ZZ*(termA(i)+ZCVG))/(DEN+(1.0-ZZ)) + endif + + if(sigmab(i)>0.)then + sigmab(i)=MIN(sigmab(i),sigmamax(i)) + sigmab(i)=MAX(sigmab(i),0.01) + endif + + if(flag_deep)then + !ca_turb(i)=ZCVG + ca_shal(i)=termC(i) + else + ca_rad(i)=ZCVG + ca1(i)=termC(i) + endif + !ca3(i)=sigmab(i) + + endif!cnvflg + enddo + + do k=1,km + do i=1,im + if(cnvflg(i))then + sigmaout(i,k)=sigmab(i) + endif + enddo + enddo + + end subroutine progsigma_calc +!> @} +!! @} + + + From 842eae33944aafb5100ccb7660a04027a7452147 Mon Sep 17 00:00:00 2001 From: Lisa Bengtsson Date: Mon, 18 Apr 2022 15:05:40 +0000 Subject: [PATCH 03/20] Ensuring the moisture budget is correct via PBL, microphysics coupling --- physics/GFS_MP_generic.F90 | 26 +++++-- physics/GFS_MP_generic.meta | 23 ++++++ physics/progsigma_calc.f90 | 33 +------- physics/samfdeepcnv.f | 149 +++++++++++++++++++++++++++++++----- physics/samfdeepcnv.meta | 86 ++++++++++++++++++++- physics/satmedmfvdifq.F | 5 +- physics/satmedmfvdifq.meta | 8 ++ 7 files changed, 274 insertions(+), 56 deletions(-) diff --git a/physics/GFS_MP_generic.F90 b/physics/GFS_MP_generic.F90 index e106cb908..dbf2d15fa 100644 --- a/physics/GFS_MP_generic.F90 +++ b/physics/GFS_MP_generic.F90 @@ -88,14 +88,15 @@ end subroutine GFS_MP_generic_post_init !> @{ subroutine GFS_MP_generic_post_run( & im, levs, kdt, nrcm, nncl, ntcw, ntrac, imp_physics, imp_physics_gfdl, imp_physics_thompson, imp_physics_nssl, & - imp_physics_mg, imp_physics_fer_hires, cal_pre, cplflx, cplchm, con_g, rainmin, dtf, frain, rainc, & + imp_physics_mg, imp_physics_fer_hires, cal_pre, cplflx, cplchm, progsigma, con_g, rainmin, dtf, frain, rainc, & rain1, rann, xlat, xlon, gt0, gq0, prsl, prsi, phii, tsfc, ice, snow, graupel, save_t, save_q, rain0, ice0, snow0,& graupel0, del, rain, domr_diag, domzr_diag, domip_diag, doms_diag, tprcp, srflag, sr, cnvprcp, totprcp, totice, & totsnw, totgrp, cnvprcpb, totprcpb, toticeb, totsnwb, totgrpb, rain_cpl, rainc_cpl, snow_cpl, pwat, & drain_cpl, dsnow_cpl, lsm, lsm_ruc, lsm_noahmp, raincprv, rainncprv, iceprv, snowprv, & graupelprv, draincprv, drainncprv, diceprv, dsnowprv, dgraupelprv, dtp, dfi_radar_max_intervals, & - dtend, dtidx, index_of_temperature, index_of_process_mp,ldiag3d, qdiag3d, lssav, num_dfi_radar, fh_dfi_radar, & - index_of_process_dfi_radar, ix_dfi_radar, dfi_radar_tten, radar_tten_limits, fhour, errmsg, errflg) + dtend, dtidx, index_of_temperature, index_of_process_mp,ldiag3d, qdiag3d,dqdt_qmicro, lssav, num_dfi_radar, & + fh_dfi_radar,index_of_process_dfi_radar, ix_dfi_radar, dfi_radar_tten, radar_tten_limits, fhour, qgrs_dsave, & + errmsg, errflg) ! use machine, only: kind_phys @@ -104,7 +105,7 @@ subroutine GFS_MP_generic_post_run( integer, intent(in) :: im, levs, kdt, nrcm, nncl, ntcw, ntrac, num_dfi_radar, index_of_process_dfi_radar integer, intent(in) :: imp_physics, imp_physics_gfdl, imp_physics_thompson, imp_physics_mg, imp_physics_fer_hires integer, intent(in) :: imp_physics_nssl - logical, intent(in) :: cal_pre, lssav, ldiag3d, qdiag3d, cplflx, cplchm + logical, intent(in) :: cal_pre, lssav, ldiag3d, qdiag3d, cplflx, cplchm, progsigma integer, intent(in) :: index_of_temperature,index_of_process_mp integer :: dfi_radar_max_intervals @@ -148,7 +149,8 @@ subroutine GFS_MP_generic_post_run( real(kind=kind_phys), dimension(:), intent(inout) :: diceprv real(kind=kind_phys), dimension(:), intent(inout) :: dsnowprv real(kind=kind_phys), dimension(:), intent(inout) :: dgraupelprv - + real(kind=kind_phys), dimension(:,:), intent(out) :: dqdt_qmicro + real(kind=kind_phys), dimension(:,:), intent(out) :: qgrs_dsave real(kind=kind_phys), intent(in) :: dtp ! CCPP error handling @@ -420,6 +422,15 @@ subroutine GFS_MP_generic_post_run( endif if_tendency_diagnostics endif if_save_fields + !If prognostic updraft area fraction is used in saSAS + if(progsigma)then + do k=1,levs + do i=1,im + dqdt_qmicro(i,k)=(gq0(i,k,1)-save_q(i,k,1))/dtp + enddo + enddo + endif + if (cplflx .or. cplchm) then do i = 1, im dsnow_cpl(i)= max(zero, rain(i) * srflag(i)) @@ -455,6 +466,11 @@ subroutine GFS_MP_generic_post_run( pwat(i) = pwat(i) * onebg enddo + do k = 1, levs + do i=1, im + qgrs_dsave(i,k) = gq0(i,k,1) + enddo + enddo end subroutine GFS_MP_generic_post_run !> @} diff --git a/physics/GFS_MP_generic.meta b/physics/GFS_MP_generic.meta index 6177b1344..f8cb0acae 100644 --- a/physics/GFS_MP_generic.meta +++ b/physics/GFS_MP_generic.meta @@ -248,6 +248,13 @@ dimensions = () type = logical intent = in +[progsigma] + standard_name = flag_for_prognostic_sigma + long_name = flag for prognostic sigma + units = flag + dimensions = () + type = logical + intent = in [con_g] standard_name = gravitational_acceleration long_name = gravitational acceleration @@ -844,6 +851,22 @@ dimensions = () type = logical intent = in +[dqdt_qmicro] + standard_name = instantanious_moisture_tendency_due_to_microphysics + long_name = moisture tendency due to microphysics + units = kg kg-1 s-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = out +[qgrs_dsave] + standard_name = tracer_concentration_dsave + long_name = model layer mean tracer concentration dsave + units = kg kg-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = out [lssav] standard_name = flag_for_diagnostics long_name = logical flag for storing diagnostics diff --git a/physics/progsigma_calc.f90 b/physics/progsigma_calc.f90 index 378d43ef4..6772bc8d4 100644 --- a/physics/progsigma_calc.f90 +++ b/physics/progsigma_calc.f90 @@ -12,12 +12,11 @@ !!\section progsigma General Algorithm !> @{ - subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_deep, & + subroutine progsigma_calc (im,km,flag_init,flag_restart, & del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt, & qgrs_dsave,q,kbcon1,ktcon,cnvflg,gdx, & do_ca, ca_closure, ca_entr, ca_trigger, nthresh, ca_deep, & - ca_turb,ca_micro,ca_shal,ca_rad,convcount,ca1,ca2,ca3,ca4, & - sigmain,sigmaout,sigmab,errmsg,errflg) + ca_micro,sigmain,sigmaout,sigmab,errmsg,errflg) ! ! use machine, only : kind_phys @@ -31,12 +30,10 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_deep, & real, intent(in) :: qgrs_dsave(im,km), q(im,km),del(im,km), & qmicro(im,km),tmf(im,km),dbyo1(im,km),zdqca(im,km), & omega_u(im,km),zeta(im,km),gdx(im) - logical, intent(in) :: flag_init,flag_restart,flag_deep,cnvflg(im) + logical, intent(in) :: flag_init,flag_restart,cnvflg(im) real(kind=kind_phys), intent(in) :: nthresh real(kind=kind_phys), intent(in) :: ca_deep(im) - real(kind=kind_phys), intent(out):: ca_turb(im), & - ca_micro(im),ca_rad(im),ca_shal(im),convcount(im),ca1(im), & - ca2(im),ca3(im),ca4(im) + real(kind=kind_phys), intent(out):: ca_micro(im) logical, intent(in) :: do_ca,ca_closure,ca_entr,ca_trigger real(kind=kind_phys), intent(in) :: sigmain(im,km) @@ -87,19 +84,6 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_deep, & mcons(i)=0. enddo - !Temporary Initialization output: - do i = 1,im - if(flag_deep)then - !ca_turb(i)=0. - ca_shal(i)=0. - endif - if(.not. flag_deep)then - ca_rad(i)=0. - convcount(i)=0. - ca1(i)=0. - endif - enddo - !Initial computations, place maximum sigmain in sigmab do k=2,km @@ -232,15 +216,6 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_deep, & sigmab(i)=MAX(sigmab(i),0.01) endif - if(flag_deep)then - !ca_turb(i)=ZCVG - ca_shal(i)=termC(i) - else - ca_rad(i)=ZCVG - ca1(i)=termC(i) - endif - !ca3(i)=sigmab(i) - endif!cnvflg enddo diff --git a/physics/samfdeepcnv.f b/physics/samfdeepcnv.f index ea92fda7f..7be31a04a 100644 --- a/physics/samfdeepcnv.f +++ b/physics/samfdeepcnv.f @@ -75,17 +75,18 @@ end subroutine samfdeepcnv_finalize !! !! \section samfdeep_detailed GFS samfdeepcnv Detailed Algorithm !! @{ - subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & + subroutine samfdeepcnv_run (im,km,first_time_step,restart, & + & tmf,qmicro,itc,ntc,cliq,cp,cvap, & & eps,epsm1,fv,grav,hvap,rd,rv, & & t0c,delt,ntk,ntr,delp, & - & prslp,psp,phil,qtr,q1,t1,u1,v1,fscav,hwrf_samfdeep, & - & cldwrk,rn,kbot,ktop,kcnv,islimsk,garea, & - & dot,ncloud,hpbl,ud_mf,dd_mf,dt_mf,cnvw,cnvc, & + & prslp,psp,phil,qtr,qgrs_dsave,q,q1,t1,u1,v1,fscav, & + & hwrf_samfdeep,progsigma,wclosureflg,cldwrk,rn,kbot,ktop,kcnv, & + & islimsk,garea,dot,ncloud,hpbl,ud_mf,dd_mf,dt_mf,cnvw,cnvc, & & QLCN, QICN, w_upi, cf_upi, CNV_MFD, & & CNV_DQLDT,CLCN,CNV_FICE,CNV_NDROP,CNV_NICE,mp_phys,mp_phys_mg,& & clam,c0s,c1,betal,betas,evef,pgcon,asolfac, & & do_ca, ca_closure, ca_entr, ca_trigger, nthresh, ca_deep, & - & rainevap, & + & rainevap, sigmain, sigmaout, ca_micro, & & errmsg,errflg) ! use machine , only : kind_phys @@ -101,10 +102,14 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & real(kind=kind_phys), intent(in) :: psp(:), delp(:,:), & & prslp(:,:), garea(:), hpbl(:), dot(:,:), phil(:,:) real(kind=kind_phys), dimension(:), intent(in) :: fscav - logical, intent(in) :: hwrf_samfdeep + logical, intent(in) :: first_time_step,restart,hwrf_samfdeep, & + & progsigma, wclosureflg real(kind=kind_phys), intent(in) :: nthresh real(kind=kind_phys), intent(in) :: ca_deep(:) - real(kind=kind_phys), intent(out) :: rainevap(:) + real(kind=kind_phys), intent(in) :: sigmain(:,:),qmicro(:,:), & + & tmf(:,:),q(:,:), qgrs_dsave(:,:) + real(kind=kind_phys), intent(out) :: rainevap(:),ca_micro(:) + real(kind=kind_phys), intent(out) :: sigmaout(:,:) logical, intent(in) :: do_ca,ca_closure,ca_entr,ca_trigger integer, intent(inout) :: kcnv(:) @@ -208,6 +213,10 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & ! & bb1, bb2 & bb1, bb2, wucb ! +! parameters for prognostic sigma closure + real(kind=kind_phys) omega_u(im,km),zdqca(im,km),qlks(im,km), + & omegac(im),zeta(im,km),dbyo1(im,km),sigmab(im) + c physical parameters ! parameter(grav=grav,asolfac=0.958) ! parameter(elocp=hvap/cp,el2orc=hvap*hvap/(rv*cp)) @@ -368,6 +377,7 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & vshear(i) = 0. advfac(i) = 0. rainevap(i) = 0. + omegac(i)=0. gdx(i) = sqrt(garea(i)) enddo @@ -570,6 +580,11 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & buo(i,k) = 0. drag(i,k) = 0. cnvwt(i,k)= 0. + dbyo1(i,k)=0. + zdqca(i,k)=0. + qlks(i,k)=0. + omega_u(i,k)=0. + zeta(i,k)=1.0 endif enddo enddo @@ -1497,6 +1512,7 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & pwavo(i) = pwavo(i) + pwo(i,k) ! cnvwt(i,k) = (etah*qlk + pwo(i,k)) * grav / dp cnvwt(i,k) = etah * qlk * grav / dp + qlks(i,k)=qlk endif ! ! compute buoyancy and drag for updraft velocity @@ -1569,6 +1585,7 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & dz1 = zo(i,k+1) - zo(i,k) ! aa1(i) = aa1(i) + buo(i,k) * dz1 * eta(i,k) aa1(i) = aa1(i) + buo(i,k) * dz1 + dbyo1(i,k) = hcko(i,k) - heso(i,k) endif endif enddo @@ -1669,6 +1686,7 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & pwavo(i) = pwavo(i) + pwo(i,k) ! cnvwt(i,k) = (etah*qlk + pwo(i,k)) * grav / dp cnvwt(i,k) = etah * qlk * grav / dp + qlks(i,k)=qlk endif endif endif @@ -1710,6 +1728,20 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo + + if(progsigma)then + do k = 2, km1 + do i = 1, im + if (cnvflg(i)) then + if(k > kbcon1(i) .and. k < ktcon(i)) then + rho = po(i,k)*100. / (rd * to(i,k)) + omega_u(i,k)=-1.0*sqrt(wu2(i,k))*rho*grav + omega_u(i,k)=MAX(omega_u(i,k),-80.) + endif + endif + enddo + enddo + endif ! ! compute updraft velocity average over the whole cumulus ! @@ -1742,6 +1774,54 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & endif enddo c + +!> - Calculate the mean updraft velocity within the cloud (wc),cast in pressure coordinates. + if(progsigma)then + + do i = 1, im + omegac(i) = 0. + sumx(i) = 0. + enddo + do k = 2, km1 + do i = 1, im + if (cnvflg(i)) then + if(k > kbcon1(i) .and. k < ktcon(i)) then + dp = 1000. * del(i,k) + tem = 0.5 * (omega_u(i,k) + omega_u(i,k-1)) + omegac(i) = omegac(i) + tem * dp + sumx(i) = sumx(i) + dp + endif + endif + enddo + enddo + do i = 1, im + if(cnvflg(i)) then + if(sumx(i) == 0.) then + cnvflg(i)=.false. + else + omegac(i) = omegac(i) / sumx(i) + endif + val = -1.2 + if (omegac(i) > val) cnvflg(i)=.false. + endif + enddo + +!> - Calculate the xi term in Bengtsson et al. 2022 (equation 8) + do k = 2, km1 + do i = 1, im + if (cnvflg(i)) then + if(k >= kbcon1(i) .and. k < ktcon(i)) then + zeta(i,k)=eta(i,k)*(omegac(i)/omega_u(i,k)) + zeta(i,k)=MAX(0.,zeta(i,k)) + zeta(i,k)=MIN(1.,zeta(i,k)) + endif + endif + enddo + enddo + + + endif !if progsigma + c exchange ktcon with ktcon1 c !> - Swap the indices of the convective cloud top (ktcon) and the overshooting convection top (ktcon1) to use the same cloud top level in the calculations of \f$A^+\f$ and \f$A^*\f$. @@ -1773,11 +1853,26 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & if(dq > 0.) then qlko_ktcon(i) = dq qcko(i,k) = qrch + qlks(i,k) = qlko_ktcon(i) endif endif enddo endif c + +c store term needed for "termC" in prognostic area fraction closure + do k = 2, km1 + do i = 1, im + dp = 1000. * del(i,k) + if (cnvflg(i)) then + if(k > kbcon(i) .and. k < ktcon(i)) then + zdqca(i,k)=((qlks(i,k)-qlks(i,k-1)) + + & pwo(i,k)+dellal(i,k)) + endif + endif + enddo + enddo + ccccc if(lat.==.latd.and.lon.==.lond.and.cnvflg(i)) then ccccc print *, ' aa1(i) before dwndrft =', aa1(i) ccccc endif @@ -2375,6 +2470,14 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & asqecflg(i) = .false. endif enddo + +!> - If wclosureflg is true, then quasi-equilibrium closure of Arakawa-Schubert is not used any longer, regardless of resolution + if(wclosureflg)then + do i = 1, im + asqecflg(i) = .false. + enddo + endif + ! !> - If grid size is larger than the threshold value (i.e., asqecflg=.true.), the quasi-equilibrium assumption is used to obtain the cloud base mass flux. To begin with, calculate the change in the temperature and moisture profiles per unit cloud base mass flux. do k = 1, km @@ -2784,13 +2887,27 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & advfac(i) = min(advfac(i), 1.) endif enddo + +!> - From Bengtsson et al. (2022) Prognostic closure scheme, equation 8, compute updraft area fraction based on a moisture budget + if(progsigma)then + call progsigma_calc(im,km,first_time_step,restart, + & del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt, + & qgrs_dsave,q,kbcon1,ktcon,cnvflg,gdx, + & do_ca, ca_closure, ca_entr, ca_trigger, nthresh, ca_deep, + & ca_micro,sigmain,sigmaout,sigmab,errmsg,errflg) + 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 for the grid sizes where the quasi-equilibrium assumption of Arakawa-Schubert is not valid any longer. !! As discussed in Han et al. (2017) \cite han_et_al_2017 , when dtconv is larger than tauadv, the convective mixing is not fully conducted before the cumulus cloud is advected out of the grid cell. In this case, therefore, the cloud base mass flux is further reduced in proportion to the ratio of tauadv to dtconv. do i= 1, im if(cnvflg(i) .and. .not.asqecflg(i)) then k = kbcon(i) rho = po(i,k)*100. / (rd*to(i,k)) - xmb(i) = advfac(i)*betaw*rho*wc(i) + if(progsigma)then + xmb(i) = sigmab(i)*((-1.0*omegac(i))/grav) + else + xmb(i) = advfac(i)*betaw*rho*wc(i) + endif endif enddo !> - For the cases where the quasi-equilibrium assumption of Arakawa-Schubert is valid, first calculate the large scale destabilization as in equation 5 of Pan and Wu (1995) \cite pan_and_wu_1995 : @@ -2859,7 +2976,11 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & do i = 1, im if(cnvflg(i)) then if (gdx(i) < dxcrtuf) then - scaldfunc(i) = (1.-sigmagfm(i)) * (1.-sigmagfm(i)) + if(progsigma)then + scaldfunc(i)=(1.-sigmab(i))*(1.-sigmab(i)) + else + scaldfunc(i) = (1.-sigmagfm(i)) * (1.-sigmagfm(i)) + endif scaldfunc(i) = max(min(scaldfunc(i), 1.0), 0.) else scaldfunc(i) = 1.0 @@ -2869,16 +2990,6 @@ subroutine samfdeepcnv_run (im,km,itc,ntc,cliq,cp,cvap, & endif enddo ! - if (do_ca .and. ca_closure)then - do i = 1, im - if(cnvflg(i)) then - if (ca_deep(i) > nthresh) then - xmb(i) = xmb(i)*1.25 - endif - endif - enddo - endif - !> - Transport aerosols if present ! ! if (do_aerosols) diff --git a/physics/samfdeepcnv.meta b/physics/samfdeepcnv.meta index baf01fb8e..3eb330551 100644 --- a/physics/samfdeepcnv.meta +++ b/physics/samfdeepcnv.meta @@ -1,7 +1,7 @@ [ccpp-table-properties] name = samfdeepcnv type = scheme - dependencies = funcphys.f90,machine.F,samfaerosols.F + dependencies = funcphys.f90,machine.F,samfaerosols.F,progsigma_calc.f90 ######################################################################## [ccpp-arg-table] @@ -55,6 +55,36 @@ dimensions = () type = integer intent = in +[first_time_step] + standard_name = flag_for_first_timestep + long_name = flag for first time step for time integration loop (cold/warmstart) + units = flag + dimensions = () + type = logical + intent = in +[restart] + standard_name = flag_for_restart + long_name = flag for restart (warmstart) or coldstart + units = flag + dimensions = () + type = logical + intent = in +[tmf] + standard_name = turbulence_moisture_flux + long_name = turbulence_moisture_flux + units = kg kg-1 s-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in +[qmicro] + standard_name = instantanious_moisture_tendency_due_to_microphysics + long_name = moisture tendency due to microphysics + units = kg kg-1 s-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in [itc] standard_name = index_of_first_chemical_tracer_for_convection long_name = index of first chemical tracer transported/scavenged by convection @@ -219,6 +249,22 @@ type = real kind = kind_phys intent = inout +[qgrs_dsave] + standard_name = tracer_concentration_dsave + long_name = model layer mean tracer concentration dsave + units = kg kg-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in +[q] + standard_name = specific_humidity + long_name = water vapor specific humidity + units = kg kg-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in [q1] standard_name = specific_humidity_of_new_state long_name = updated vapor specific humidity @@ -266,6 +312,28 @@ dimensions = () type = logical intent = in +[wclosureflg] + standard_name = flag_for_wclosure + long_name = flag for vertical velocity closure + units = flag + dimensions = () + type = logical + intent = in +[progsigma] + standard_name = flag_for_prognostic_sigma + long_name = flag for prognostic sigma + units = flag + dimensions = () + type = logical + intent = in +[ca_micro] + standard_name = output_prognostic_sigma_two + long_name = output of prognostic area fraction two + units = frac + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out [cldwrk] standard_name = cloud_work_function long_name = cloud work function @@ -381,6 +449,22 @@ type = real kind = kind_phys intent = inout +[sigmain] + standard_name = updraft_area_fraction + long_name = convective updraft area fraction + units = frac + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = inout +[sigmaout] + standard_name = updraft_area_fraction_updated_by_physics + long_name = convective updraft area fraction updated by physics + units = frac + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = out [qlcn] standard_name = mass_fraction_of_convective_cloud_liquid_water long_name = mass fraction of convective cloud liquid water diff --git a/physics/satmedmfvdifq.F b/physics/satmedmfvdifq.F index eb2b7ad1c..b0367b4d4 100644 --- a/physics/satmedmfvdifq.F +++ b/physics/satmedmfvdifq.F @@ -69,7 +69,7 @@ end subroutine satmedmfvdifq_finalize !! @{ subroutine satmedmfvdifq_run(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, & + & dv,du,tdt,rtg,tmf,u1,v1,t1,q1,swh,hlw,xmu,garea,zvfun, & & psk,rbsoil,zorl,u10m,v10m,fm,fh, & & tsea,heat,evap,stress,spd1,kpbl, & & prsi,del,prsl,prslk,phii,phil,delt, & @@ -121,7 +121,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & real(kind=kind_phys), intent(out) :: & & dusfc(:), dvsfc(:), & & dtsfc(:), dqsfc(:), & - & hpbl(:) + & hpbl(:), tmf(:,:) real(kind=kind_phys), intent(out) :: & & dkt(:,:), dku(:,:) ! @@ -2114,6 +2114,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & qtend = (f2(i,k)-q1(i,k,1))*rdt tdt(i,k) = tdt(i,k)+ttend rtg(i,k,1) = rtg(i,k,1)+qtend + tmf(i,k) = qtend ! dtsfc(i) = dtsfc(i)+cont*del(i,k)*ttend ! dqsfc(i) = dqsfc(i)+conq*del(i,k)*qtend enddo diff --git a/physics/satmedmfvdifq.meta b/physics/satmedmfvdifq.meta index db89f488d..211bbaec6 100644 --- a/physics/satmedmfvdifq.meta +++ b/physics/satmedmfvdifq.meta @@ -201,6 +201,14 @@ type = real kind = kind_phys intent = inout +[tmf] + standard_name = turbulence_moisture_flux + long_name = turbulence_moisture_flux + units = kg kg-1 s-1 + dimensions = (horizontal_dimension,vertical_layer_dimension) + type = real + kind = kind_phys + intent = out [u1] standard_name = x_wind long_name = x component of layer wind From 4f84ed749af88d8d833d78ec352c2a02dfb52589 Mon Sep 17 00:00:00 2001 From: Lisa Bengtsson Date: Tue, 19 Apr 2022 19:22:53 +0000 Subject: [PATCH 04/20] add shallow convection closure updates, add ntsigma in generic files --- physics/GFS_DCNV_generic.F90 | 14 ++-- physics/GFS_DCNV_generic.meta | 14 ++++ physics/GFS_MP_generic.meta | 2 +- physics/GFS_SCNV_generic.F90 | 12 +-- physics/GFS_SCNV_generic.meta | 14 ++++ physics/progsigma_calc.f90 | 73 +++++++++--------- physics/samfdeepcnv.f | 31 ++++---- physics/samfdeepcnv.meta | 10 +-- physics/samfshalcnv.f | 139 +++++++++++++++++++++++++++++++--- physics/samfshalcnv.meta | 79 ++++++++++++++++++- physics/satmedmfvdifq.F | 5 +- physics/satmedmfvdifq.meta | 6 +- 12 files changed, 315 insertions(+), 84 deletions(-) diff --git a/physics/GFS_DCNV_generic.F90 b/physics/GFS_DCNV_generic.F90 index a9e0ba7e0..07defcba1 100644 --- a/physics/GFS_DCNV_generic.F90 +++ b/physics/GFS_DCNV_generic.F90 @@ -18,8 +18,8 @@ end subroutine GFS_DCNV_generic_pre_finalize subroutine GFS_DCNV_generic_pre_run (im, levs, ldiag3d, qdiag3d, do_cnvgwd, cplchm, & gu0, gv0, gt0, gq0, nsamftrac, ntqv, & save_u, save_v, save_t, save_q, clw, & - ntcw,ntiw,ntclamt,ntrw,ntsw,ntrnc,ntsnc,ntgl, & - ntgnc, nthl, nthnc, nthv, ntgv, & + ntcw,ntiw,ntclamt,ntsigma,ntrw,ntsw,ntrnc,ntsnc,& + ntgl,ntgnc, nthl, nthnc, nthv, ntgv, & cscnv, satmedmf, trans_trac, ras, ntrac, & dtidx, index_of_process_dcnv, errmsg, errflg) @@ -28,7 +28,7 @@ subroutine GFS_DCNV_generic_pre_run (im, levs, ldiag3d, qdiag3d, do_cnvgwd, cplc implicit none integer, intent(in) :: im, levs, nsamftrac, ntqv, index_of_process_dcnv, dtidx(:,:), & - ntcw,ntiw,ntclamt,ntrw,ntsw,ntrnc,ntsnc,ntgl,ntrac,ntgnc,nthl,nthnc,nthv,ntgv + ntcw,ntiw,ntclamt,ntrw,ntsw,ntrnc,ntsnc,ntgl,ntrac,ntgnc,nthl,nthnc,nthv,ntgv,ntsigma logical, intent(in) :: ldiag3d, qdiag3d, do_cnvgwd, cplchm real(kind=kind_phys), dimension(:,:), intent(in) :: gu0 real(kind=kind_phys), dimension(:,:), intent(in) :: gv0 @@ -74,7 +74,7 @@ subroutine GFS_DCNV_generic_pre_run (im, levs, ldiag3d, qdiag3d, do_cnvgwd, cplc n /= ntrw .and. n /= ntsw .and. n /= ntrnc .and. & n /= ntsnc .and. n /= ntgl .and. n /= ntgnc .and. & n /= nthl .and. n /= nthnc .and. n /= nthv .and. & - n /= ntgv ) then + n /= ntgv .and. n /= ntsigma) then tracers = tracers + 1 if(dtidx(100+n,index_of_process_dcnv)>0) then save_q(:,:,n) = clw(:,:,tracers) @@ -114,7 +114,7 @@ subroutine GFS_DCNV_generic_post_run (im, levs, lssav, ldiag3d, qdiag3d, ras, & rainc, cldwrk, upd_mf, dwn_mf, det_mf, dtend, dtidx, index_of_process_dcnv, & index_of_temperature, index_of_x_wind, index_of_y_wind, ntqv, gq0, save_q, & cnvw, cnvc, cnvw_phy_f3d, cnvc_phy_f3d, flag_for_dcnv_generic_tend, & - ntcw,ntiw,ntclamt,ntrw,ntsw,ntrnc,ntsnc,ntgl, & + ntcw,ntiw,ntclamt,ntsigma,ntrw,ntsw,ntrnc,ntsnc,ntgl, & ntgnc, nthl, nthnc, nthv, ntgv, ntrac,clw, & satmedmf, trans_trac, errmsg, errflg) @@ -145,7 +145,7 @@ subroutine GFS_DCNV_generic_post_run (im, levs, lssav, ldiag3d, qdiag3d, ras, & integer, intent(in) :: dtidx(:,:), index_of_process_dcnv, index_of_temperature, & index_of_x_wind, index_of_y_wind, ntqv integer, intent(in) :: ntcw,ntiw,ntclamt,ntrw,ntsw,ntrnc,ntsnc,ntgl, & - ntgnc, nthl, nthnc, nthv, ntgv, ntrac + ntgnc, nthl, nthnc, nthv, ntgv, ntsigma, ntrac real(kind=kind_phys), dimension(:,:,:), intent(in) :: clw @@ -212,7 +212,7 @@ subroutine GFS_DCNV_generic_post_run (im, levs, lssav, ldiag3d, qdiag3d, ras, & n /= ntrw .and. n /= ntsw .and. n /= ntrnc .and. & n /= ntsnc .and. n /= ntgl .and. n /= ntgnc .and. & n /= nthl .and. n /= nthnc .and. n /= nthv .and. & - n /= ntgv ) then + n /= ntgv .and. n /= ntsigma) then tracers = tracers + 1 idtend = dtidx(100+n,index_of_process_dcnv) if(idtend>0) then diff --git a/physics/GFS_DCNV_generic.meta b/physics/GFS_DCNV_generic.meta index e15acaf1c..f095259d4 100644 --- a/physics/GFS_DCNV_generic.meta +++ b/physics/GFS_DCNV_generic.meta @@ -190,6 +190,13 @@ dimensions = () type = integer intent = in +[ntsigma] + standard_name = index_of_updraft_area_fraction_in_tracer_concentration_array + long_name = tracer index of updraft_area_fraction + units = index + dimensions = () + type = integer + intent = in [ntrw] standard_name = index_of_rain_mixing_ratio_in_tracer_concentration_array long_name = tracer index for rain water @@ -670,6 +677,13 @@ dimensions = () type = integer intent = in +[ntsigma] + standard_name = index_of_updraft_area_fraction_in_tracer_concentration_array + long_name = tracer index of updraft_area_fraction + units = index + dimensions = () + type = integer + intent = in [ntrw] standard_name = index_of_rain_mixing_ratio_in_tracer_concentration_array long_name = tracer index for rain water diff --git a/physics/GFS_MP_generic.meta b/physics/GFS_MP_generic.meta index f8cb0acae..763cad85a 100644 --- a/physics/GFS_MP_generic.meta +++ b/physics/GFS_MP_generic.meta @@ -852,7 +852,7 @@ type = logical intent = in [dqdt_qmicro] - standard_name = instantanious_moisture_tendency_due_to_microphysics + standard_name = instantaneous_moisture_tendency_due_to_microphysics long_name = moisture tendency due to microphysics units = kg kg-1 s-1 dimensions = (horizontal_loop_extent,vertical_layer_dimension) diff --git a/physics/GFS_SCNV_generic.F90 b/physics/GFS_SCNV_generic.F90 index 58447f6bf..cbef02bb0 100644 --- a/physics/GFS_SCNV_generic.F90 +++ b/physics/GFS_SCNV_generic.F90 @@ -17,14 +17,14 @@ end subroutine GFS_SCNV_generic_pre_finalize subroutine GFS_SCNV_generic_pre_run (im, levs, ldiag3d, qdiag3d, gu0, gv0, gt0, gq0, & save_u, save_v, save_t, save_q, ntqv, nsamftrac, flag_for_scnv_generic_tend, & dtidx, index_of_process_scnv, ntcw,ntiw,ntclamt,ntrw,ntsw,ntrnc,ntsnc,ntgl,ntgnc, & - cscnv, satmedmf, trans_trac, ras, ntrac, clw, errmsg, errflg) + ntsigma, cscnv, satmedmf, trans_trac, ras, ntrac, clw, errmsg, errflg) use machine, only: kind_phys implicit none integer, intent(in) :: im, levs, ntqv, nsamftrac, index_of_process_scnv, dtidx(:,:) - integer, intent(in) :: ntcw,ntiw,ntclamt,ntrw,ntsw,ntrnc,ntsnc,ntgl,ntgnc,ntrac + integer, intent(in) :: ntcw,ntiw,ntclamt,ntsigma,ntrw,ntsw,ntrnc,ntsnc,ntgl,ntgnc,ntrac logical, intent(in) :: ldiag3d, qdiag3d, flag_for_scnv_generic_tend real(kind=kind_phys), dimension(:,:), intent(in) :: gu0, gv0, gt0 real(kind=kind_phys), dimension(:,:,:), intent(in) :: gq0 @@ -55,7 +55,7 @@ subroutine GFS_SCNV_generic_pre_run (im, levs, ldiag3d, qdiag3d, gu0, gv0, gt0, do n=2,ntrac if ( n /= ntcw .and. n /= ntiw .and. n /= ntclamt .and. & n /= ntrw .and. n /= ntsw .and. n /= ntrnc .and. & - n /= ntsnc .and. n /= ntgl .and. n /= ntgnc) then + n /= ntsnc .and. n /= ntgl .and. n /= ntgnc .and. n /= ntsigma) then tracers = tracers + 1 if(dtidx(100+n,index_of_process_scnv)>0) then save_q(:,:,n) = clw(:,:,tracers) @@ -97,7 +97,7 @@ subroutine GFS_SCNV_generic_post_run (im, levs, nn, lssav, ldiag3d, qdiag3d, & rainc, cnvprcp, cnvprcpb, cnvw_phy_f3d, cnvc_phy_f3d, & dtend, dtidx, index_of_temperature, index_of_x_wind, index_of_y_wind, & index_of_process_scnv, ntqv, flag_for_scnv_generic_tend, & - ntcw,ntiw,ntclamt,ntrw,ntsw,ntrnc,ntsnc,ntgl,ntgnc, & + ntcw,ntiw,ntclamt,ntsigma,ntrw,ntsw,ntrnc,ntsnc,ntgl,ntgnc, & imfshalcnv, imfshalcnv_sas, imfshalcnv_samf, ntrac, & cscnv, satmedmf, trans_trac, ras, errmsg, errflg) @@ -106,7 +106,7 @@ subroutine GFS_SCNV_generic_post_run (im, levs, nn, lssav, ldiag3d, qdiag3d, & implicit none integer, intent(in) :: im, levs, nn, ntqv, nsamftrac - integer, intent(in) :: ntcw,ntiw,ntclamt,ntrw,ntsw,ntrnc,ntsnc,ntgl,ntgnc,ntrac + integer, intent(in) :: ntcw,ntiw,ntclamt,ntsigma,ntrw,ntsw,ntrnc,ntsnc,ntgl,ntgnc,ntrac logical, intent(in) :: lssav, ldiag3d, qdiag3d, flag_for_scnv_generic_tend real(kind=kind_phys), intent(in) :: frain real(kind=kind_phys), dimension(:,:), intent(in) :: gu0, gv0, gt0 @@ -186,7 +186,7 @@ subroutine GFS_SCNV_generic_post_run (im, levs, nn, lssav, ldiag3d, qdiag3d, & do n=2,ntrac if ( n /= ntcw .and. n /= ntiw .and. n /= ntclamt .and. & n /= ntrw .and. n /= ntsw .and. n /= ntrnc .and. & - n /= ntsnc .and. n /= ntgl .and. n /= ntgnc) then + n /= ntsnc .and. n /= ntgl .and. n /= ntgnc .and. n/= ntsigma) then tracers = tracers + 1 idtend = dtidx(100+n,index_of_process_scnv) if(idtend>0) then diff --git a/physics/GFS_SCNV_generic.meta b/physics/GFS_SCNV_generic.meta index 5cbda127c..4fd189948 100644 --- a/physics/GFS_SCNV_generic.meta +++ b/physics/GFS_SCNV_generic.meta @@ -183,6 +183,13 @@ dimensions = () type = integer intent = in +[ntsigma] + standard_name = index_of_updraft_area_fraction_in_tracer_concentration_array + long_name = tracer index of updraft_area_fraction + units = index + dimensions = () + type = integer + intent = in [ntrw] standard_name = index_of_rain_mixing_ratio_in_tracer_concentration_array long_name = tracer index for rain water @@ -614,6 +621,13 @@ dimensions = () type = integer intent = in +[ntsigma] + standard_name = index_of_updraft_area_fraction_in_tracer_concentration_array + long_name = tracer index of updraft_area_fraction + units = index + dimensions = () + type = integer + intent = in [ntrw] standard_name = index_of_rain_mixing_ratio_in_tracer_concentration_array long_name = tracer index for rain water diff --git a/physics/progsigma_calc.f90 b/physics/progsigma_calc.f90 index 6772bc8d4..ca05f6778 100644 --- a/physics/progsigma_calc.f90 +++ b/physics/progsigma_calc.f90 @@ -1,7 +1,8 @@ !>\file progsigma !! This file contains the subroutine that calculates the prognostic !! updraft area fraction that is used for closure computations in -!! saSAS deep and shallow convection. +!! saSAS deep and shallow convection, based on a moisture budget +!! as described in Bengtsson et al. 2022. !>\ingroup samfdeepcnv !! This subroutine computes a prognostic updraft area fraction @@ -13,9 +14,8 @@ !> @{ subroutine progsigma_calc (im,km,flag_init,flag_restart, & - del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt, & - qgrs_dsave,q,kbcon1,ktcon,cnvflg,gdx, & - do_ca, ca_closure, ca_entr, ca_trigger, nthresh, ca_deep, & + flag_shallow,del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap, & + delt,qgrs_dsave,q,kbcon1,ktcon,cnvflg,gdx, & ca_micro,sigmain,sigmaout,sigmab,errmsg,errflg) ! ! @@ -30,12 +30,8 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & real, intent(in) :: qgrs_dsave(im,km), q(im,km),del(im,km), & qmicro(im,km),tmf(im,km),dbyo1(im,km),zdqca(im,km), & omega_u(im,km),zeta(im,km),gdx(im) - logical, intent(in) :: flag_init,flag_restart,cnvflg(im) - real(kind=kind_phys), intent(in) :: nthresh - real(kind=kind_phys), intent(in) :: ca_deep(im) + logical, intent(in) :: flag_init,flag_restart,cnvflg(im),flag_shallow real(kind=kind_phys), intent(out):: ca_micro(im) - logical, intent(in) :: do_ca,ca_closure,ca_entr,ca_trigger - real(kind=kind_phys), intent(in) :: sigmain(im,km) ! intent out @@ -47,28 +43,29 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & ! Local variables integer :: i,k,km1 real(kind=kind_phys) :: termA(im),termB(im),termC(im),termD(im), & - mcons(im),zfdqa(im),zform(im,km), & + mcons(im),fdqa(im),form(im,km), & qadv(im,km),sigmamax(im) - real(kind=kind_phys) :: gcvalmx,ZEPS7,ZZ,ZCVG,mcon,buy2, & - zfdqb,dtdyn,dxlim,rmulacvg,dp,tem, & - alpha,DEN + real(kind=kind_phys) :: gcvalmx,epsilon,ZZ,cvg,mcon,buy2, & + fdqb,dtdyn,dxlim,rmulacvg,dp,tem, & + alpha,DEN,betascu integer :: inbu(im,km) !Parameters gcvalmx = 0.1 rmulacvg=10. - ZEPS7=1.E-11 + epsilon=1.E-11 km1=km-1 alpha=7000. + betascu = 3.0 !Initialization 2D do k = 1,km do i = 1,im sigmaout(i,k)=0. inbu(i,k)=0 - zform(i,k)=0. + form(i,k)=0. enddo enddo @@ -80,8 +77,9 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & termB(i)=0. termC(i)=0. termD(i)=0. - zfdqa(i)=0. + fdqa(i)=0. mcons(i)=0. + ca_micro(i)=0. enddo !Initial computations, place maximum sigmain in sigmab @@ -94,9 +92,6 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & endif else if(cnvflg(i))then - !if(sigmain(i,k)<1.E-5)then - ! sigmain(i,k)=0. - !endif if(sigmain(i,k)>sigmab(i))then sigmab(i)=sigmain(i,k) endif @@ -107,7 +102,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & do i=1,im if(sigmab(i) < 1.E-5)then !after advection - sigmab(i)=0. + sigmab(i)=0. endif enddo @@ -180,11 +175,11 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & do i = 1,im if(cnvflg(i))then dp = 1000. * del(i,k) - zform(i,k)=-1.0*float(inbu(i,k))*(omega_u(i,k)*delt) - zfdqb=0.5*((zform(i,k)*zdqca(i,k))) + form(i,k)=-1.0*float(inbu(i,k))*(omega_u(i,k)*delt) + fdqb=0.5*((form(i,k)*zdqca(i,k))) termC(i)=termC(i)+(float(inbu(i,k))* & - (zfdqb+zfdqa(i))*hvap*zeta(i,k)) - zfdqa(i)=zfdqb + (fdqb+fdqa(i))*hvap*zeta(i,k)) + fdqa(i)=fdqb endif enddo enddo @@ -193,29 +188,26 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & do i = 1,im if(cnvflg(i))then - DEN=MIN(termC(i)+termB(i),1.E8) !1.E8 - !DEN=MAX(termC(i)+termB(i),1.E7) !1.E7 - - ZCVG=termD(i)*delt - + DEN=MIN(termC(i)+termB(i),1.E8) + cvg=termD(i)*delt ZZ=MAX(0.0,SIGN(1.0,termA(i))) & *MAX(0.0,SIGN(1.0,termB(i))) & - *MAX(0.0,SIGN(1.0,termC(i)-ZEPS7)) + *MAX(0.0,SIGN(1.0,termC(i)-epsilon)) - ZCVG=MAX(0.0,ZCVG) + cvg=MAX(0.0,cvg) - if(flag_init)then + if(flag_init .and. .not. flag_restart)then sigmab(i)=0.03 else - sigmab(i)=(ZZ*(termA(i)+ZCVG))/(DEN+(1.0-ZZ)) + sigmab(i)=(ZZ*(termA(i)+cvg))/(DEN+(1.0-ZZ)) endif if(sigmab(i)>0.)then sigmab(i)=MIN(sigmab(i),sigmamax(i)) sigmab(i)=MAX(sigmab(i),0.01) endif - + ca_micro(i)=sigmab(i) endif!cnvflg enddo @@ -226,7 +218,20 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & endif enddo enddo + + !Since updraft velocity is much lower in shallow cu region, termC becomes small in shallow cu application, thus the area fraction + !in this regime becomes too large compared with the deep cu region. To address this simply apply a scaling factor for shallow cu + !before computing the massflux to reduce the total strength of the SC MF: + if(flag_shallow)then + do i= 1, im + if(cnvflg(i)) then + sigmab(i)=sigmab(i)/betascu + endif + enddo + endif + + end subroutine progsigma_calc !> @} !! @} diff --git a/physics/samfdeepcnv.f b/physics/samfdeepcnv.f index 7be31a04a..35aea0eb1 100644 --- a/physics/samfdeepcnv.f +++ b/physics/samfdeepcnv.f @@ -216,7 +216,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & ! parameters for prognostic sigma closure real(kind=kind_phys) omega_u(im,km),zdqca(im,km),qlks(im,km), & omegac(im),zeta(im,km),dbyo1(im,km),sigmab(im) - + logical flag_shallow c physical parameters ! parameter(grav=grav,asolfac=0.958) ! parameter(elocp=hvap/cp,el2orc=hvap*hvap/(rv*cp)) @@ -1729,7 +1729,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & enddo enddo - if(progsigma)then + if(progsigma)then do k = 2, km1 do i = 1, im if (cnvflg(i)) then @@ -1776,8 +1776,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & c !> - Calculate the mean updraft velocity within the cloud (wc),cast in pressure coordinates. - if(progsigma)then - + if(progsigma)then do i = 1, im omegac(i) = 0. sumx(i) = 0. @@ -1861,17 +1860,19 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & c c store term needed for "termC" in prognostic area fraction closure - do k = 2, km1 - do i = 1, im - dp = 1000. * del(i,k) - if (cnvflg(i)) then - if(k > kbcon(i) .and. k < ktcon(i)) then - zdqca(i,k)=((qlks(i,k)-qlks(i,k-1)) + - & pwo(i,k)+dellal(i,k)) + if(progsigma)then + do k = 2, km1 + do i = 1, im + dp = 1000. * del(i,k) + if (cnvflg(i)) then + if(k > kbcon(i) .and. k < ktcon(i)) then + zdqca(i,k)=((qlks(i,k)-qlks(i,k-1)) + + & pwo(i,k)+dellal(i,k)) + endif endif - endif + enddo enddo - enddo + endif ccccc if(lat.==.latd.and.lon.==.lond.and.cnvflg(i)) then ccccc print *, ' aa1(i) before dwndrft =', aa1(i) @@ -2890,10 +2891,10 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & !> - From Bengtsson et al. (2022) Prognostic closure scheme, equation 8, compute updraft area fraction based on a moisture budget if(progsigma)then - call progsigma_calc(im,km,first_time_step,restart, + flag_shallow = .false. + call progsigma_calc(im,km,first_time_step,restart,flag_shallow, & del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt, & qgrs_dsave,q,kbcon1,ktcon,cnvflg,gdx, - & do_ca, ca_closure, ca_entr, ca_trigger, nthresh, ca_deep, & ca_micro,sigmain,sigmaout,sigmab,errmsg,errflg) endif diff --git a/physics/samfdeepcnv.meta b/physics/samfdeepcnv.meta index 3eb330551..71f9b87a5 100644 --- a/physics/samfdeepcnv.meta +++ b/physics/samfdeepcnv.meta @@ -70,15 +70,15 @@ type = logical intent = in [tmf] - standard_name = turbulence_moisture_flux - long_name = turbulence_moisture_flux + standard_name = turbulence_moisture_flux_for_coupling_to_convection + long_name = turbulence_moisture_flux_for_coupling_to_convection units = kg kg-1 s-1 dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys intent = in [qmicro] - standard_name = instantanious_moisture_tendency_due_to_microphysics + standard_name = instantaneous_moisture_tendency_due_to_microphysics long_name = moisture tendency due to microphysics units = kg kg-1 s-1 dimensions = (horizontal_loop_extent,vertical_layer_dimension) @@ -450,13 +450,13 @@ kind = kind_phys intent = inout [sigmain] - standard_name = updraft_area_fraction + standard_name = prognostic_updraft_area_fraction_in_convection long_name = convective updraft area fraction units = frac dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys - intent = inout + intent = in [sigmaout] standard_name = updraft_area_fraction_updated_by_physics long_name = convective updraft area fraction updated by physics diff --git a/physics/samfshalcnv.f b/physics/samfshalcnv.f index 24e01b040..6a682e9eb 100644 --- a/physics/samfshalcnv.f +++ b/physics/samfshalcnv.f @@ -57,11 +57,13 @@ end subroutine samfshalcnv_finalize !! @{ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & eps,epsm1,fv,grav,hvap,rd,rv, & - & t0c,delt,ntk,ntr,delp, & - & prslp,psp,phil,qtr,q1,t1,u1,v1,fscav, & + & t0c,delt,ntk,ntr,delp,first_time_step,restart, & + & tmf,qmicro,progsigma, & + & prslp,psp,phil,qtr,qgrs_dsave,q,q1,t1,u1,v1,fscav, & & rn,kbot,ktop,kcnv,islimsk,garea, & & dot,ncloud,hpbl,ud_mf,dt_mf,cnvw,cnvc, & - & clam,c0s,c1,evef,pgcon,asolfac,hwrf_samfshal,errmsg,errflg) + & clam,c0s,c1,evef,pgcon,asolfac,hwrf_samfshal, + & ca_micro,sigmain,sigmaout,errmsg,errflg) ! use machine , only : kind_phys use funcphys , only : fpvs @@ -74,7 +76,8 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & eps, epsm1, fv, grav, hvap, rd, rv, t0c real(kind=kind_phys), intent(in) :: delt real(kind=kind_phys), intent(in) :: psp(:), delp(:,:), & - & prslp(:,:), garea(:), hpbl(:), dot(:,:), phil(:,:) + & prslp(:,:), garea(:), hpbl(:), dot(:,:), phil(:,:), & + & qmicro(:,:),tmf(:,:),qgrs_dsave(:,:),q(:,:),sigmain(:,:) ! real(kind=kind_phys), dimension(:), intent(in) :: fscav integer, intent(inout) :: kcnv(:) @@ -83,12 +86,13 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & q1(:,:), t1(:,:), u1(:,:), v1(:,:) ! integer, intent(out) :: kbot(:), ktop(:) - real(kind=kind_phys), intent(out) :: rn(:), & - & cnvw(:,:), cnvc(:,:), ud_mf(:,:), dt_mf(:,:) + real(kind=kind_phys), intent(out) :: rn(:), ca_micro(:), & + & cnvw(:,:), cnvc(:,:), ud_mf(:,:), dt_mf(:,:), sigmaout(:,:) ! real(kind=kind_phys), intent(in) :: clam, c0s, c1, & & asolfac, evef, pgcon - logical, intent(in) :: hwrf_samfshal + logical, intent(in) :: hwrf_samfshal,first_time_step, & + & restart,progsigma character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg ! @@ -155,6 +159,13 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & bb1, bb2, wucb cc + +! parameters for prognostic sigma closure + real(kind=kind_phys) omega_u(im,km),zdqca(im,km),qlks(im,km), + & omegac(im),zeta(im,km),dbyo1(im,km), + & sigmab(im) + logical flag_shallow + c physical parameters ! parameter(g=grav,asolfac=0.89) ! parameter(g=grav) @@ -323,6 +334,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & ! vshear(i) = 0. gdx(i) = sqrt(garea(i)) xmb(i) = 0. + ca_micro(i) = 0. enddo endif !! @@ -498,6 +510,21 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & endif enddo enddo + + + do i = 1,im + omegac(i)=0. + enddo + + do k = 1, km + do i = 1, im + dbyo1(i,k)=0. + zdqca(i,k)=0. + qlks(i,k)=0. + omega_u(i,k)=0. + zeta(i,k)=1.0 + enddo + enddo ! ! initialize tracer variables ! @@ -1237,6 +1264,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & qcko(i,k)= qlk + qrch pwo(i,k) = etah * c0t(i,k) * dz * qlk cnvwt(i,k) = etah * qlk * grav / dp + qlks(i,k)=qlk endif ! ! compute buoyancy and drag for updraft velocity @@ -1304,6 +1332,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & if(k >= kbcon(i) .and. k < ktcon(i)) then dz1 = zo(i,k+1) - zo(i,k) aa1(i) = aa1(i) + buo(i,k) * dz1 + dbyo1(i,k) = hcko(i,k) - heso(i,k) endif endif enddo @@ -1402,6 +1431,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & qcko(i,k) = qlk + qrch pwo(i,k) = etah * c0t(i,k) * dz * qlk cnvwt(i,k) = etah * qlk * grav / dp + qlks(i,k)=qlk endif endif endif @@ -1444,6 +1474,20 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & enddo enddo ! + if(progsigma)then + do k = 2, km1 + do i = 1, im + if (cnvflg(i)) then + if(k > kbcon1(i) .and. k < ktcon(i)) then + rho = po(i,k)*100. / (rd * to(i,k)) + omega_u(i,k)=-1.0*sqrt(wu2(i,k))*rho*grav + omega_u(i,k)=MAX(omega_u(i,k),-80.) + endif + endif + enddo + enddo + endif + ! compute updraft velocity averaged over the whole cumulus ! !> - Calculate the mean updraft velocity within the cloud (wc). @@ -1475,6 +1519,50 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & endif enddo c +!> - Calculate the mean updraft velocity in pressure coordinates within the cloud (wc). + if(progsigma)then + do i = 1, im + omegac(i) = 0. + sumx(i) = 0. + enddo + do k = 2, km1 + do i = 1, im + if (cnvflg(i)) then + if(k > kbcon1(i) .and. k < ktcon(i)) then + dp = 1000. * del(i,k) + tem = 0.5 * (omega_u(i,k) + omega_u(i,k-1)) + omegac(i) = omegac(i) + tem * dp + sumx(i) = sumx(i) + dp + endif + endif + enddo + enddo + do i = 1, im + if(cnvflg(i)) then + if(sumx(i) == 0.) then + cnvflg(i)=.false. + else + omegac(i) = omegac(i) / sumx(i) + endif + val = -1.2 + if (omegac(i) > val) cnvflg(i)=.false. + endif + enddo +c +! > - Calculate the mean updraft velocity within the cloud (omega). + do k = 2, km1 + do i = 1, im + if (cnvflg(i)) then + if(k > kbcon1(i) .and. k < ktcon(i)) then + zeta(i,k)=eta(i,k)*(omegac(i)/omega_u(i,k)) + zeta(i,k)=MAX(0.,zeta(i,k)) + zeta(i,k)=MIN(1.,zeta(i,k)) + endif + endif + enddo + enddo + endif !if progsigma + c exchange ktcon with ktcon1 c do i = 1, im @@ -1505,11 +1593,25 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & if(dq > 0.) then qlko_ktcon(i) = dq qcko(i,k) = qrch + qlks(i,k) = qlko_ktcon(i) endif endif enddo endif c + + do k = 2, km1 + do i = 1, im + dp = 1000. * del(i,k) + if (cnvflg(i)) then + if(k > kbcon(i) .and. k < ktcon(i)) then + zdqca(i,k)=((qlks(i,k)-qlks(i,k-1)) + + & pwo(i,k)+dellal(i,k)) + endif + endif + enddo + enddo + c--- compute precipitation efficiency in terms of windshear c !! - Calculate the wind shear and precipitation efficiency according to equation 58 in Fritsch and Chappell (1980) \cite fritsch_and_chappell_1980 : @@ -1824,13 +1926,26 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & c compute cloud base mass flux as a function of the mean c updraft velcoity c +c Prognostic closure + if(progsigma)then + flag_shallow = .true. + call progsigma_calc(im,km,first_time_step,restart,flag_shallow, + & del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt, + & qgrs_dsave,q,kbcon1,ktcon,cnvflg,gdx, + & ca_micro,sigmain,sigmaout,sigmab,errmsg,errflg) + 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. !! As discussed in Han et al. (2017) \cite han_et_al_2017 , when dtconv is larger than tauadv, the convective mixing is not fully conducted before the cumulus cloud is advected out of the grid cell. In this case, therefore, the cloud base mass flux is further reduced in proportion to the ratio of tauadv to dtconv. do i= 1, im if(cnvflg(i)) then k = kbcon(i) rho = po(i,k)*100. / (rd*to(i,k)) - xmb(i) = advfac(i)*betaw*rho*wc(i) + if(progsigma)then + xmb(i) = sigmab(i)*((-1.0*omegac(i))/grav) + else + xmb(i) = advfac(i)*betaw*rho*wc(i) + endif endif enddo ! @@ -1850,8 +1965,12 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & do i = 1, im if(cnvflg(i)) then if (gdx(i) < dxcrt) then - scaldfunc(i) = (1.-sigmagfm(i)) * (1.-sigmagfm(i)) - scaldfunc(i) = max(min(scaldfunc(i), 1.0), 0.) + if(progsigma)then + scaldfunc(i)=(1.-sigmab(i))*(1.-sigmab(i)) + else + scaldfunc(i) = (1.-sigmagfm(i)) * (1.-sigmagfm(i)) + endif + scaldfunc(i) = max(min(scaldfunc(i), 1.0), 0.) else scaldfunc(i) = 1.0 endif diff --git a/physics/samfshalcnv.meta b/physics/samfshalcnv.meta index d768d4451..4383b8a67 100644 --- a/physics/samfshalcnv.meta +++ b/physics/samfshalcnv.meta @@ -1,7 +1,7 @@ [ccpp-table-properties] name = samfshalcnv type = scheme - dependencies = funcphys.f90,machine.F,samfaerosols.F + dependencies = funcphys.f90,machine.F,samfaerosols.F,progsigma_calc.f90 ######################################################################## [ccpp-arg-table] @@ -55,6 +55,36 @@ dimensions = () type = integer intent = in +[first_time_step] + standard_name = flag_for_first_timestep + long_name = flag for first time step for time integration loop (cold/warmstart) + units = flag + dimensions = () + type = logical + intent = in +[restart] + standard_name = flag_for_restart + long_name = flag for restart (warmstart) or coldstart + units = flag + dimensions = () + type = logical + intent = in +[tmf] + standard_name = turbulence_moisture_flux_for_coupling_to_convection + long_name = turbulence_moisture_flux_for_coupling_to_convection + units = kg kg-1 s-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in +[qmicro] + standard_name = instantaneous_moisture_tendency_due_to_microphysics + long_name = moisture tendency due to microphysics + units = kg kg-1 s-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in [itc] standard_name = index_of_first_chemical_tracer_for_convection long_name = index of first chemical tracer transported/scavenged by convection @@ -219,6 +249,22 @@ type = real kind = kind_phys intent = inout +[qgrs_dsave] + standard_name = tracer_concentration_dsave + long_name = model layer mean tracer concentration dsave + units = kg kg-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in +[q] + standard_name = specific_humidity + long_name = water vapor specific humidity + units = kg kg-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in [q1] standard_name = specific_humidity_of_new_state long_name = updated vapor specific humidity @@ -413,6 +459,37 @@ dimensions = () type = logical intent = in +[progsigma] + standard_name = flag_for_prognostic_sigma + long_name = flag for prognostic sigma + units = flag + dimensions = () + type = logical + intent = in +[ca_micro] + standard_name = output_prognostic_sigma_two + long_name = output of prognostic area fraction two + units = frac + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out +[sigmain] + standard_name = prognostic_updraft_area_fraction_in_convection + long_name = convective updraft area fraction + units = frac + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in +[sigmaout] + standard_name = updraft_area_fraction_updated_by_physics + long_name = convective updraft area fraction updated by physics + units = frac + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = out [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP diff --git a/physics/satmedmfvdifq.F b/physics/satmedmfvdifq.F index b0367b4d4..0fce7dd9a 100644 --- a/physics/satmedmfvdifq.F +++ b/physics/satmedmfvdifq.F @@ -121,9 +121,9 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & real(kind=kind_phys), intent(out) :: & & dusfc(:), dvsfc(:), & & dtsfc(:), dqsfc(:), & - & hpbl(:), tmf(:,:) + & hpbl(:) real(kind=kind_phys), intent(out) :: & - & dkt(:,:), dku(:,:) + & dkt(:,:), dku(:,:), tmf(:,:) ! logical, intent(in) :: dspheat character(len=*), intent(out) :: errmsg @@ -299,6 +299,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & xmfd(i,k) = 0. buou(i,k) = 0. buod(i,k) = 0. + tmf(i,k) = 0. ckz(i,k) = ck1 chz(i,k) = ch1 rlmnz(i,k) = rlmn0 diff --git a/physics/satmedmfvdifq.meta b/physics/satmedmfvdifq.meta index 211bbaec6..9b803e4a5 100644 --- a/physics/satmedmfvdifq.meta +++ b/physics/satmedmfvdifq.meta @@ -202,10 +202,10 @@ kind = kind_phys intent = inout [tmf] - standard_name = turbulence_moisture_flux - long_name = turbulence_moisture_flux + standard_name = turbulence_moisture_flux_for_coupling_to_convection + long_name = turbulence_moisture_flux_for_coupling_to_convection units = kg kg-1 s-1 - dimensions = (horizontal_dimension,vertical_layer_dimension) + dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys intent = out From b530db11801ca67a10e2757aef04ed9492a06e7e Mon Sep 17 00:00:00 2001 From: Lisa Bengtsson Date: Wed, 20 Apr 2022 01:29:13 +0000 Subject: [PATCH 05/20] cleaning some diagnostics --- physics/progsigma_calc.f90 | 5 +---- physics/samfdeepcnv.f | 8 +++----- physics/samfdeepcnv.meta | 8 -------- physics/samfshalcnv.f | 7 +++---- physics/samfshalcnv.meta | 8 -------- 5 files changed, 7 insertions(+), 29 deletions(-) diff --git a/physics/progsigma_calc.f90 b/physics/progsigma_calc.f90 index ca05f6778..7673602b6 100644 --- a/physics/progsigma_calc.f90 +++ b/physics/progsigma_calc.f90 @@ -16,7 +16,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & flag_shallow,del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap, & delt,qgrs_dsave,q,kbcon1,ktcon,cnvflg,gdx, & - ca_micro,sigmain,sigmaout,sigmab,errmsg,errflg) + sigmain,sigmaout,sigmab,errmsg,errflg) ! ! use machine, only : kind_phys @@ -31,7 +31,6 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & qmicro(im,km),tmf(im,km),dbyo1(im,km),zdqca(im,km), & omega_u(im,km),zeta(im,km),gdx(im) logical, intent(in) :: flag_init,flag_restart,cnvflg(im),flag_shallow - real(kind=kind_phys), intent(out):: ca_micro(im) real(kind=kind_phys), intent(in) :: sigmain(im,km) ! intent out @@ -79,7 +78,6 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & termD(i)=0. fdqa(i)=0. mcons(i)=0. - ca_micro(i)=0. enddo !Initial computations, place maximum sigmain in sigmab @@ -207,7 +205,6 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & sigmab(i)=MIN(sigmab(i),sigmamax(i)) sigmab(i)=MAX(sigmab(i),0.01) endif - ca_micro(i)=sigmab(i) endif!cnvflg enddo diff --git a/physics/samfdeepcnv.f b/physics/samfdeepcnv.f index 35aea0eb1..45cbf70e1 100644 --- a/physics/samfdeepcnv.f +++ b/physics/samfdeepcnv.f @@ -86,7 +86,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & & CNV_DQLDT,CLCN,CNV_FICE,CNV_NDROP,CNV_NICE,mp_phys,mp_phys_mg,& & clam,c0s,c1,betal,betas,evef,pgcon,asolfac, & & do_ca, ca_closure, ca_entr, ca_trigger, nthresh, ca_deep, & - & rainevap, sigmain, sigmaout, ca_micro, & + & rainevap, sigmain, sigmaout, & & errmsg,errflg) ! use machine , only : kind_phys @@ -108,7 +108,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & real(kind=kind_phys), intent(in) :: ca_deep(:) real(kind=kind_phys), intent(in) :: sigmain(:,:),qmicro(:,:), & & tmf(:,:),q(:,:), qgrs_dsave(:,:) - real(kind=kind_phys), intent(out) :: rainevap(:),ca_micro(:) + real(kind=kind_phys), intent(out) :: rainevap(:) real(kind=kind_phys), intent(out) :: sigmaout(:,:) logical, intent(in) :: do_ca,ca_closure,ca_entr,ca_trigger @@ -919,8 +919,6 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & enddo if(totflg) return !! -! -!Lisa: at this point only trigger criteria is set ! turbulent entrainment rate assumed to be proportional ! to subcloud mean TKE @@ -2895,7 +2893,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & call progsigma_calc(im,km,first_time_step,restart,flag_shallow, & del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt, & qgrs_dsave,q,kbcon1,ktcon,cnvflg,gdx, - & ca_micro,sigmain,sigmaout,sigmab,errmsg,errflg) + & sigmain,sigmaout,sigmab,errmsg,errflg) 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 for the grid sizes where the quasi-equilibrium assumption of Arakawa-Schubert is not valid any longer. diff --git a/physics/samfdeepcnv.meta b/physics/samfdeepcnv.meta index 71f9b87a5..71c78036d 100644 --- a/physics/samfdeepcnv.meta +++ b/physics/samfdeepcnv.meta @@ -326,14 +326,6 @@ dimensions = () type = logical intent = in -[ca_micro] - standard_name = output_prognostic_sigma_two - long_name = output of prognostic area fraction two - units = frac - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = out [cldwrk] standard_name = cloud_work_function long_name = cloud work function diff --git a/physics/samfshalcnv.f b/physics/samfshalcnv.f index 6a682e9eb..343691279 100644 --- a/physics/samfshalcnv.f +++ b/physics/samfshalcnv.f @@ -63,7 +63,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & rn,kbot,ktop,kcnv,islimsk,garea, & & dot,ncloud,hpbl,ud_mf,dt_mf,cnvw,cnvc, & & clam,c0s,c1,evef,pgcon,asolfac,hwrf_samfshal, - & ca_micro,sigmain,sigmaout,errmsg,errflg) + & sigmain,sigmaout,errmsg,errflg) ! use machine , only : kind_phys use funcphys , only : fpvs @@ -86,7 +86,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & q1(:,:), t1(:,:), u1(:,:), v1(:,:) ! integer, intent(out) :: kbot(:), ktop(:) - real(kind=kind_phys), intent(out) :: rn(:), ca_micro(:), & + real(kind=kind_phys), intent(out) :: rn(:), & & cnvw(:,:), cnvc(:,:), ud_mf(:,:), dt_mf(:,:), sigmaout(:,:) ! real(kind=kind_phys), intent(in) :: clam, c0s, c1, & @@ -334,7 +334,6 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & ! vshear(i) = 0. gdx(i) = sqrt(garea(i)) xmb(i) = 0. - ca_micro(i) = 0. enddo endif !! @@ -1932,7 +1931,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & call progsigma_calc(im,km,first_time_step,restart,flag_shallow, & del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt, & qgrs_dsave,q,kbcon1,ktcon,cnvflg,gdx, - & ca_micro,sigmain,sigmaout,sigmab,errmsg,errflg) + & sigmain,sigmaout,sigmab,errmsg,errflg) 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. diff --git a/physics/samfshalcnv.meta b/physics/samfshalcnv.meta index 4383b8a67..895460ffd 100644 --- a/physics/samfshalcnv.meta +++ b/physics/samfshalcnv.meta @@ -466,14 +466,6 @@ dimensions = () type = logical intent = in -[ca_micro] - standard_name = output_prognostic_sigma_two - long_name = output of prognostic area fraction two - units = frac - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = out [sigmain] standard_name = prognostic_updraft_area_fraction_in_convection long_name = convective updraft area fraction From fc7e7a0e226663e0fec80729097a7dc69ed5551d Mon Sep 17 00:00:00 2001 From: Lisa Bengtsson Date: Fri, 22 Apr 2022 04:00:31 +0000 Subject: [PATCH 06/20] addressing some review comments --- physics/GFS_MP_generic.F90 | 16 +++--- physics/GFS_MP_generic.meta | 12 ++--- physics/progsigma_calc.f90 | 103 ++++++++++++++++++------------------ physics/samfdeepcnv.f | 10 ++-- physics/samfdeepcnv.meta | 14 ++--- physics/samfshalcnv.f | 12 +++-- physics/samfshalcnv.meta | 14 ++--- physics/satmedmfvdifq.meta | 4 +- 8 files changed, 95 insertions(+), 90 deletions(-) diff --git a/physics/GFS_MP_generic.F90 b/physics/GFS_MP_generic.F90 index dbf2d15fa..f8dd8ab6c 100644 --- a/physics/GFS_MP_generic.F90 +++ b/physics/GFS_MP_generic.F90 @@ -95,7 +95,7 @@ subroutine GFS_MP_generic_post_run( drain_cpl, dsnow_cpl, lsm, lsm_ruc, lsm_noahmp, raincprv, rainncprv, iceprv, snowprv, & graupelprv, draincprv, drainncprv, diceprv, dsnowprv, dgraupelprv, dtp, dfi_radar_max_intervals, & dtend, dtidx, index_of_temperature, index_of_process_mp,ldiag3d, qdiag3d,dqdt_qmicro, lssav, num_dfi_radar, & - fh_dfi_radar,index_of_process_dfi_radar, ix_dfi_radar, dfi_radar_tten, radar_tten_limits, fhour, qgrs_dsave, & + fh_dfi_radar,index_of_process_dfi_radar, ix_dfi_radar, dfi_radar_tten, radar_tten_limits, fhour, prevsq, & errmsg, errflg) ! use machine, only: kind_phys @@ -150,7 +150,7 @@ subroutine GFS_MP_generic_post_run( real(kind=kind_phys), dimension(:), intent(inout) :: dsnowprv real(kind=kind_phys), dimension(:), intent(inout) :: dgraupelprv real(kind=kind_phys), dimension(:,:), intent(out) :: dqdt_qmicro - real(kind=kind_phys), dimension(:,:), intent(out) :: qgrs_dsave + real(kind=kind_phys), dimension(:,:), intent(out) :: prevsq real(kind=kind_phys), intent(in) :: dtp ! CCPP error handling @@ -466,11 +466,13 @@ subroutine GFS_MP_generic_post_run( pwat(i) = pwat(i) * onebg enddo - do k = 1, levs - do i=1, im - qgrs_dsave(i,k) = gq0(i,k,1) - enddo - enddo + if(progsigma)then + do k = 1, levs + do i=1, im + prevsq(i,k) = gq0(i,k,1) + enddo + enddo + endif end subroutine GFS_MP_generic_post_run !> @} diff --git a/physics/GFS_MP_generic.meta b/physics/GFS_MP_generic.meta index 763cad85a..eb0f17fa8 100644 --- a/physics/GFS_MP_generic.meta +++ b/physics/GFS_MP_generic.meta @@ -249,8 +249,8 @@ type = logical intent = in [progsigma] - standard_name = flag_for_prognostic_sigma - long_name = flag for prognostic sigma + standard_name = do_prognostic_updraft_area_fraction + long_name = flag for prognostic sigma in cumulus scheme units = flag dimensions = () type = logical @@ -852,16 +852,16 @@ type = logical intent = in [dqdt_qmicro] - standard_name = instantaneous_moisture_tendency_due_to_microphysics + standard_name = instantaneous_tendency_of_specific_humidity_due_to_microphysics long_name = moisture tendency due to microphysics units = kg kg-1 s-1 dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys intent = out -[qgrs_dsave] - standard_name = tracer_concentration_dsave - long_name = model layer mean tracer concentration dsave +[prevsq] + standard_name = specific_humidity_on_previous_timestep + long_name = specific_humidity_on_previous_timestep units = kg kg-1 dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real diff --git a/physics/progsigma_calc.f90 b/physics/progsigma_calc.f90 index 7673602b6..58a6fc0ef 100644 --- a/physics/progsigma_calc.f90 +++ b/physics/progsigma_calc.f90 @@ -1,4 +1,4 @@ -!>\file progsigma +!>\file progsigma_calc.f90 !! This file contains the subroutine that calculates the prognostic !! updraft area fraction that is used for closure computations in !! saSAS deep and shallow convection, based on a moisture budget @@ -15,7 +15,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & flag_shallow,del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap, & - delt,qgrs_dsave,q,kbcon1,ktcon,cnvflg,gdx, & + delt,prevsq,q,kbcon1,ktcon,cnvflg,gdx, & sigmain,sigmaout,sigmab,errmsg,errflg) ! ! @@ -27,7 +27,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & ! intent in integer, intent(in) :: im,km,kbcon1(im),ktcon(im) real, intent(in) :: hvap,delt - real, intent(in) :: qgrs_dsave(im,km), q(im,km),del(im,km), & + real, intent(in) :: prevsq(im,km), q(im,km),del(im,km), & qmicro(im,km),tmf(im,km),dbyo1(im,km),zdqca(im,km), & omega_u(im,km),zeta(im,km),gdx(im) logical, intent(in) :: flag_init,flag_restart,cnvflg(im),flag_shallow @@ -43,33 +43,32 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & integer :: i,k,km1 real(kind=kind_phys) :: termA(im),termB(im),termC(im),termD(im), & mcons(im),fdqa(im),form(im,km), & - qadv(im,km),sigmamax(im) + qadv(im,km),sigmamax(im),dp(im),inbu(im,km) real(kind=kind_phys) :: gcvalmx,epsilon,ZZ,cvg,mcon,buy2, & - fdqb,dtdyn,dxlim,rmulacvg,dp,tem, & - alpha,DEN,betascu - integer :: inbu(im,km) + fdqb,dtdyn,dxlim,rmulacvg,tem, & + alpha,DEN,betascu,dp1 !Parameters - gcvalmx = 0.1 - rmulacvg=10. - epsilon=1.E-11 - km1=km-1 - alpha=7000. - betascu = 3.0 + gcvalmx = 0.1 + rmulacvg=10. + epsilon=1.E-11 + km1=km-1 + alpha=7000. + betascu = 3.0 !Initialization 2D - do k = 1,km - do i = 1,im - sigmaout(i,k)=0. - inbu(i,k)=0 - form(i,k)=0. - enddo - enddo + do k = 1,km + do i = 1,im + sigmaout(i,k)=0. + inbu(i,k)=0. + form(i,k)=0. + enddo + enddo !Initialization 1D - do i=1,im + do i=1,im sigmab(i)=0. sigmamax(i)=0.95 termA(i)=0. @@ -80,23 +79,32 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & mcons(i)=0. enddo - !Initial computations, place maximum sigmain in sigmab + do k = 2,km1 + do i = 1,im + if(cnvflg(i))then + dp(i) = 1000. * del(i,k) + endif + enddo + enddo - do k=2,km - do i=1,im - if(flag_init .and. .not. flag_restart)then + !Initial computations, place maximum sigmain in sigmab + if(flag_init .and. .not. flag_restart)then + do i=1,im if(cnvflg(i))then sigmab(i)=0.03 endif - else + enddo + else + do i=1,im if(cnvflg(i))then - if(sigmain(i,k)>sigmab(i))then - sigmab(i)=sigmain(i,k) - endif + do k=2,km + if(sigmain(i,k)>sigmab(i))then + sigmab(i)=sigmain(i,k) + endif + enddo endif - endif - enddo - enddo + enddo + endif do i=1,im if(sigmab(i) < 1.E-5)then !after advection @@ -116,7 +124,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & if(flag_init .and. .not.flag_restart)then qadv(i,k)=0. else - qadv(i,k)=(q(i,k) - qgrs_dsave(i,k))/delt + qadv(i,k)=(q(i,k) - prevsq(i,k))/delt endif enddo enddo @@ -125,22 +133,21 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & !buoyant layers with positive moisture convergence (accumulated from the surface). !Lowest level: do i = 1,im - dp = 1000. * del(i,1) - mcons(i)=(hvap*(qadv(i,1)+tmf(i,1)+qmicro(i,1))*dp) + dp1 = 1000. * del(i,1) + mcons(i)=(hvap*(qadv(i,1)+tmf(i,1)+qmicro(i,1))*dp1) enddo !Levels above: do k = 2,km1 do i = 1,im - dp = 1000. * del(i,k) if(cnvflg(i))then - mcon = (hvap*(qadv(i,k)+tmf(i,k)+qmicro(i,k))*dp) + mcon = (hvap*(qadv(i,k)+tmf(i,k)+qmicro(i,k))*dp(i)) buy2 = termD(i)+mcon+mcons(i) ! Do the integral over buoyant layers with positive mcon acc from surface if(k > kbcon1(i) .and. k < ktcon(i) .and. buy2 > 0.)then - inbu(i,k)=1 + inbu(i,k)=1. endif inbu(i,k-1)=MAX(inbu(i,k-1),inbu(i,k)) - termD(i) = termD(i) + float(inbu(i,k-1))*mcons(i) + termD(i) = termD(i) + inbu(i,k-1)*mcons(i) mcons(i)=mcon endif enddo @@ -149,9 +156,8 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & !termA do k = 2,km1 do i = 1,im - dp = 1000. * del(i,k) if(cnvflg(i))then - tem=(sigmab(i)*zeta(i,k)*float(inbu(i,k))*dbyo1(i,k))*dp + tem=(sigmab(i)*zeta(i,k)*inbu(i,k)*dbyo1(i,k))*dp(i) termA(i)=termA(i)+tem endif enddo @@ -160,9 +166,8 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & !termB do k = 2,km1 do i = 1,im - dp = 1000. * del(i,k) if(cnvflg(i))then - tem=(dbyo1(i,k)*float(inbu(i,k)))*dp + tem=(dbyo1(i,k)*inbu(i,k))*dp(i) termB(i)=termB(i)+tem endif enddo @@ -172,10 +177,9 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & do k = 2,km1 do i = 1,im if(cnvflg(i))then - dp = 1000. * del(i,k) - form(i,k)=-1.0*float(inbu(i,k))*(omega_u(i,k)*delt) + form(i,k)=-1.0*inbu(i,k)*(omega_u(i,k)*delt) fdqb=0.5*((form(i,k)*zdqca(i,k))) - termC(i)=termC(i)+(float(inbu(i,k))* & + termC(i)=termC(i)+inbu(i,k)* & (fdqb+fdqa(i))*hvap*zeta(i,k)) fdqa(i)=fdqb endif @@ -185,22 +189,17 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & !sigmab do i = 1,im if(cnvflg(i))then - DEN=MIN(termC(i)+termB(i),1.E8) cvg=termD(i)*delt ZZ=MAX(0.0,SIGN(1.0,termA(i))) & *MAX(0.0,SIGN(1.0,termB(i))) & - *MAX(0.0,SIGN(1.0,termC(i)-epsilon)) - - + *MAX(0.0,SIGN(1.0,termC(i)-epsilon)) cvg=MAX(0.0,cvg) - if(flag_init .and. .not. flag_restart)then sigmab(i)=0.03 else sigmab(i)=(ZZ*(termA(i)+cvg))/(DEN+(1.0-ZZ)) endif - if(sigmab(i)>0.)then sigmab(i)=MIN(sigmab(i),sigmamax(i)) sigmab(i)=MAX(sigmab(i),0.01) diff --git a/physics/samfdeepcnv.f b/physics/samfdeepcnv.f index 45cbf70e1..02b2dcb83 100644 --- a/physics/samfdeepcnv.f +++ b/physics/samfdeepcnv.f @@ -79,7 +79,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & & tmf,qmicro,itc,ntc,cliq,cp,cvap, & & eps,epsm1,fv,grav,hvap,rd,rv, & & t0c,delt,ntk,ntr,delp, & - & prslp,psp,phil,qtr,qgrs_dsave,q,q1,t1,u1,v1,fscav, & + & prslp,psp,phil,qtr,prevsq,q,q1,t1,u1,v1,fscav, & & hwrf_samfdeep,progsigma,wclosureflg,cldwrk,rn,kbot,ktop,kcnv, & & islimsk,garea,dot,ncloud,hpbl,ud_mf,dd_mf,dt_mf,cnvw,cnvc, & & QLCN, QICN, w_upi, cf_upi, CNV_MFD, & @@ -107,7 +107,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & real(kind=kind_phys), intent(in) :: nthresh real(kind=kind_phys), intent(in) :: ca_deep(:) real(kind=kind_phys), intent(in) :: sigmain(:,:),qmicro(:,:), & - & tmf(:,:),q(:,:), qgrs_dsave(:,:) + & tmf(:,:),q(:,:), prevsq(:,:) real(kind=kind_phys), intent(out) :: rainevap(:) real(kind=kind_phys), intent(out) :: sigmaout(:,:) logical, intent(in) :: do_ca,ca_closure,ca_entr,ca_trigger @@ -217,6 +217,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & real(kind=kind_phys) omega_u(im,km),zdqca(im,km),qlks(im,km), & omegac(im),zeta(im,km),dbyo1(im,km),sigmab(im) logical flag_shallow + real(kind=kind_phys) gravinv c physical parameters ! parameter(grav=grav,asolfac=0.958) ! parameter(elocp=hvap/cp,el2orc=hvap*hvap/(rv*cp)) @@ -309,6 +310,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & errmsg = '' errflg = 0 + gravinv = 1./grav elocp = hvap/cp el2orc = hvap*hvap/(rv*cp) @@ -2892,7 +2894,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & flag_shallow = .false. call progsigma_calc(im,km,first_time_step,restart,flag_shallow, & del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt, - & qgrs_dsave,q,kbcon1,ktcon,cnvflg,gdx, + & prevsq,q,kbcon1,ktcon,cnvflg,gdx, & sigmain,sigmaout,sigmab,errmsg,errflg) endif @@ -2903,7 +2905,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & k = kbcon(i) rho = po(i,k)*100. / (rd*to(i,k)) if(progsigma)then - xmb(i) = sigmab(i)*((-1.0*omegac(i))/grav) + xmb(i) = advfac(i)*sigmab(i)*((-1.0*omegac(i))*gravinv) else xmb(i) = advfac(i)*betaw*rho*wc(i) endif diff --git a/physics/samfdeepcnv.meta b/physics/samfdeepcnv.meta index 71c78036d..e956d24ed 100644 --- a/physics/samfdeepcnv.meta +++ b/physics/samfdeepcnv.meta @@ -70,8 +70,8 @@ type = logical intent = in [tmf] - standard_name = turbulence_moisture_flux_for_coupling_to_convection - long_name = turbulence_moisture_flux_for_coupling_to_convection + standard_name = instantaneous_tendency_of_specific_humidity_due_to_PBL + long_name = instantaneous_tendency_of_specific_humidity_due_to_PBL units = kg kg-1 s-1 dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real @@ -249,9 +249,9 @@ type = real kind = kind_phys intent = inout -[qgrs_dsave] - standard_name = tracer_concentration_dsave - long_name = model layer mean tracer concentration dsave +[prevsq] + standard_name = specific_humidity_on_previous_timestep + long_name = specific_humidity_on_previous_timestep units = kg kg-1 dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real @@ -320,8 +320,8 @@ type = logical intent = in [progsigma] - standard_name = flag_for_prognostic_sigma - long_name = flag for prognostic sigma + standard_name = do_prognostic_updraft_area_fraction + long_name = flag for prognostic sigma in cumuls scheme units = flag dimensions = () type = logical diff --git a/physics/samfshalcnv.f b/physics/samfshalcnv.f index 343691279..c3bb842b3 100644 --- a/physics/samfshalcnv.f +++ b/physics/samfshalcnv.f @@ -59,7 +59,7 @@ subroutine samfshalcnv_run(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, & - & prslp,psp,phil,qtr,qgrs_dsave,q,q1,t1,u1,v1,fscav, & + & prslp,psp,phil,qtr,prevsq,q,q1,t1,u1,v1,fscav, & & rn,kbot,ktop,kcnv,islimsk,garea, & & dot,ncloud,hpbl,ud_mf,dt_mf,cnvw,cnvc, & & clam,c0s,c1,evef,pgcon,asolfac,hwrf_samfshal, @@ -77,7 +77,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & real(kind=kind_phys), intent(in) :: delt real(kind=kind_phys), intent(in) :: psp(:), delp(:,:), & & prslp(:,:), garea(:), hpbl(:), dot(:,:), phil(:,:), & - & qmicro(:,:),tmf(:,:),qgrs_dsave(:,:),q(:,:),sigmain(:,:) + & qmicro(:,:),tmf(:,:),prevsq(:,:),q(:,:),sigmain(:,:) ! real(kind=kind_phys), dimension(:), intent(in) :: fscav integer, intent(inout) :: kcnv(:) @@ -165,6 +165,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & omegac(im),zeta(im,km),dbyo1(im,km), & sigmab(im) logical flag_shallow + real(kind=kind_phys) gravinv c physical parameters ! parameter(g=grav,asolfac=0.89) @@ -249,6 +250,8 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & errmsg = '' errflg = 0 + gravinv = 1./grav + elocp = hvap/cp el2orc = hvap*hvap/(rv*cp) @@ -1601,7 +1604,6 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & do k = 2, km1 do i = 1, im - dp = 1000. * del(i,k) if (cnvflg(i)) then if(k > kbcon(i) .and. k < ktcon(i)) then zdqca(i,k)=((qlks(i,k)-qlks(i,k-1)) + @@ -1930,7 +1932,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & flag_shallow = .true. call progsigma_calc(im,km,first_time_step,restart,flag_shallow, & del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt, - & qgrs_dsave,q,kbcon1,ktcon,cnvflg,gdx, + & prevsq,q,kbcon1,ktcon,cnvflg,gdx, & sigmain,sigmaout,sigmab,errmsg,errflg) endif @@ -1941,7 +1943,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & k = kbcon(i) rho = po(i,k)*100. / (rd*to(i,k)) if(progsigma)then - xmb(i) = sigmab(i)*((-1.0*omegac(i))/grav) + xmb(i) = advfac(i)*sigmab(i)*((-1.0*omegac(i))*gravinv) else xmb(i) = advfac(i)*betaw*rho*wc(i) endif diff --git a/physics/samfshalcnv.meta b/physics/samfshalcnv.meta index 895460ffd..a4cca64b8 100644 --- a/physics/samfshalcnv.meta +++ b/physics/samfshalcnv.meta @@ -70,8 +70,8 @@ type = logical intent = in [tmf] - standard_name = turbulence_moisture_flux_for_coupling_to_convection - long_name = turbulence_moisture_flux_for_coupling_to_convection + standard_name = instantaneous_tendency_of_specific_humidity_due_to_PBL + long_name = instantaneous_tendency_of_specific_humidity_due_to_PBL units = kg kg-1 s-1 dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real @@ -249,9 +249,9 @@ type = real kind = kind_phys intent = inout -[qgrs_dsave] - standard_name = tracer_concentration_dsave - long_name = model layer mean tracer concentration dsave +[prevsq] + standard_name = specific_humidity_on_previous_timestep + long_name = specific_humidity_on_previous_timestep units = kg kg-1 dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real @@ -460,8 +460,8 @@ type = logical intent = in [progsigma] - standard_name = flag_for_prognostic_sigma - long_name = flag for prognostic sigma + standard_name = do_prognostic_updraft_area_fraction + long_name = flag for prognostic sigma in cumulus scheme units = flag dimensions = () type = logical diff --git a/physics/satmedmfvdifq.meta b/physics/satmedmfvdifq.meta index 9b803e4a5..fa30cd9f7 100644 --- a/physics/satmedmfvdifq.meta +++ b/physics/satmedmfvdifq.meta @@ -202,8 +202,8 @@ kind = kind_phys intent = inout [tmf] - standard_name = turbulence_moisture_flux_for_coupling_to_convection - long_name = turbulence_moisture_flux_for_coupling_to_convection + standard_name = instantaneous_tendency_of_specific_humidity_due_to_PBL + long_name = instantaneous_tendency_of_specific_humidity_due_to_PBL units = kg kg-1 s-1 dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real From 0200e2d05757af34a070e52fe8558ecbaa73a42a Mon Sep 17 00:00:00 2001 From: Lisa Bengtsson Date: Wed, 27 Apr 2022 18:48:38 +0000 Subject: [PATCH 07/20] addressing some review comments --- physics/progsigma_calc.f90 | 2 +- physics/samfdeepcnv.f | 2 +- physics/samfdeepcnv.meta | 2 +- physics/samfshalcnv.f | 2 +- physics/samfshalcnv.meta | 2 +- physics/satmedmfvdifq.F | 25 ++++++++++++++++++++----- physics/satmedmfvdifq.meta | 7 +++++++ 7 files changed, 32 insertions(+), 10 deletions(-) diff --git a/physics/progsigma_calc.f90 b/physics/progsigma_calc.f90 index 58a6fc0ef..55f6b5e3a 100644 --- a/physics/progsigma_calc.f90 +++ b/physics/progsigma_calc.f90 @@ -180,7 +180,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & form(i,k)=-1.0*inbu(i,k)*(omega_u(i,k)*delt) fdqb=0.5*((form(i,k)*zdqca(i,k))) termC(i)=termC(i)+inbu(i,k)* & - (fdqb+fdqa(i))*hvap*zeta(i,k)) + (fdqb+fdqa(i))*hvap*zeta(i,k) fdqa(i)=fdqb endif enddo diff --git a/physics/samfdeepcnv.f b/physics/samfdeepcnv.f index f8a60a5e3..5fd54a2ec 100644 --- a/physics/samfdeepcnv.f +++ b/physics/samfdeepcnv.f @@ -79,7 +79,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & & tmf,qmicro,itc,ntc,cliq,cp,cvap, & & eps,epsm1,fv,grav,hvap,rd,rv, & & t0c,delt,ntk,ntr,delp, & - & prslp,psp,phil,qtr,prevsq,q,q1,t1,u1,v1,fscav, & + & prslp,psp,phil,qtr,prevsq,q,q1,t1,u1,v1,fscav, & & hwrf_samfdeep,progsigma,wclosureflg,cldwrk,rn,kbot,ktop,kcnv, & & islimsk,garea,dot,ncloud,hpbl,ud_mf,dd_mf,dt_mf,cnvw,cnvc, & & QLCN, QICN, w_upi, cf_upi, CNV_MFD, & diff --git a/physics/samfdeepcnv.meta b/physics/samfdeepcnv.meta index e956d24ed..2b2942812 100644 --- a/physics/samfdeepcnv.meta +++ b/physics/samfdeepcnv.meta @@ -78,7 +78,7 @@ kind = kind_phys intent = in [qmicro] - standard_name = instantaneous_moisture_tendency_due_to_microphysics + standard_name = instantaneous_tendency_of_specific_humidity_due_to_microphysics long_name = moisture tendency due to microphysics units = kg kg-1 s-1 dimensions = (horizontal_loop_extent,vertical_layer_dimension) diff --git a/physics/samfshalcnv.f b/physics/samfshalcnv.f index 56571457a..325566877 100644 --- a/physics/samfshalcnv.f +++ b/physics/samfshalcnv.f @@ -59,7 +59,7 @@ subroutine samfshalcnv_run(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, & - & prslp,psp,phil,qtr,prevsq,q,q1,t1,u1,v1,fscav, & + & prslp,psp,phil,qtr,prevsq,q,q1,t1,u1,v1,fscav, & & rn,kbot,ktop,kcnv,islimsk,garea, & & dot,ncloud,hpbl,ud_mf,dt_mf,cnvw,cnvc, & & clam,c0s,c1,evef,pgcon,asolfac,hwrf_samfshal, diff --git a/physics/samfshalcnv.meta b/physics/samfshalcnv.meta index a4cca64b8..8c9735c32 100644 --- a/physics/samfshalcnv.meta +++ b/physics/samfshalcnv.meta @@ -78,7 +78,7 @@ kind = kind_phys intent = in [qmicro] - standard_name = instantaneous_moisture_tendency_due_to_microphysics + standard_name = instantaneous_tendency_of_specific_humidity_due_to_microphysics long_name = moisture tendency due to microphysics units = kg kg-1 s-1 dimensions = (horizontal_loop_extent,vertical_layer_dimension) diff --git a/physics/satmedmfvdifq.F b/physics/satmedmfvdifq.F index 0fce7dd9a..c7a6fadc9 100644 --- a/physics/satmedmfvdifq.F +++ b/physics/satmedmfvdifq.F @@ -67,8 +67,8 @@ end subroutine satmedmfvdifq_finalize !! (mfscuq.f). !! \section detail_satmedmfvidfq GFS satmedmfvdifq Detailed Algorithm !! @{ - subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & - & grav,rd,cp,rv,hvap,hfus,fv,eps,epsm1, & + subroutine satmedmfvdifq_run(im,km,progsigma,ntrac,ntcw,ntrw, & + & ntiw,ntke,grav,rd,cp,rv,hvap,hfus,fv,eps,epsm1, & & dv,du,tdt,rtg,tmf,u1,v1,t1,q1,swh,hlw,xmu,garea,zvfun, & & psk,rbsoil,zorl,u10m,v10m,fm,fh, & & tsea,heat,evap,stress,spd1,kpbl, & @@ -91,7 +91,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & integer, intent(in) :: sfc_rlm integer, intent(in) :: kinver(:) integer, intent(out) :: kpbl(:) - logical, intent(in) :: gen_tend,ldiag3d + logical, intent(in) :: gen_tend,ldiag3d,progsigma ! real(kind=kind_phys), intent(in) :: grav,rd,cp,rv,hvap,hfus,fv, & & eps,epsm1 @@ -299,7 +299,6 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & xmfd(i,k) = 0. buou(i,k) = 0. buod(i,k) = 0. - tmf(i,k) = 0. ckz(i,k) = ck1 chz(i,k) = ch1 rlmnz(i,k) = rlmn0 @@ -313,6 +312,14 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & zm(i,k) = zi(i,k+1) enddo enddo +!> - Initialize variables needed for prognostic cumulus closure + if(progsigma)then + do k=1,km + do i=1,im + tmf(i,k) = 0. + enddo + enddo + endif !> - Compute horizontal grid size (\p gdx) do i=1,im gdx(i) = sqrt(garea(i)) @@ -2115,11 +2122,19 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw,ntiw,ntke, & qtend = (f2(i,k)-q1(i,k,1))*rdt tdt(i,k) = tdt(i,k)+ttend rtg(i,k,1) = rtg(i,k,1)+qtend - tmf(i,k) = qtend ! dtsfc(i) = dtsfc(i)+cont*del(i,k)*ttend ! dqsfc(i) = dqsfc(i)+conq*del(i,k)*qtend enddo enddo + + if(progsigma)then + do k = 1,km + do i = 1,im + tmf(i,k)=(f2(i,k)-q1(i,k,1))*rdt + enddo + enddo + endif + do i = 1,im dtsfc(i) = rho_a(i) * cp * heat(i) dqsfc(i) = rho_a(i) * hvap * evap(i) diff --git a/physics/satmedmfvdifq.meta b/physics/satmedmfvdifq.meta index fa30cd9f7..88ab676b8 100644 --- a/physics/satmedmfvdifq.meta +++ b/physics/satmedmfvdifq.meta @@ -62,6 +62,13 @@ dimensions = () type = integer intent = in +[progsigma] + standard_name = do_prognostic_updraft_area_fraction + long_name = flag for prognostic sigma in cumuls scheme + units = flag + dimensions = () + type = logical + intent = in [ntrac] standard_name = number_of_vertical_diffusion_tracers long_name = number of tracers to diffuse vertically From e2d5a2a6310a3aa8313ba9c65ff29dc16bf6b954 Mon Sep 17 00:00:00 2001 From: Lisa Bengtsson Date: Wed, 27 Apr 2022 22:40:40 +0000 Subject: [PATCH 08/20] cleaning out some print statements --- physics/progsigma_calc.f90 | 12 ++++++------ physics/samfdeepcnv.f | 36 +++++++++++++++++++++++------------- physics/samfdeepcnv.meta | 15 ++++++++------- 3 files changed, 37 insertions(+), 26 deletions(-) diff --git a/physics/progsigma_calc.f90 b/physics/progsigma_calc.f90 index 55f6b5e3a..21612fd6c 100644 --- a/physics/progsigma_calc.f90 +++ b/physics/progsigma_calc.f90 @@ -43,7 +43,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & integer :: i,k,km1 real(kind=kind_phys) :: termA(im),termB(im),termC(im),termD(im), & mcons(im),fdqa(im),form(im,km), & - qadv(im,km),sigmamax(im),dp(im),inbu(im,km) + qadv(im,km),sigmamax(im),dp(im,km),inbu(im,km) real(kind=kind_phys) :: gcvalmx,epsilon,ZZ,cvg,mcon,buy2, & @@ -82,7 +82,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & do k = 2,km1 do i = 1,im if(cnvflg(i))then - dp(i) = 1000. * del(i,k) + dp(i,k) = 1000. * del(i,k) endif enddo enddo @@ -128,7 +128,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & endif enddo enddo - + !compute termD "The vertical integral of the latent heat convergence is limited to the !buoyant layers with positive moisture convergence (accumulated from the surface). !Lowest level: @@ -140,7 +140,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & do k = 2,km1 do i = 1,im if(cnvflg(i))then - mcon = (hvap*(qadv(i,k)+tmf(i,k)+qmicro(i,k))*dp(i)) + mcon = (hvap*(qadv(i,k)+tmf(i,k)+qmicro(i,k))*dp(i,k)) buy2 = termD(i)+mcon+mcons(i) ! Do the integral over buoyant layers with positive mcon acc from surface if(k > kbcon1(i) .and. k < ktcon(i) .and. buy2 > 0.)then @@ -157,7 +157,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & do k = 2,km1 do i = 1,im if(cnvflg(i))then - tem=(sigmab(i)*zeta(i,k)*inbu(i,k)*dbyo1(i,k))*dp(i) + tem=(sigmab(i)*zeta(i,k)*inbu(i,k)*dbyo1(i,k))*dp(i,k) termA(i)=termA(i)+tem endif enddo @@ -167,7 +167,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & do k = 2,km1 do i = 1,im if(cnvflg(i))then - tem=(dbyo1(i,k)*inbu(i,k))*dp(i) + tem=(dbyo1(i,k)*inbu(i,k))*dp(i,k) termB(i)=termB(i)+tem endif enddo diff --git a/physics/samfdeepcnv.f b/physics/samfdeepcnv.f index 5fd54a2ec..bf48d2035 100644 --- a/physics/samfdeepcnv.f +++ b/physics/samfdeepcnv.f @@ -80,13 +80,13 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & & eps,epsm1,fv,grav,hvap,rd,rv, & & t0c,delt,ntk,ntr,delp, & & prslp,psp,phil,qtr,prevsq,q,q1,t1,u1,v1,fscav, & - & hwrf_samfdeep,progsigma,wclosureflg,cldwrk,rn,kbot,ktop,kcnv, & + & hwrf_samfdeep,progsigma,cldwrk,rn,kbot,ktop,kcnv, & & islimsk,garea,dot,ncloud,hpbl,ud_mf,dd_mf,dt_mf,cnvw,cnvc, & & QLCN, QICN, w_upi, cf_upi, CNV_MFD, & & CNV_DQLDT,CLCN,CNV_FICE,CNV_NDROP,CNV_NICE,mp_phys,mp_phys_mg,& & clam,c0s,c1,betal,betas,evef,pgcon,asolfac, & & do_ca, ca_closure, ca_entr, ca_trigger, nthresh, ca_deep, & - & rainevap, sigmain, sigmaout, & + & rainevap, sigmain, sigmaout, ca_micro, & & errmsg,errflg) ! use machine , only : kind_phys @@ -103,12 +103,12 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & & prslp(:,:), garea(:), hpbl(:), dot(:,:), phil(:,:) real(kind=kind_phys), dimension(:), intent(in) :: fscav logical, intent(in) :: first_time_step,restart,hwrf_samfdeep, & - & progsigma, wclosureflg + & progsigma real(kind=kind_phys), intent(in) :: nthresh real(kind=kind_phys), intent(in) :: ca_deep(:) real(kind=kind_phys), intent(in) :: sigmain(:,:),qmicro(:,:), & & tmf(:,:),q(:,:), prevsq(:,:) - real(kind=kind_phys), intent(out) :: rainevap(:) + real(kind=kind_phys), intent(out) :: rainevap(:), ca_micro(:) real(kind=kind_phys), intent(out) :: sigmaout(:,:) logical, intent(in) :: do_ca,ca_closure,ca_entr,ca_trigger @@ -243,7 +243,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & ! parameter(cinacrmx=-120.,cinacrmn=-120.) parameter(cinacrmx=-120.,cinacrmn=-80.) parameter(bet1=1.875,cd1=.506,f1=2.0,gam1=.5) - parameter(betaw=.03,dxcrtas=8.e3,dxcrtuf=15.e3) + parameter(betaw=.03,dxcrtuf=15.e3) ! ! local variables and arrays @@ -380,6 +380,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & advfac(i) = 0. rainevap(i) = 0. omegac(i)=0. + ca_micro(i)=0. gdx(i) = sqrt(garea(i)) enddo @@ -2456,8 +2457,15 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & c c------- final changed variable per unit mass flux c -!> - If grid size is less than a threshold value (dxcrtas: currently 8km), the quasi-equilibrium assumption of Arakawa-Schubert is not used any longer. +!> - If grid size is less than a threshold value (dxcrtas: currently 8km if progsigma is not used and 30km if progsigma is used), the quasi-equilibrium assumption of Arakawa-Schubert is not used any longer. ! + if(progsigma)then + dxcrtas=30.e3 + else + dxcrtas=8.e3 + endif + + do i = 1, im asqecflg(i) = cnvflg(i) if(asqecflg(i) .and. gdx(i) < dxcrtas) then @@ -2465,13 +2473,6 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & endif enddo -!> - If wclosureflg is true, then quasi-equilibrium closure of Arakawa-Schubert is not used any longer, regardless of resolution - if(wclosureflg)then - do i = 1, im - asqecflg(i) = .false. - enddo - endif - ! !> - If grid size is larger than the threshold value (i.e., asqecflg=.true.), the quasi-equilibrium assumption is used to obtain the cloud base mass flux. To begin with, calculate the change in the temperature and moisture profiles per unit cloud base mass flux. do k = 1, km @@ -2884,6 +2885,15 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & !> - 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 for the grid sizes where the quasi-equilibrium assumption of Arakawa-Schubert is not valid any longer. !! As discussed in Han et al. (2017) \cite han_et_al_2017 , when dtconv is larger than tauadv, the convective mixing is not fully conducted before the cumulus cloud is advected out of the grid cell. In this case, therefore, the cloud base mass flux is further reduced in proportion to the ratio of tauadv to dtconv. + + if(progsigma)then + do i= 1, im + if(cnvflg(i))then + ca_micro(i)=sigmab(i) + endif + enddo + endif + do i= 1, im if(cnvflg(i) .and. .not.asqecflg(i)) then k = kbcon(i) diff --git a/physics/samfdeepcnv.meta b/physics/samfdeepcnv.meta index 2b2942812..5e589b318 100644 --- a/physics/samfdeepcnv.meta +++ b/physics/samfdeepcnv.meta @@ -312,13 +312,6 @@ dimensions = () type = logical intent = in -[wclosureflg] - standard_name = flag_for_wclosure - long_name = flag for vertical velocity closure - units = flag - dimensions = () - type = logical - intent = in [progsigma] standard_name = do_prognostic_updraft_area_fraction long_name = flag for prognostic sigma in cumuls scheme @@ -667,6 +660,14 @@ type = real kind = kind_phys intent = out +[ca_micro] + standard_name = output_prognostic_sigma_two + long_name = output of prognostic area fraction two + units = frac + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP From 8b815e026ddcd8e660e437dc94bb89058e6f80de Mon Sep 17 00:00:00 2001 From: Lisa Bengtsson Date: Fri, 29 Apr 2022 03:03:46 +0000 Subject: [PATCH 09/20] address some bugs caught by debug flag --- physics/progsigma_calc.f90 | 11 ++++++----- physics/samfdeepcnv.f | 34 +++++++++++++++++----------------- physics/samfdeepcnv.meta | 8 -------- physics/samfshalcnv.f | 9 ++++++--- 4 files changed, 29 insertions(+), 33 deletions(-) diff --git a/physics/progsigma_calc.f90 b/physics/progsigma_calc.f90 index 21612fd6c..c05af3003 100644 --- a/physics/progsigma_calc.f90 +++ b/physics/progsigma_calc.f90 @@ -26,8 +26,8 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & ! intent in integer, intent(in) :: im,km,kbcon1(im),ktcon(im) - real, intent(in) :: hvap,delt - real, intent(in) :: prevsq(im,km), q(im,km),del(im,km), & + real(kind=kind_phys), intent(in) :: hvap,delt + real(kind=kind_phys), intent(in) :: prevsq(im,km), q(im,km),del(im,km), & qmicro(im,km),tmf(im,km),dbyo1(im,km),zdqca(im,km), & omega_u(im,km),zeta(im,km),gdx(im) logical, intent(in) :: flag_init,flag_restart,cnvflg(im),flag_shallow @@ -63,7 +63,8 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & do i = 1,im sigmaout(i,k)=0. inbu(i,k)=0. - form(i,k)=0. + form(i,k)=0. + dp(i,k)=0. enddo enddo @@ -157,7 +158,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & do k = 2,km1 do i = 1,im if(cnvflg(i))then - tem=(sigmab(i)*zeta(i,k)*inbu(i,k)*dbyo1(i,k))*dp(i,k) + tem=sigmab(i)*zeta(i,k)*inbu(i,k)*dbyo1(i,k)*dp(i,k) termA(i)=termA(i)+tem endif enddo @@ -167,7 +168,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & do k = 2,km1 do i = 1,im if(cnvflg(i))then - tem=(dbyo1(i,k)*inbu(i,k))*dp(i,k) + tem=zeta(i,k)*dbyo1(i,k)*inbu(i,k)*dp(i,k) termB(i)=termB(i)+tem endif enddo diff --git a/physics/samfdeepcnv.f b/physics/samfdeepcnv.f index bf48d2035..071bf0557 100644 --- a/physics/samfdeepcnv.f +++ b/physics/samfdeepcnv.f @@ -86,8 +86,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & & CNV_DQLDT,CLCN,CNV_FICE,CNV_NDROP,CNV_NICE,mp_phys,mp_phys_mg,& & clam,c0s,c1,betal,betas,evef,pgcon,asolfac, & & do_ca, ca_closure, ca_entr, ca_trigger, nthresh, ca_deep, & - & rainevap, sigmain, sigmaout, ca_micro, & - & errmsg,errflg) + & rainevap, sigmain, sigmaout, errmsg,errflg) ! use machine , only : kind_phys use funcphys , only : fpvs @@ -108,7 +107,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & real(kind=kind_phys), intent(in) :: ca_deep(:) real(kind=kind_phys), intent(in) :: sigmain(:,:),qmicro(:,:), & & tmf(:,:),q(:,:), prevsq(:,:) - real(kind=kind_phys), intent(out) :: rainevap(:), ca_micro(:) + real(kind=kind_phys), intent(out) :: rainevap(:) real(kind=kind_phys), intent(out) :: sigmaout(:,:) logical, intent(in) :: do_ca,ca_closure,ca_entr,ca_trigger @@ -380,7 +379,6 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & advfac(i) = 0. rainevap(i) = 0. omegac(i)=0. - ca_micro(i)=0. gdx(i) = sqrt(garea(i)) enddo @@ -583,14 +581,20 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & buo(i,k) = 0. drag(i,k) = 0. cnvwt(i,k)= 0. + endif + enddo + enddo + + do k = 1, km + do i = 1, im dbyo1(i,k)=0. zdqca(i,k)=0. qlks(i,k)=0. omega_u(i,k)=0. zeta(i,k)=1.0 - endif - enddo + enddo enddo + ! ! initialize tracer variables ! @@ -1811,9 +1815,13 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & do i = 1, im if (cnvflg(i)) then if(k >= kbcon1(i) .and. k < ktcon(i)) then - zeta(i,k)=eta(i,k)*(omegac(i)/omega_u(i,k)) - zeta(i,k)=MAX(0.,zeta(i,k)) - zeta(i,k)=MIN(1.,zeta(i,k)) + if(omega_u(i,k) .ne. 0.)then + zeta(i,k)=eta(i,k)*(omegac(i)/omega_u(i,k)) + else + zeta(i,k)=0. + endif + zeta(i,k)=MAX(0.,zeta(i,k)) + zeta(i,k)=MIN(1.,zeta(i,k)) endif endif enddo @@ -2886,14 +2894,6 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & !> - 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 for the grid sizes where the quasi-equilibrium assumption of Arakawa-Schubert is not valid any longer. !! As discussed in Han et al. (2017) \cite han_et_al_2017 , when dtconv is larger than tauadv, the convective mixing is not fully conducted before the cumulus cloud is advected out of the grid cell. In this case, therefore, the cloud base mass flux is further reduced in proportion to the ratio of tauadv to dtconv. - if(progsigma)then - do i= 1, im - if(cnvflg(i))then - ca_micro(i)=sigmab(i) - endif - enddo - endif - do i= 1, im if(cnvflg(i) .and. .not.asqecflg(i)) then k = kbcon(i) diff --git a/physics/samfdeepcnv.meta b/physics/samfdeepcnv.meta index 5e589b318..3f28035b6 100644 --- a/physics/samfdeepcnv.meta +++ b/physics/samfdeepcnv.meta @@ -660,14 +660,6 @@ type = real kind = kind_phys intent = out -[ca_micro] - standard_name = output_prognostic_sigma_two - long_name = output of prognostic area fraction two - units = frac - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = out [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP diff --git a/physics/samfshalcnv.f b/physics/samfshalcnv.f index 325566877..ef0366b84 100644 --- a/physics/samfshalcnv.f +++ b/physics/samfshalcnv.f @@ -513,7 +513,6 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & enddo enddo - do i = 1,im omegac(i)=0. enddo @@ -1551,12 +1550,16 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & endif enddo c -! > - Calculate the mean updraft velocity within the cloud (omega). +c Compute zeta for prog closure do k = 2, km1 do i = 1, im if (cnvflg(i)) then if(k > kbcon1(i) .and. k < ktcon(i)) then - zeta(i,k)=eta(i,k)*(omegac(i)/omega_u(i,k)) + if(omega_u(i,k) .ne. 0.)then + zeta(i,k)=eta(i,k)*(omegac(i)/omega_u(i,k)) + else + zeta(i,k)=0. + endif zeta(i,k)=MAX(0.,zeta(i,k)) zeta(i,k)=MIN(1.,zeta(i,k)) endif From 527e1b976bd74dc0214a13f91f804ec2334d862c Mon Sep 17 00:00:00 2001 From: "samuel.trahan" Date: Mon, 2 May 2022 22:11:42 +0000 Subject: [PATCH 10/20] Pass -DCCPP_SINGLE_PRECISION from cmake to -DSINGLE_PREC in cpp --- CMakeLists.txt | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index 60531b9a5..691b283f2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -29,6 +29,13 @@ if(CMAKE_BUILD_TYPE STREQUAL "Debug") add_definitions(-DDEBUG) endif() +if(CCPP_SINGLE_PREC) + message(STATUS "CCPP Single Precision Mode activated.") + add_definitions(SINGLE_PREC) +else(CCPP_SINGLE_PREC) + message(STATUS "CCPP Double Precision Mode activated.") +endif(CCPP_SINGLE_PREC) + #------------------------------------------------------------------------------ # Request a static build option(BUILD_SHARED_LIBS "Build a shared library" OFF) From 6871a936a9df8054fa2b4b34c6e52ad5d5cce738 Mon Sep 17 00:00:00 2001 From: Samuel Trahan Date: Wed, 4 May 2022 17:32:24 +0000 Subject: [PATCH 11/20] Changes needed for 32-bit physics --- CMakeLists.txt | 7 ---- physics/GFS_rrtmgp_cloud_overlap.F90 | 4 +-- physics/GFS_suite_interstitial_4.F90 | 10 +++--- physics/calpreciptype.f90 | 9 ++++- physics/machine.F | 2 +- physics/maximum_hourly_diagnostics.F90 | 14 ++++---- physics/mersenne_twister.f | 46 ++++++++++++++------------ physics/module_sf_mynn.F90 | 4 +-- physics/module_sf_noahmplsm.f90 | 4 +-- physics/module_sf_ruclsm.F90 | 5 +-- physics/module_soil_pre.F90 | 24 ++++++++------ physics/radiation_gases.f | 2 +- physics/radlw_main.meta | 2 +- physics/radsw_main.F90 | 2 +- physics/radsw_main.meta | 2 +- 15 files changed, 72 insertions(+), 65 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 691b283f2..60531b9a5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -29,13 +29,6 @@ if(CMAKE_BUILD_TYPE STREQUAL "Debug") add_definitions(-DDEBUG) endif() -if(CCPP_SINGLE_PREC) - message(STATUS "CCPP Single Precision Mode activated.") - add_definitions(SINGLE_PREC) -else(CCPP_SINGLE_PREC) - message(STATUS "CCPP Double Precision Mode activated.") -endif(CCPP_SINGLE_PREC) - #------------------------------------------------------------------------------ # Request a static build option(BUILD_SHARED_LIBS "Build a shared library" OFF) diff --git a/physics/GFS_rrtmgp_cloud_overlap.F90 b/physics/GFS_rrtmgp_cloud_overlap.F90 index 13794641b..c1a6c4763 100644 --- a/physics/GFS_rrtmgp_cloud_overlap.F90 +++ b/physics/GFS_rrtmgp_cloud_overlap.F90 @@ -99,7 +99,7 @@ subroutine GFS_rrtmgp_cloud_overlap_run(nCol, nLev, yearlen, doSWrad, doLWrad, ! Cloud overlap parameter ! if (iovr == iovr_dcorr .or. iovr == iovr_exp .or. iovr == iovr_exprand) then - call get_alpha_exper(nCol, nLev, iovr, iovr_exprand, deltaZc*0.001, de_lgth, cld_frac, cloud_overlap_param) + call get_alpha_exper(nCol, nLev, iovr, iovr_exprand, deltaZc*0.001_kind_phys, de_lgth, cld_frac, cloud_overlap_param) else de_lgth(:) = 0. cloud_overlap_param(:,:) = 0. @@ -110,7 +110,7 @@ subroutine GFS_rrtmgp_cloud_overlap_run(nCol, nLev, yearlen, doSWrad, doLWrad, ! if (imfdeepcnv == imfdeepcnv_samf .or. imfdeepcnv == imfdeepcnv_gf) then if (iovr_convcld == iovr_dcorr .or. iovr_convcld == iovr_exp .or. iovr_convcld == iovr_exprand) then - call get_alpha_exper(nCol, nLev, iovr_convcld, iovr_exprand, deltaZc*0.001, de_lgth, cld_cnv_frac, cnv_cloud_overlap_param) + call get_alpha_exper(nCol, nLev, iovr_convcld, iovr_exprand, deltaZc*0.001_kind_phys, de_lgth, cld_cnv_frac, cnv_cloud_overlap_param) else de_lgth(:) = 0. cnv_cloud_overlap_param(:,:) = 0. diff --git a/physics/GFS_suite_interstitial_4.F90 b/physics/GFS_suite_interstitial_4.F90 index cbabb991b..18fcfda09 100644 --- a/physics/GFS_suite_interstitial_4.F90 +++ b/physics/GFS_suite_interstitial_4.F90 @@ -224,7 +224,7 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, tracers_total, ntr qc_mp(i,k) = (clw(i,k,2)-save_qc(i,k)) / (one-spechum(i,k)) !> - Convert number concentration from moist to dry nc_mp(i,k) = gq0(i,k,ntlnc) / (one-spechum(i,k)) - nc_mp(i,k) = max(zero, nc_mp(i,k) + make_DropletNumber(qc_mp(i,k) * rho, nwfa(i,k)*rho) * orho) + nc_mp(i,k) = max(zero, nc_mp(i,k) + make_DropletNumber(real(qc_mp(i,k) * rho), real(nwfa(i,k)*rho)) * orho) !> - Convert number concentrations from dry to moist gq0(i,k,ntlnc) = nc_mp(i,k) / (one+qv_mp(i,k)) endif @@ -233,7 +233,7 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, tracers_total, ntr qi_mp(i,k) = (clw(i,k,1)-save_qi(i,k)) / (one-spechum(i,k)) !> - Convert number concentration from moist to dry ni_mp(i,k) = gq0(i,k,ntinc) / (one-spechum(i,k)) - ni_mp(i,k) = max(zero, ni_mp(i,k) + make_IceNumber(qi_mp(i,k) * rho, save_tcp(i,k)) * orho) + ni_mp(i,k) = max(zero, ni_mp(i,k) + make_IceNumber(real(qi_mp(i,k) * rho), real(save_tcp(i,k)) * orho)) !> - Convert number concentrations from dry to moist gq0(i,k,ntinc) = ni_mp(i,k) / (one+qv_mp(i,k)) endif @@ -249,13 +249,13 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, tracers_total, ntr !> - Update cloud water mixing ratio qc_mp(i,k) = (clw(i,k,2)-save_qc(i,k)) !> - Update cloud water number concentration - gq0(i,k,ntlnc) = max(zero, gq0(i,k,ntlnc) + make_DropletNumber(qc_mp(i,k) * rho, nwfa(i,k)*rho) * orho) + gq0(i,k,ntlnc) = max(zero, gq0(i,k,ntlnc) + make_DropletNumber(real(qc_mp(i,k) * rho), real(nwfa(i,k)*rho)) * orho) endif if (ntinc>0) then !> - Update cloud ice mixing ratio qi_mp(i,k) = (clw(i,k,1)-save_qi(i,k)) !> - Update cloud ice number concentration - gq0(i,k,ntinc) = max(zero, gq0(i,k,ntinc) + make_IceNumber(qi_mp(i,k) * rho, save_tcp(i,k)) * orho) + gq0(i,k,ntinc) = max(zero, gq0(i,k,ntinc) + make_IceNumber(real(qi_mp(i,k) * rho), real(save_tcp(i,k))) * orho) endif enddo enddo @@ -290,4 +290,4 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, tracers_total, ntr end subroutine GFS_suite_interstitial_4_run - end module GFS_suite_interstitial_4 \ No newline at end of file + end module GFS_suite_interstitial_4 diff --git a/physics/calpreciptype.f90 b/physics/calpreciptype.f90 index d3fbb253b..956ed8c55 100644 --- a/physics/calpreciptype.f90 +++ b/physics/calpreciptype.f90 @@ -509,7 +509,14 @@ subroutine calwxt_ramer(lm,lp1,t,q,pmid,rh,td,pint,ptyp) real(kind=kind_phys) rhmax,twmax,ptop,dpdrh,twtop,rhtop,wgt1,wgt2, & rhavg,dtavg,dpk,ptw,pbot ! real(kind=kind_phys) b,qtmp,rate,qc -! real(kind=kind_phys),external :: xmytw (now inside the module) +! + interface + function xmytw(t,td,p) + use machine , only : kind_phys + implicit none + real(kind=kind_phys) t, td, p, xmytw + end function xmytw + end interface ! ! initialize. icefrac = -9999. diff --git a/physics/machine.F b/physics/machine.F index 2ee7fb865..9b09d235c 100644 --- a/physics/machine.F +++ b/physics/machine.F @@ -33,7 +33,7 @@ module machine # endif &, kind_rad = 4 & &, kind_phys = 4 ,kind_taum=4 & - &, kind_grid = 4 & + &, kind_grid = 8 &! atmos_cubed_sphere requres kind_grid=8 &, kind_REAL = 4 &! used in cmp_comm &, kind_LOGICAL = 4 & &, kind_INTEGER = 4 ! -,,- diff --git a/physics/maximum_hourly_diagnostics.F90 b/physics/maximum_hourly_diagnostics.F90 index 6beae0da2..ddbff5725 100644 --- a/physics/maximum_hourly_diagnostics.F90 +++ b/physics/maximum_hourly_diagnostics.F90 @@ -144,11 +144,11 @@ subroutine max_fields(phil,ref3D,grav,im,levs,refd,tk,refd263k) real (kind=kind_phys), intent(in) :: grav real (kind=kind_phys), intent(in),dimension(:,:) :: phil,ref3D,tk integer :: i,k,ll,ipt,kpt - real :: dbz1avg,zmidp1,zmidloc,refl,fact - real, dimension(im,levs) :: z - real, dimension(im) :: zintsfc - real, dimension(:), intent(inout) :: refd,refd263k - REAL :: dbz1(2),dbzk,dbzk1 + real(kind_phys) :: dbz1avg,zmidp1,zmidloc,refl,fact + real(kind_phys), dimension(im,levs) :: z + real(kind_phys), dimension(im) :: zintsfc + real(kind_phys), dimension(:), intent(inout) :: refd,refd263k + REAL(kind_phys) :: dbz1(2),dbzk,dbzk1 logical :: counter do i=1,im do k=1,levs @@ -185,7 +185,7 @@ subroutine max_fields(phil,ref3D,grav,im,levs,refd,tk,refd263k) dbz1avg=dbz1(2)+(dbz1(2)-dbz1(1))*fact !-- Convert to dBZ (10*logZ) as the last step if (dbz1avg>0.01) then - dbz1avg=10.*alog10(dbz1avg) + dbz1avg=10.*log10(dbz1avg) else dbz1avg=-35. endif @@ -214,7 +214,7 @@ subroutine max_fields(phil,ref3D,grav,im,levs,refd,tk,refd263k) dbz1avg=maxval(dbz1) !-- Convert to dBZ (10*logZ) as the last step if (dbz1avg>0.01) then - dbz1avg=10.*alog10(dbz1avg) + dbz1avg=10.*log10(dbz1avg) else dbz1avg=-35. endif diff --git a/physics/mersenne_twister.f b/physics/mersenne_twister.f index 8cc6bd5e5..58bf43487 100644 --- a/physics/mersenne_twister.f +++ b/physics/mersenne_twister.f @@ -160,6 +160,7 @@ ! !$$$ module mersenne_twister + use machine, only: kind_dbl_prec private ! Public declarations public random_stat @@ -188,7 +189,7 @@ module mersenne_twister integer:: mti=n+1 integer:: mt(0:n-1) integer:: iset - real:: gset + real(kind_dbl_prec):: gset end type ! Saved data type(random_stat),save:: sstat @@ -300,8 +301,8 @@ subroutine random_setseed_t(inseed,stat) !> This function generates random numbers in functional mode. function random_number_f() result(harvest) implicit none - real:: harvest - real h(1) + real(kind_dbl_prec):: harvest + real(kind_dbl_prec) :: h(1) if(sstat%mti.eq.n+1) call random_setseed_t(iseed,sstat) call random_number_t(h,sstat) harvest=h(1) @@ -310,7 +311,7 @@ function random_number_f() result(harvest) !> This subroutine generates random numbers in interactive mode. subroutine random_number_i(harvest,inseed) implicit none - real,intent(out):: harvest(:) + real(kind_dbl_prec),intent(out):: harvest(:) integer,intent(in):: inseed type(random_stat) stat call random_setseed_t(inseed,stat) @@ -320,7 +321,7 @@ subroutine random_number_i(harvest,inseed) !> This subroutine generates random numbers in saved mode; overloads Fortran 90 standard. subroutine random_number_s(harvest) implicit none - real,intent(out):: harvest(:) + real(kind_dbl_prec),intent(out):: harvest(:) if(sstat%mti.eq.n+1) call random_setseed_t(iseed,sstat) call random_number_t(harvest,sstat) end subroutine @@ -328,7 +329,7 @@ subroutine random_number_s(harvest) !> This subroutine generates random numbers in thread-safe mode. subroutine random_number_t(harvest,stat) implicit none - real,intent(out):: harvest(:) + real(kind_dbl_prec),intent(out):: harvest(:) type(random_stat),intent(inout):: stat integer j,kk,y integer tshftu,tshfts,tshftt,tshftl @@ -359,9 +360,12 @@ subroutine random_number_t(harvest,stat) y=ieor(y,iand(tshftt(y),tmaskc)) y=ieor(y,tshftl(y)) if(y.lt.0) then - harvest(j)=(real(y)+2.0**32)/(2.0**32-1.0) + harvest(j)=(real(y,kind=kind_dbl_prec)+ & + & 2.0_kind_dbl_prec**32)/ & + & (2.0_kind_dbl_prec**32-1.0_kind_dbl_prec) else - harvest(j)=real(y)/(2.0**32-1.0) + harvest(j)=real(y)/(2.0_kind_dbl_prec**32- & + & 1.0_kind_dbl_prec) endif stat%mti=stat%mti+1 enddo @@ -370,8 +374,8 @@ subroutine random_number_t(harvest,stat) !> This subrouitne generates Gaussian random numbers in functional mode. function random_gauss_f() result(harvest) implicit none - real:: harvest - real h(1) + real(kind_dbl_prec):: harvest + real(kind_dbl_prec) :: h(1) if(sstat%mti.eq.n+1) call random_setseed_t(iseed,sstat) call random_gauss_t(h,sstat) harvest=h(1) @@ -380,7 +384,7 @@ function random_gauss_f() result(harvest) !> This subrouitne generates Gaussian random numbers in interactive mode. subroutine random_gauss_i(harvest,inseed) implicit none - real,intent(out):: harvest(:) + real(kind_dbl_prec),intent(out):: harvest(:) integer,intent(in):: inseed type(random_stat) stat call random_setseed_t(inseed,stat) @@ -390,7 +394,7 @@ subroutine random_gauss_i(harvest,inseed) !> This subroutine generates Gaussian random numbers in saved mode. subroutine random_gauss_s(harvest) implicit none - real,intent(out):: harvest(:) + real(kind_dbl_prec),intent(out):: harvest(:) if(sstat%mti.eq.n+1) call random_setseed_t(iseed,sstat) call random_gauss_t(harvest,sstat) end subroutine @@ -398,10 +402,10 @@ subroutine random_gauss_s(harvest) !> This subroutine generates Gaussian random numbers in thread-safe mode. subroutine random_gauss_t(harvest,stat) implicit none - real,intent(out):: harvest(:) + real(kind_dbl_prec),intent(out):: harvest(:) type(random_stat),intent(inout):: stat integer mx,my,mz,j - real r2(2),r,g1,g2 + real(kind_dbl_prec) :: r2(2),r,g1,g2 mz=size(harvest) if(mz.le.0) return mx=0 @@ -436,14 +440,14 @@ subroutine random_gauss_t(harvest,stat) contains !> This subroutine contains numerical Recipes algorithm to generate Gaussian random numbers. subroutine rgauss(r1,r2,r,g1,g2) - real,intent(in):: r1,r2 - real,intent(out):: r,g1,g2 - real v1,v2,fac - v1=2.*r1-1. - v2=2.*r2-1. + real(kind_dbl_prec),intent(in):: r1,r2 + real(kind_dbl_prec),intent(out):: r,g1,g2 + real(kind_dbl_prec) :: v1,v2,fac + v1=2._kind_dbl_prec*r1-1._kind_dbl_prec + v2=2._kind_dbl_prec*r2-1._kind_dbl_prec r=v1**2+v2**2 if(r.lt.1.) then - fac=sqrt(-2.*log(r)/r) + fac=sqrt(-2._kind_dbl_prec*log(r)/r) g1=v1*fac g2=v2*fac endif @@ -489,7 +493,7 @@ subroutine random_index_t(imax,iharvest,stat) type(random_stat),intent(inout):: stat integer,parameter:: mh=n integer i1,i2,mz - real h(mh) + real(kind_dbl_prec) :: h(mh) mz=size(iharvest) do i1=1,mz,mh i2=min((i1-1)+mh,mz) diff --git a/physics/module_sf_mynn.F90 b/physics/module_sf_mynn.F90 index 5f227750a..bc874ace6 100644 --- a/physics/module_sf_mynn.F90 +++ b/physics/module_sf_mynn.F90 @@ -2804,8 +2804,8 @@ SUBROUTINE znot_m_v6(uref, znotm) ! znotm(meter): areodynamical roughness scale over water ! - REAL(kind=kind_phys), INTENT(IN) :: uref - REAL(kind=kind_phys), INTENT(OUT):: znotm + REAL, INTENT(IN) :: uref + REAL, INTENT(OUT):: znotm real(kind=kind_phys), parameter :: p13 = -1.296521881682694e-02,& & p12 = 2.855780863283819e-01, p11 = -1.597898515251717e+00,& & p10 = -8.396975715683501e+00, & diff --git a/physics/module_sf_noahmplsm.f90 b/physics/module_sf_noahmplsm.f90 index 1c899e4bd..61b92990b 100644 --- a/physics/module_sf_noahmplsm.f90 +++ b/physics/module_sf_noahmplsm.f90 @@ -681,7 +681,7 @@ subroutine noahmp_sflx (parameters, & logical :: dveg_active !< flag to run dynamic vegetation logical :: crop_active !< flag to run crop model ! add canopy heat storage (C.He added based on GY Niu's communication) - real :: canhs ! canopy heat storage change w/m2 + real (kind=kind_phys) :: canhs ! canopy heat storage change w/m2 ! maximum lai/sai used for some parameterizations based on plant growthi @@ -4494,7 +4494,7 @@ subroutine bare_flux (parameters,nsnow ,nsoil ,isnow ,dt ,sag , & if(opt_sfc == 3) then call sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur , & !in zlvl ,tgb ,thsfc_loc,prslkix,prsik1x ,prslk1x ,z0m , & !in - zpd ,snowh,shdfac ,garea1 ,.false. ,0.0,ivgtyp , & !in + zpd ,snowh,shdfac ,garea1 ,.false. ,0.0_kind_phys,ivgtyp , & !in ustarx ,fm ,fh ,fm2 ,fh2 , & !inout z0h ,fv ,csigmaf0,cm ,ch ) !out diff --git a/physics/module_sf_ruclsm.F90 b/physics/module_sf_ruclsm.F90 index b39610bc8..a27d0f287 100644 --- a/physics/module_sf_ruclsm.F90 +++ b/physics/module_sf_ruclsm.F90 @@ -7603,10 +7603,11 @@ END SUBROUTINE SOILIN !>\ingroup lsm_ruc_group !> This function calculates the liquid saturation vapor mixing ratio as !! a function of temperature and pressure (from Thompson scheme). - REAL FUNCTION RSLF(P,T) + FUNCTION RSLF(P,T) IMPLICIT NONE - REAL, INTENT(IN):: P, T + REAL(kind_phys), INTENT(IN):: P, T + REAL(kind_phys) :: RSLF REAL:: ESL,X REAL, PARAMETER:: C0= .611583699E03 REAL, PARAMETER:: C1= .444606896E02 diff --git a/physics/module_soil_pre.F90 b/physics/module_soil_pre.F90 index 8eb5a5775..149f87a1c 100644 --- a/physics/module_soil_pre.F90 +++ b/physics/module_soil_pre.F90 @@ -5,6 +5,8 @@ module module_soil_pre !tgs Initialize RUC LSM levels, soil temp/moisture + use machine, only: kind_phys + implicit none private @@ -26,8 +28,8 @@ SUBROUTINE init_soil_depth_3 ( zs , dzs , num_soil_levels ) INTEGER, INTENT(IN) :: num_soil_levels - REAL, DIMENSION(1:num_soil_levels), INTENT(OUT) :: zs, dzs - REAL, DIMENSION(1:num_soil_levels) :: zs2 + REAL(kind_phys), DIMENSION(1:num_soil_levels), INTENT(OUT) :: zs, dzs + REAL(kind_phys), DIMENSION(1:num_soil_levels) :: zs2 INTEGER :: l @@ -90,21 +92,21 @@ SUBROUTINE init_soil_3_real ( tsk , tmn , smois , tslb , & INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input - REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input - REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input - REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst + REAL(kind_phys) , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input + REAL(kind_phys) , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input + REAL(kind_phys) , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst - REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn - REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk - REAL , DIMENSION(num_soil_layers) :: zs , dzs + REAL(kind_phys) , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn + REAL(kind_phys) , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk + REAL(kind_phys) , DIMENSION(num_soil_layers) :: zs , dzs - REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois + REAL(kind_phys) , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois - REAL , ALLOCATABLE , DIMENSION(:) :: zhave + REAL(kind_phys) , ALLOCATABLE , DIMENSION(:) :: zhave logical :: debug_print = .false. INTEGER :: i , j , l , lout , lin , lwant , lhave, k - REAL :: temp + REAL(kind_phys) :: temp ! Allocate the soil layer array used for interpolating. diff --git a/physics/radiation_gases.f b/physics/radiation_gases.f index 157da8e09..d6f1d7259 100644 --- a/physics/radiation_gases.f +++ b/physics/radiation_gases.f @@ -371,7 +371,7 @@ subroutine gas_init & endif do k = 1, LOZ - pkstr(k) = fpkapx(pstr(k)*100.0) + pkstr(k) = fpkapx(pstr(k)*100.0_kind_phys) enddo endif ! end if_ioznflg_block diff --git a/physics/radlw_main.meta b/physics/radlw_main.meta index df1a368c5..9286c45cb 100644 --- a/physics/radlw_main.meta +++ b/physics/radlw_main.meta @@ -1,7 +1,7 @@ [ccpp-table-properties] name = rrtmg_lw type = scheme - dependencies = mersenne_twister.f,physcons.F90,physparam.f,radlw_datatb.f,radlw_param.f + dependencies = machine.F,mersenne_twister.f,physcons.F90,physparam.f,radlw_datatb.f,radlw_param.f ######################################################################## [ccpp-arg-table] diff --git a/physics/radsw_main.F90 b/physics/radsw_main.F90 index ae2fa7ad5..5d7d62dcc 100644 --- a/physics/radsw_main.F90 +++ b/physics/radsw_main.F90 @@ -2040,7 +2040,7 @@ subroutine mcica_subcol & real (kind=kind_phys) :: cdfunc(nlay,ngptsw), tem1, & & fac_lcf(nlay), & & cdfun2(nlay,ngptsw) - real (kind=kind_dbl_prec) :: rand2d(nlay*ngptsw), rand1d(ngptsw) + real (kind=kind_dbl_prec) :: rand2d(nlay*ngptsw), rand1d(ngptsw) ! must be default real kind to match mersenne twister code type (random_stat) :: stat ! for thread safe random generator diff --git a/physics/radsw_main.meta b/physics/radsw_main.meta index 70bc136f3..506e2edf0 100644 --- a/physics/radsw_main.meta +++ b/physics/radsw_main.meta @@ -1,7 +1,7 @@ [ccpp-table-properties] name = rrtmg_sw type = scheme - dependencies = mersenne_twister.f,physcons.F90,physparam.f,radsw_datatb.f,radsw_param.f + dependencies = machine.F,mersenne_twister.f,physcons.F90,physparam.f,radsw_datatb.f,radsw_param.f ######################################################################## [ccpp-arg-table] From e7c42c7d57740d6f8c3852ce3d9dfbab720e6c86 Mon Sep 17 00:00:00 2001 From: "samuel.trahan" Date: Thu, 5 May 2022 00:43:21 +0000 Subject: [PATCH 12/20] Move some code to modules --- physics/GFS_MP_generic_post.F90 | 2 +- physics/calpreciptype.f90 | 10 +++------- physics/cires_orowam2017.f | 3 +++ physics/cires_ugwp.F90 | 5 +++++ physics/cires_ugwp_triggers.F90 | 3 +++ physics/cires_ugwpv1_oro.F90 | 2 +- physics/cires_ugwpv1_sporo.F90 | 4 +++- physics/hedmf.f | 3 ++- physics/lsm_noah.f | 1 + physics/mfpbl.f | 4 +++- physics/mfpblt.f | 4 +++- physics/mfpbltq.f | 4 +++- physics/mfscu.f | 4 +++- physics/mfscuq.f | 4 +++- physics/module_bl_mynn.F90 | 3 +++ physics/moninshoc.f | 3 +++ physics/satmedmfvdif.F | 4 ++++ physics/satmedmfvdifq.F | 4 +++- physics/sflx.f | 12 +++++++++++- physics/tridi.f | 4 +++- physics/ugwp_driver_v0.F | 5 ++++- physics/unified_ugwp.F90 | 3 ++- 22 files changed, 70 insertions(+), 21 deletions(-) diff --git a/physics/GFS_MP_generic_post.F90 b/physics/GFS_MP_generic_post.F90 index a7be0ab4c..97deec10f 100644 --- a/physics/GFS_MP_generic_post.F90 +++ b/physics/GFS_MP_generic_post.F90 @@ -30,7 +30,7 @@ subroutine GFS_MP_generic_post_run( index_of_process_dfi_radar, ix_dfi_radar, dfi_radar_tten, radar_tten_limits, fhour, errmsg, errflg) ! use machine, only: kind_phys - + use calpreciptype_mod, only: calpreciptype implicit none integer, intent(in) :: im, levs, kdt, nrcm, nncl, ntcw, ntrac, num_dfi_radar, index_of_process_dfi_radar diff --git a/physics/calpreciptype.f90 b/physics/calpreciptype.f90 index 956ed8c55..54e8fa2b9 100644 --- a/physics/calpreciptype.f90 +++ b/physics/calpreciptype.f90 @@ -1,6 +1,8 @@ !>\file calpreciptype.f90 !! This file contains the subroutines that calculates dominant precipitation type. +module calpreciptype_mod +contains !>\ingroup gfs_calpreciptype !! Foure algorithms are called to calculate dominant precipitation type, and the !!tallies are sumed in calwxt_dominant(). @@ -510,13 +512,6 @@ subroutine calwxt_ramer(lm,lp1,t,q,pmid,rh,td,pint,ptyp) rhavg,dtavg,dpk,ptw,pbot ! real(kind=kind_phys) b,qtmp,rate,qc ! - interface - function xmytw(t,td,p) - use machine , only : kind_phys - implicit none - real(kind=kind_phys) t, td, p, xmytw - end function xmytw - end interface ! ! initialize. icefrac = -9999. @@ -1391,3 +1386,4 @@ subroutine calwxt_dominant(nalg,rain,freezr,sleet,snow, & return end !! @} +end module calpreciptype_mod diff --git a/physics/cires_orowam2017.f b/physics/cires_orowam2017.f index c20f98f42..ae5f355d3 100644 --- a/physics/cires_orowam2017.f +++ b/physics/cires_orowam2017.f @@ -1,3 +1,5 @@ + module cires_orowam2017 + contains subroutine oro_wam_2017(im, levs,npt,ipt, kref,kdt,me,master, & dtp,dxres, taub, u1, v1, t1, xn, yn, bn2, rho, prsi, prsL, & del, sigma, hprime, gamma, theta, @@ -384,3 +386,4 @@ subroutine ugwpv0_tofd1d(levs, sigflt, elvmax, zsurf, enddo ! end subroutine ugwpv0_tofd1d + end module cires_orowam2017 diff --git a/physics/cires_ugwp.F90 b/physics/cires_ugwp.F90 index c4f0a255d..2d8eafc19 100644 --- a/physics/cires_ugwp.F90 +++ b/physics/cires_ugwp.F90 @@ -16,9 +16,14 @@ module cires_ugwp use machine, only: kind_phys use cires_ugwpv0_module, only: knob_ugwp_version, cires_ugwpv0_mod_init, cires_ugwpv0_mod_finalize + use ugwp_driver_v0 use gwdps, only: gwdps_run + use cires_ugwp_triggers + + use ugwp_driver_v0 + implicit none private diff --git a/physics/cires_ugwp_triggers.F90 b/physics/cires_ugwp_triggers.F90 index 4a8b97590..82f762c56 100644 --- a/physics/cires_ugwp_triggers.F90 +++ b/physics/cires_ugwp_triggers.F90 @@ -1,3 +1,5 @@ + module cires_ugwp_triggers + contains ! subroutine slat_geos5_tamp_v0(im, tau_amp, xlatdeg, tau_gw) !================= @@ -97,3 +99,4 @@ subroutine init_nazdir_v0(naz, xaz, yaz) yaz(4) =-1.0 !S endif end subroutine init_nazdir_v0 + end module cires_ugwp_triggers diff --git a/physics/cires_ugwpv1_oro.F90 b/physics/cires_ugwpv1_oro.F90 index 247112bf1..66d0e472c 100644 --- a/physics/cires_ugwpv1_oro.F90 +++ b/physics/cires_ugwpv1_oro.F90 @@ -1,5 +1,5 @@ module cires_ugwpv1_oro - + use cires_ugwpv1_sporo contains subroutine orogw_v1 (im, km, imx, me, master, dtp, kdt, do_tofd, & diff --git a/physics/cires_ugwpv1_sporo.F90 b/physics/cires_ugwpv1_sporo.F90 index c840b49d8..fbd3eaa0b 100644 --- a/physics/cires_ugwpv1_sporo.F90 +++ b/physics/cires_ugwpv1_sporo.F90 @@ -1,4 +1,5 @@ - + module cires_ugwpv1_sporo + contains subroutine oro_spectral_solver(im, levs,npt,ipt, kref,kdt,me,master, & dtp,dxres, taub, u1, v1, t1, xn, yn, bn2, rho, prsi, prsL, & del, sigma, hprime, gamma, theta, & @@ -349,3 +350,4 @@ subroutine oro_meanflow(nz, nzi, u1, v1, t1, pint, pmid, & end subroutine oro_meanflow + end module cires_ugwpv1_sporo diff --git a/physics/hedmf.f b/physics/hedmf.f index 83d0fe1b0..604483e53 100644 --- a/physics/hedmf.f +++ b/physics/hedmf.f @@ -5,7 +5,8 @@ !> This module contains the CCPP-compliant hybrid eddy-diffusivity mass-flux !! scheme. module hedmf - + use tridi_mod + use mfpbl_mod contains !> \section arg_table_hedmf_init Argument Table diff --git a/physics/lsm_noah.f b/physics/lsm_noah.f index d519dcda5..7a8e17bf8 100644 --- a/physics/lsm_noah.f +++ b/physics/lsm_noah.f @@ -7,6 +7,7 @@ module lsm_noah use machine, only: kind_phys use set_soilveg_mod, only: set_soilveg use namelist_soilveg + use sflx implicit none diff --git a/physics/mfpbl.f b/physics/mfpbl.f index 2df84945b..dac548711 100644 --- a/physics/mfpbl.f +++ b/physics/mfpbl.f @@ -1,6 +1,7 @@ !> \file mfpbl.f !! This file contains the subroutine that calculates the updraft properties and mass flux for use in the Hybrid EDMF PBL scheme. - + module mfpbl_mod + contains !> \ingroup HEDMF !! \brief This subroutine is used for calculating the mass flux and updraft properties. !! @@ -396,3 +397,4 @@ subroutine mfpbl(im,ix,km,ntrac,delt,cnvflg, & return end !> @} + end module mfpbl_mod diff --git a/physics/mfpblt.f b/physics/mfpblt.f index bd0baf558..67e554b92 100644 --- a/physics/mfpblt.f +++ b/physics/mfpblt.f @@ -2,7 +2,8 @@ !! This file contains the subroutine that calculates mass flux and !! updraft parcel properties for thermals driven by surface heating !! for use in the TKE-EDMF PBL scheme. - + module mfpblt_mod + contains !>\ingroup satmedmf !! This subroutine computes mass flux and updraft parcel properties for !! thermals driven by surface heating. @@ -452,3 +453,4 @@ subroutine mfpblt(im,ix,km,kmpbl,ntcw,ntrac1,delt, & return end !> @} + end module mfpblt_mod diff --git a/physics/mfpbltq.f b/physics/mfpbltq.f index c4333290b..4555af268 100644 --- a/physics/mfpbltq.f +++ b/physics/mfpbltq.f @@ -2,7 +2,8 @@ !! This file contains the subroutine that calculates mass flux and !! updraft parcel properties for thermals driven by surface heating !! for use in the TKE-EDMF PBL scheme (updated version). - + module mfpbltq_mod + contains !>\ingroup satmedmfvdifq !! This subroutine computes mass flux and updraft parcel properties for !! thermals driven by surface heating. @@ -477,3 +478,4 @@ subroutine mfpbltq(im,ix,km,kmpbl,ntcw,ntrac1,delt, return end !> @} + end module mfpbltq_mod diff --git a/physics/mfscu.f b/physics/mfscu.f index 9128c7c10..e0c184139 100644 --- a/physics/mfscu.f +++ b/physics/mfscu.f @@ -1,7 +1,8 @@ !>\file mfscu.f !! This file contains the mass flux and downdraft parcel preperties !! parameterization for stratocumulus-top-driven turbulence. - + module mfscu_mod + contains !>\ingroup satmedmf !! This subroutine computes mass flux and downdraft parcel properties !! for stratocumulus-top-driven turbulence. @@ -554,3 +555,4 @@ subroutine mfscu(im,ix,km,kmscu,ntcw,ntrac1,delt, & return end !> @} + end module mfscu_mod diff --git a/physics/mfscuq.f b/physics/mfscuq.f index 3c54b0bda..ca4819956 100644 --- a/physics/mfscuq.f +++ b/physics/mfscuq.f @@ -1,7 +1,8 @@ !>\file mfscuq.f !! This file contains the mass flux and downdraft parcel preperties !! parameterization for stratocumulus-top-driven turbulence (updated version). - + module mfscuq_mod + contains !>\ingroup satmedmfvdifq !! This subroutine computes mass flux and downdraft parcel properties !! for stratocumulus-top-driven turbulence. @@ -557,3 +558,4 @@ subroutine mfscuq(im,ix,km,kmscu,ntcw,ntrac1,delt, return end !> @} + end module mfscuq_mod diff --git a/physics/module_bl_mynn.F90 b/physics/module_bl_mynn.F90 index 7b98b1c93..334d1db4c 100644 --- a/physics/module_bl_mynn.F90 +++ b/physics/module_bl_mynn.F90 @@ -1384,8 +1384,11 @@ SUBROUTINE boulac_length0(k,kts,kte,zw,dz,qtke,theta,lb1,lb2) dld = min(dld,zw(k+1))!not used in PBL anyway, only free atmos lb1 = min(dlu,dld) !minimum !JOE-fight floating point errors +#ifdef SINGLE_PREC + !JM: keep up the fight, JOE dlu=MAX(0.1,MIN(dlu,1000.)) dld=MAX(0.1,MIN(dld,1000.)) +#endif lb2 = sqrt(dlu*dld) !average - biased towards smallest !lb2 = 0.5*(dlu+dld) !average diff --git a/physics/moninshoc.f b/physics/moninshoc.f index 4e9e60b46..ee4715e81 100644 --- a/physics/moninshoc.f +++ b/physics/moninshoc.f @@ -4,6 +4,9 @@ !> This module contains the CCPP-compliant SHOC scheme. module moninshoc + use mfpbl_mod + use tridi_mod + contains subroutine moninshoc_init (do_shoc, errmsg, errflg) diff --git a/physics/satmedmfvdif.F b/physics/satmedmfvdif.F index feb4ef870..f791a2de4 100644 --- a/physics/satmedmfvdif.F +++ b/physics/satmedmfvdif.F @@ -5,6 +5,10 @@ module satmedmfvdif + use tridi_mod + use mfscu_mod + use mfpblt_mod + contains !> \section arg_table_satmedmfvdif_init Argument Table diff --git a/physics/satmedmfvdifq.F b/physics/satmedmfvdifq.F index eb2b7ad1c..9c5ad4029 100644 --- a/physics/satmedmfvdifq.F +++ b/physics/satmedmfvdifq.F @@ -4,7 +4,9 @@ !! eddy-diffusion mass-flux (TKE-EDMF) parameterization (by Jongil Han). module satmedmfvdifq - + use mfpbltq_mod + use tridi_mod + use mfscuq_mod contains !> \defgroup satmedmfvdifq GFS Scale-aware TKE-based Moist Eddy-Diffusivity Mass-flux (TKE-EDMF, updated version) Scheme Module diff --git a/physics/sflx.f b/physics/sflx.f index 61fe015cc..026e2b854 100644 --- a/physics/sflx.f +++ b/physics/sflx.f @@ -1,6 +1,7 @@ !>\file sflx.f !! This file is the entity of GFS Noah LSM Model(Version 2.7). - + module sflx + contains !>\ingroup Noah_LSM !!\brief This is the entity of GFS Noah LSM model of physics subroutines. !! It is a soil/veg/snowpack land-surface model to update soil moisture, soil @@ -906,7 +907,15 @@ subroutine gfssflx &! --- input eta = etp endif +#ifdef SINGLE_PREC + IF (ETP == 0.0) THEN + BETA = 0.0 + ELSE + BETA = ETA/ETP + ENDIF +#else beta = eta / etp +#endif !> - Convert the sign of soil heat flux so that: !! - ssoil>0: warm the surface (night time) @@ -5801,3 +5810,4 @@ end subroutine wdfcnd end subroutine gfssflx !! @} !----------------------------------- + end module sflx diff --git a/physics/tridi.f b/physics/tridi.f index 0103b388f..13202512f 100644 --- a/physics/tridi.f +++ b/physics/tridi.f @@ -1,6 +1,7 @@ !>\file tridi.f !! These subroutines are originally internal subroutines in moninedmf.f - + module tridi_mod + contains !>\ingroup HEDMF !!\brief Routine to solve the tridiagonal system to calculate !!temperature and moisture at \f$ t + \Delta t \f$; part of two-part @@ -220,3 +221,4 @@ subroutine tridit(l,n,nt,cl,cm,cu,rt,au,at) return end subroutine tridit !> @} + end module tridi_mod diff --git a/physics/ugwp_driver_v0.F b/physics/ugwp_driver_v0.F index 844acf722..cd19f5f71 100644 --- a/physics/ugwp_driver_v0.F +++ b/physics/ugwp_driver_v0.F @@ -1,5 +1,7 @@ !>\file ugwp_driver_v0.F - + module ugwp_driver_v0 + use cires_orowam2017 + contains ! !===================================================================== ! @@ -1485,3 +1487,4 @@ subroutine fv3_ugwp_solv2_v0(klon, klev, dtime, end subroutine fv3_ugwp_solv2_v0 + end module ugwp_driver_v0 diff --git a/physics/unified_ugwp.F90 b/physics/unified_ugwp.F90 index 9e93bd5fc..0b45d680d 100644 --- a/physics/unified_ugwp.F90 +++ b/physics/unified_ugwp.F90 @@ -37,7 +37,8 @@ module unified_ugwp ! use cires_ugwp_module, only: knob_ugwp_version, cires_ugwp_mod_init, cires_ugwp_mod_finalize use cires_ugwpv0_module, only: knob_ugwp_version, cires_ugwpv0_mod_init, cires_ugwpv0_mod_finalize use gwdps, only: gwdps_run - + use cires_ugwp_triggers + use ugwp_driver_v0 use drag_suite, only: drag_suite_run implicit none From be534e7be9c335c9d8736fe0d5bf8fe83e14cce9 Mon Sep 17 00:00:00 2001 From: Lisa Bengtsson Date: Thu, 5 May 2022 04:08:06 +0000 Subject: [PATCH 13/20] addressing some review comments --- physics/progsigma_calc.f90 | 114 +++++++++++++++++++------------------ physics/samfdeepcnv.f | 16 +++++- physics/samfdeepcnv.meta | 8 +++ 3 files changed, 80 insertions(+), 58 deletions(-) diff --git a/physics/progsigma_calc.f90 b/physics/progsigma_calc.f90 index c05af3003..fe74dc0c1 100644 --- a/physics/progsigma_calc.f90 +++ b/physics/progsigma_calc.f90 @@ -48,7 +48,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & real(kind=kind_phys) :: gcvalmx,epsilon,ZZ,cvg,mcon,buy2, & fdqb,dtdyn,dxlim,rmulacvg,tem, & - alpha,DEN,betascu,dp1 + alpha,DEN,betascu,dp1,invdelt !Parameters gcvalmx = 0.1 @@ -57,11 +57,11 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & km1=km-1 alpha=7000. betascu = 3.0 + invdelt = 1./delt !Initialization 2D do k = 1,km do i = 1,im - sigmaout(i,k)=0. inbu(i,k)=0. form(i,k)=0. dp(i,k)=0. @@ -70,7 +70,9 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & !Initialization 1D do i=1,im - sigmab(i)=0. + if(cnvflg(i))then + sigmab(i)=0. + endif sigmamax(i)=0.95 termA(i)=0. termB(i)=0. @@ -89,24 +91,16 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & enddo !Initial computations, place maximum sigmain in sigmab - if(flag_init .and. .not. flag_restart)then - do i=1,im - if(cnvflg(i))then - sigmab(i)=0.03 - endif - enddo - else - do i=1,im - if(cnvflg(i))then - do k=2,km - if(sigmain(i,k)>sigmab(i))then - sigmab(i)=sigmain(i,k) - endif - enddo - endif - enddo - endif - + do i=1,im + if(cnvflg(i))then + do k=2,km + if(sigmain(i,k)>sigmab(i))then + sigmab(i)=sigmain(i,k) + endif + enddo + endif + enddo + do i=1,im if(sigmab(i) < 1.E-5)then !after advection sigmab(i)=0. @@ -120,15 +114,20 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & enddo !Initial computations, dynamic q-tendency - do k = 1,km - do i = 1,im - if(flag_init .and. .not.flag_restart)then + if(flag_init .and. .not.flag_restart)then + do k = 1,km + do i = 1,im qadv(i,k)=0. - else - qadv(i,k)=(q(i,k) - prevsq(i,k))/delt - endif + enddo enddo - enddo + else + do k = 1,km + do i = 1,im + qadv(i,k)=(q(i,k) - prevsq(i,k))*invdelt + enddo + enddo + endif + !compute termD "The vertical integral of the latent heat convergence is limited to the !buoyant layers with positive moisture convergence (accumulated from the surface). @@ -173,7 +172,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & endif enddo enddo - + !termC do k = 2,km1 do i = 1,im @@ -188,33 +187,38 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & enddo !sigmab - do i = 1,im - if(cnvflg(i))then - DEN=MIN(termC(i)+termB(i),1.E8) - cvg=termD(i)*delt - ZZ=MAX(0.0,SIGN(1.0,termA(i))) & - *MAX(0.0,SIGN(1.0,termB(i))) & - *MAX(0.0,SIGN(1.0,termC(i)-epsilon)) - cvg=MAX(0.0,cvg) - if(flag_init .and. .not. flag_restart)then - sigmab(i)=0.03 - else - sigmab(i)=(ZZ*(termA(i)+cvg))/(DEN+(1.0-ZZ)) - endif - if(sigmab(i)>0.)then - sigmab(i)=MIN(sigmab(i),sigmamax(i)) - sigmab(i)=MAX(sigmab(i),0.01) - endif - endif!cnvflg - enddo + if(flag_init .and. .not. flag_restart)then + do i = 1,im + if(cnvflg(i))then + sigmab(i)=0.03 + endif + enddo + else + do i = 1,im + if(cnvflg(i))then + DEN=MIN(termC(i)+termB(i),1.E8) + cvg=termD(i)*delt + ZZ=MAX(0.0,SIGN(1.0,termA(i))) & + *MAX(0.0,SIGN(1.0,termB(i))) & + *MAX(0.0,SIGN(1.0,termC(i)-epsilon)) + cvg=MAX(0.0,cvg) + sigmab(i)=(ZZ*(termA(i)+cvg))/(DEN+(1.0-ZZ)) + if(sigmab(i)>0.)then + sigmab(i)=MIN(sigmab(i),sigmamax(i)) + sigmab(i)=MAX(sigmab(i),0.01) + endif + endif!cnvflg + enddo + endif - do k=1,km - do i=1,im - if(cnvflg(i))then - sigmaout(i,k)=sigmab(i) - endif - enddo - enddo + do k=1,km + do i=1,im + if(cnvflg(i))then + sigmaout(i,k)=sigmab(i) + endif + enddo + enddo + !Since updraft velocity is much lower in shallow cu region, termC becomes small in shallow cu application, thus the area fraction !in this regime becomes too large compared with the deep cu region. To address this simply apply a scaling factor for shallow cu diff --git a/physics/samfdeepcnv.f b/physics/samfdeepcnv.f index 071bf0557..8398af769 100644 --- a/physics/samfdeepcnv.f +++ b/physics/samfdeepcnv.f @@ -85,8 +85,8 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & & QLCN, QICN, w_upi, cf_upi, CNV_MFD, & & CNV_DQLDT,CLCN,CNV_FICE,CNV_NDROP,CNV_NICE,mp_phys,mp_phys_mg,& & clam,c0s,c1,betal,betas,evef,pgcon,asolfac, & - & do_ca, ca_closure, ca_entr, ca_trigger, nthresh, ca_deep, & - & rainevap, sigmain, sigmaout, errmsg,errflg) + & do_ca, ca_closure, ca_entr, ca_trigger, nthresh,ca_deep, & + & ca_micro, rainevap, sigmain, sigmaout, errmsg,errflg) ! use machine , only : kind_phys use funcphys , only : fpvs @@ -107,7 +107,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & real(kind=kind_phys), intent(in) :: ca_deep(:) real(kind=kind_phys), intent(in) :: sigmain(:,:),qmicro(:,:), & & tmf(:,:),q(:,:), prevsq(:,:) - real(kind=kind_phys), intent(out) :: rainevap(:) + real(kind=kind_phys), intent(out) :: rainevap(:), ca_micro(:) real(kind=kind_phys), intent(out) :: sigmaout(:,:) logical, intent(in) :: do_ca,ca_closure,ca_entr,ca_trigger @@ -2894,6 +2894,16 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & !> - 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 for the grid sizes where the quasi-equilibrium assumption of Arakawa-Schubert is not valid any longer. !! As discussed in Han et al. (2017) \cite han_et_al_2017 , when dtconv is larger than tauadv, the convective mixing is not fully conducted before the cumulus cloud is advected out of the grid cell. In this case, therefore, the cloud base mass flux is further reduced in proportion to the ratio of tauadv to dtconv. + do i=1,im + ca_micro(i)=0. + enddo + + do i=1,im + if(cnvflg(i))then + ca_micro(i)=sigmab(i) + endif + enddo + do i= 1, im if(cnvflg(i) .and. .not.asqecflg(i)) then k = kbcon(i) diff --git a/physics/samfdeepcnv.meta b/physics/samfdeepcnv.meta index 3f28035b6..1764a74fd 100644 --- a/physics/samfdeepcnv.meta +++ b/physics/samfdeepcnv.meta @@ -652,6 +652,14 @@ type = real kind = kind_phys intent = in +[ca_micro] + standard_name = output_prognostic_sigma_two + long_name = output of prognostic area fraction two + units = frac + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out [rainevap] standard_name = physics_field_for_coupling long_name = physics_field_for_coupling From 63020ec6a737511a46102865458b9843e340a404 Mon Sep 17 00:00:00 2001 From: "samuel.trahan" Date: Thu, 5 May 2022 22:46:10 +0000 Subject: [PATCH 14/20] Switch to another version of the code that works with 64 bit --- physics/GFS_rrtmgp_cloud_overlap.F90 | 4 ++-- physics/GFS_suite_interstitial_4.F90 | 10 +++++----- physics/cires_ugwp.F90 | 4 ---- physics/cires_ugwpv1_oro.F90 | 2 +- physics/hedmf.f | 2 ++ physics/maximum_hourly_diagnostics.F90 | 14 +++++++------- physics/module_bl_mynn.F90 | 23 ++++++++++++++++++----- physics/module_sf_mynn.F90 | 4 ++-- physics/module_sf_noahmplsm.f90 | 4 ++-- physics/module_sf_ruclsm.F90 | 5 ++--- physics/module_soil_pre.F90 | 24 +++++++++++------------- physics/radiation_gases.f | 2 +- physics/satmedmfvdif.F | 2 -- physics/surface_perturbation.F90 | 2 +- 14 files changed, 54 insertions(+), 48 deletions(-) diff --git a/physics/GFS_rrtmgp_cloud_overlap.F90 b/physics/GFS_rrtmgp_cloud_overlap.F90 index c1a6c4763..13794641b 100644 --- a/physics/GFS_rrtmgp_cloud_overlap.F90 +++ b/physics/GFS_rrtmgp_cloud_overlap.F90 @@ -99,7 +99,7 @@ subroutine GFS_rrtmgp_cloud_overlap_run(nCol, nLev, yearlen, doSWrad, doLWrad, ! Cloud overlap parameter ! if (iovr == iovr_dcorr .or. iovr == iovr_exp .or. iovr == iovr_exprand) then - call get_alpha_exper(nCol, nLev, iovr, iovr_exprand, deltaZc*0.001_kind_phys, de_lgth, cld_frac, cloud_overlap_param) + call get_alpha_exper(nCol, nLev, iovr, iovr_exprand, deltaZc*0.001, de_lgth, cld_frac, cloud_overlap_param) else de_lgth(:) = 0. cloud_overlap_param(:,:) = 0. @@ -110,7 +110,7 @@ subroutine GFS_rrtmgp_cloud_overlap_run(nCol, nLev, yearlen, doSWrad, doLWrad, ! if (imfdeepcnv == imfdeepcnv_samf .or. imfdeepcnv == imfdeepcnv_gf) then if (iovr_convcld == iovr_dcorr .or. iovr_convcld == iovr_exp .or. iovr_convcld == iovr_exprand) then - call get_alpha_exper(nCol, nLev, iovr_convcld, iovr_exprand, deltaZc*0.001_kind_phys, de_lgth, cld_cnv_frac, cnv_cloud_overlap_param) + call get_alpha_exper(nCol, nLev, iovr_convcld, iovr_exprand, deltaZc*0.001, de_lgth, cld_cnv_frac, cnv_cloud_overlap_param) else de_lgth(:) = 0. cnv_cloud_overlap_param(:,:) = 0. diff --git a/physics/GFS_suite_interstitial_4.F90 b/physics/GFS_suite_interstitial_4.F90 index 18fcfda09..cbabb991b 100644 --- a/physics/GFS_suite_interstitial_4.F90 +++ b/physics/GFS_suite_interstitial_4.F90 @@ -224,7 +224,7 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, tracers_total, ntr qc_mp(i,k) = (clw(i,k,2)-save_qc(i,k)) / (one-spechum(i,k)) !> - Convert number concentration from moist to dry nc_mp(i,k) = gq0(i,k,ntlnc) / (one-spechum(i,k)) - nc_mp(i,k) = max(zero, nc_mp(i,k) + make_DropletNumber(real(qc_mp(i,k) * rho), real(nwfa(i,k)*rho)) * orho) + nc_mp(i,k) = max(zero, nc_mp(i,k) + make_DropletNumber(qc_mp(i,k) * rho, nwfa(i,k)*rho) * orho) !> - Convert number concentrations from dry to moist gq0(i,k,ntlnc) = nc_mp(i,k) / (one+qv_mp(i,k)) endif @@ -233,7 +233,7 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, tracers_total, ntr qi_mp(i,k) = (clw(i,k,1)-save_qi(i,k)) / (one-spechum(i,k)) !> - Convert number concentration from moist to dry ni_mp(i,k) = gq0(i,k,ntinc) / (one-spechum(i,k)) - ni_mp(i,k) = max(zero, ni_mp(i,k) + make_IceNumber(real(qi_mp(i,k) * rho), real(save_tcp(i,k)) * orho)) + ni_mp(i,k) = max(zero, ni_mp(i,k) + make_IceNumber(qi_mp(i,k) * rho, save_tcp(i,k)) * orho) !> - Convert number concentrations from dry to moist gq0(i,k,ntinc) = ni_mp(i,k) / (one+qv_mp(i,k)) endif @@ -249,13 +249,13 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, tracers_total, ntr !> - Update cloud water mixing ratio qc_mp(i,k) = (clw(i,k,2)-save_qc(i,k)) !> - Update cloud water number concentration - gq0(i,k,ntlnc) = max(zero, gq0(i,k,ntlnc) + make_DropletNumber(real(qc_mp(i,k) * rho), real(nwfa(i,k)*rho)) * orho) + gq0(i,k,ntlnc) = max(zero, gq0(i,k,ntlnc) + make_DropletNumber(qc_mp(i,k) * rho, nwfa(i,k)*rho) * orho) endif if (ntinc>0) then !> - Update cloud ice mixing ratio qi_mp(i,k) = (clw(i,k,1)-save_qi(i,k)) !> - Update cloud ice number concentration - gq0(i,k,ntinc) = max(zero, gq0(i,k,ntinc) + make_IceNumber(real(qi_mp(i,k) * rho), real(save_tcp(i,k))) * orho) + gq0(i,k,ntinc) = max(zero, gq0(i,k,ntinc) + make_IceNumber(qi_mp(i,k) * rho, save_tcp(i,k)) * orho) endif enddo enddo @@ -290,4 +290,4 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, tracers_total, ntr end subroutine GFS_suite_interstitial_4_run - end module GFS_suite_interstitial_4 + end module GFS_suite_interstitial_4 \ No newline at end of file diff --git a/physics/cires_ugwp.F90 b/physics/cires_ugwp.F90 index 2d8eafc19..f2d6b3e3c 100644 --- a/physics/cires_ugwp.F90 +++ b/physics/cires_ugwp.F90 @@ -17,13 +17,9 @@ module cires_ugwp use cires_ugwpv0_module, only: knob_ugwp_version, cires_ugwpv0_mod_init, cires_ugwpv0_mod_finalize use ugwp_driver_v0 - use gwdps, only: gwdps_run - use cires_ugwp_triggers - use ugwp_driver_v0 - implicit none private diff --git a/physics/cires_ugwpv1_oro.F90 b/physics/cires_ugwpv1_oro.F90 index 66d0e472c..959bbd6c5 100644 --- a/physics/cires_ugwpv1_oro.F90 +++ b/physics/cires_ugwpv1_oro.F90 @@ -1,5 +1,5 @@ module cires_ugwpv1_oro - use cires_ugwpv1_sporo + use cires_ugwpv1_sporo contains subroutine orogw_v1 (im, km, imx, me, master, dtp, kdt, do_tofd, & diff --git a/physics/hedmf.f b/physics/hedmf.f index 604483e53..a1d8df9c3 100644 --- a/physics/hedmf.f +++ b/physics/hedmf.f @@ -5,8 +5,10 @@ !> This module contains the CCPP-compliant hybrid eddy-diffusivity mass-flux !! scheme. module hedmf + use tridi_mod use mfpbl_mod + contains !> \section arg_table_hedmf_init Argument Table diff --git a/physics/maximum_hourly_diagnostics.F90 b/physics/maximum_hourly_diagnostics.F90 index ddbff5725..6beae0da2 100644 --- a/physics/maximum_hourly_diagnostics.F90 +++ b/physics/maximum_hourly_diagnostics.F90 @@ -144,11 +144,11 @@ subroutine max_fields(phil,ref3D,grav,im,levs,refd,tk,refd263k) real (kind=kind_phys), intent(in) :: grav real (kind=kind_phys), intent(in),dimension(:,:) :: phil,ref3D,tk integer :: i,k,ll,ipt,kpt - real(kind_phys) :: dbz1avg,zmidp1,zmidloc,refl,fact - real(kind_phys), dimension(im,levs) :: z - real(kind_phys), dimension(im) :: zintsfc - real(kind_phys), dimension(:), intent(inout) :: refd,refd263k - REAL(kind_phys) :: dbz1(2),dbzk,dbzk1 + real :: dbz1avg,zmidp1,zmidloc,refl,fact + real, dimension(im,levs) :: z + real, dimension(im) :: zintsfc + real, dimension(:), intent(inout) :: refd,refd263k + REAL :: dbz1(2),dbzk,dbzk1 logical :: counter do i=1,im do k=1,levs @@ -185,7 +185,7 @@ subroutine max_fields(phil,ref3D,grav,im,levs,refd,tk,refd263k) dbz1avg=dbz1(2)+(dbz1(2)-dbz1(1))*fact !-- Convert to dBZ (10*logZ) as the last step if (dbz1avg>0.01) then - dbz1avg=10.*log10(dbz1avg) + dbz1avg=10.*alog10(dbz1avg) else dbz1avg=-35. endif @@ -214,7 +214,7 @@ subroutine max_fields(phil,ref3D,grav,im,levs,refd,tk,refd263k) dbz1avg=maxval(dbz1) !-- Convert to dBZ (10*logZ) as the last step if (dbz1avg>0.01) then - dbz1avg=10.*log10(dbz1avg) + dbz1avg=10.*alog10(dbz1avg) else dbz1avg=-35. endif diff --git a/physics/module_bl_mynn.F90 b/physics/module_bl_mynn.F90 index 334d1db4c..f16ca722a 100644 --- a/physics/module_bl_mynn.F90 +++ b/physics/module_bl_mynn.F90 @@ -1384,11 +1384,9 @@ SUBROUTINE boulac_length0(k,kts,kte,zw,dz,qtke,theta,lb1,lb2) dld = min(dld,zw(k+1))!not used in PBL anyway, only free atmos lb1 = min(dlu,dld) !minimum !JOE-fight floating point errors -#ifdef SINGLE_PREC !JM: keep up the fight, JOE dlu=MAX(0.1,MIN(dlu,1000.)) dld=MAX(0.1,MIN(dld,1000.)) -#endif lb2 = sqrt(dlu*dld) !average - biased towards smallest !lb2 = 0.5*(dlu+dld) !average @@ -1542,11 +1540,9 @@ SUBROUTINE boulac_length(kts,kte,zw,dz,qtke,theta,lb1,lb2) dld(iz) = min(dld(iz),zw(iz+1))!not used in PBL anyway, only free atmos lb1(iz) = min(dlu(iz),dld(iz)) !minimum !JOE-fight floating point errors -#ifdef SINGLE_PREC !JM: keep up the fight, JOE dlu(iz)=MAX(0.1,MIN(dlu(iz),1000.)) dld(iz)=MAX(0.1,MIN(dld(iz),1000.)) -#endif lb2(iz) = sqrt(dlu(iz)*dld(iz)) !average - biased towards smallest !lb2(iz) = 0.5*(dlu(iz)+dld(iz)) !average @@ -2955,8 +2951,12 @@ SUBROUTINE mym_condensation (kts,kte, & zagl = zagl + dz(k) !CLOUD WATER AND ICE - IF (q1k < 0.) THEN !unstaurated + IF (q1k < 0.) THEN !unsaturated +#ifdef SINGLE_PREC ql_water = sgm(k)*EXP(1.2*q1k-1.) +#else + ql_water = sgm(k)*EXP(1.2*q1k-1) +#endif ql_ice = sgm(k)*EXP(1.2*q1k-1.) !Reduce ice mixing ratios in the upper troposphere ! low_weight = MIN(MAX(p(k)-40000.0, 0.0),40000.0)/40000.0 @@ -7608,15 +7608,28 @@ FUNCTION qsat_blend(t, P, waterice) IF ((t .GE. 273.16) .OR. (wrt .EQ. 'w')) THEN ESL = J0+XC*(J1+XC*(J2+XC*(J3+XC*(J4+XC*(J5+XC*(J6+XC*(J7+XC*J8))))))) +#ifdef SINGLE_PREC qsat_blend = 0.622*ESL/max((P-ESL),1.0E-7_kind_phys) +#else + qsat_blend = 0.622*ESL/(P-ESL) +#endif ELSE IF (t .LE. 253.) THEN ESI = K0+XC*(K1+XC*(K2+XC*(K3+XC*(K4+XC*(K5+XC*(K6+XC*(K7+XC*K8))))))) +#ifdef SINGLE_PREC qsat_blend = 0.622*ESI/max((P-ESI),1.0E-7_kind_phys) +#else + qsat_blend = 0.622*ESI/(P-ESI) +#endif ELSE ESL = J0+XC*(J1+XC*(J2+XC*(J3+XC*(J4+XC*(J5+XC*(J6+XC*(J7+XC*J8))))))) ESI = K0+XC*(K1+XC*(K2+XC*(K3+XC*(K4+XC*(K5+XC*(K6+XC*(K7+XC*K8))))))) +#ifdef SINGLE_PREC RSLF = 0.622*ESL/max((P-ESL),1.0E-7_kind_phys) RSIF = 0.622*ESI/max((P-ESI),1.0E-7_kind_phys) +#else + RSLF = 0.622*ESL/(P-ESL) + RSIF = 0.622*ESI/(P-ESI) +#endif chi = (273.16-t)/20.16 qsat_blend = (1.-chi)*RSLF + chi*RSIF END IF diff --git a/physics/module_sf_mynn.F90 b/physics/module_sf_mynn.F90 index e14c23882..22b142c33 100644 --- a/physics/module_sf_mynn.F90 +++ b/physics/module_sf_mynn.F90 @@ -2804,8 +2804,8 @@ SUBROUTINE znot_m_v6(uref, znotm) ! znotm(meter): areodynamical roughness scale over water ! - REAL, INTENT(IN) :: uref - REAL, INTENT(OUT):: znotm + REAL(kind=kind_phys), INTENT(IN) :: uref + REAL(kind=kind_phys), INTENT(OUT):: znotm real(kind=kind_phys), parameter :: p13 = -1.296521881682694e-02,& & p12 = 2.855780863283819e-01, p11 = -1.597898515251717e+00,& & p10 = -8.396975715683501e+00, & diff --git a/physics/module_sf_noahmplsm.f90 b/physics/module_sf_noahmplsm.f90 index 61b92990b..1c899e4bd 100644 --- a/physics/module_sf_noahmplsm.f90 +++ b/physics/module_sf_noahmplsm.f90 @@ -681,7 +681,7 @@ subroutine noahmp_sflx (parameters, & logical :: dveg_active !< flag to run dynamic vegetation logical :: crop_active !< flag to run crop model ! add canopy heat storage (C.He added based on GY Niu's communication) - real (kind=kind_phys) :: canhs ! canopy heat storage change w/m2 + real :: canhs ! canopy heat storage change w/m2 ! maximum lai/sai used for some parameterizations based on plant growthi @@ -4494,7 +4494,7 @@ subroutine bare_flux (parameters,nsnow ,nsoil ,isnow ,dt ,sag , & if(opt_sfc == 3) then call sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur , & !in zlvl ,tgb ,thsfc_loc,prslkix,prsik1x ,prslk1x ,z0m , & !in - zpd ,snowh,shdfac ,garea1 ,.false. ,0.0_kind_phys,ivgtyp , & !in + zpd ,snowh,shdfac ,garea1 ,.false. ,0.0,ivgtyp , & !in ustarx ,fm ,fh ,fm2 ,fh2 , & !inout z0h ,fv ,csigmaf0,cm ,ch ) !out diff --git a/physics/module_sf_ruclsm.F90 b/physics/module_sf_ruclsm.F90 index a27d0f287..b39610bc8 100644 --- a/physics/module_sf_ruclsm.F90 +++ b/physics/module_sf_ruclsm.F90 @@ -7603,11 +7603,10 @@ END SUBROUTINE SOILIN !>\ingroup lsm_ruc_group !> This function calculates the liquid saturation vapor mixing ratio as !! a function of temperature and pressure (from Thompson scheme). - FUNCTION RSLF(P,T) + REAL FUNCTION RSLF(P,T) IMPLICIT NONE - REAL(kind_phys), INTENT(IN):: P, T - REAL(kind_phys) :: RSLF + REAL, INTENT(IN):: P, T REAL:: ESL,X REAL, PARAMETER:: C0= .611583699E03 REAL, PARAMETER:: C1= .444606896E02 diff --git a/physics/module_soil_pre.F90 b/physics/module_soil_pre.F90 index 149f87a1c..8eb5a5775 100644 --- a/physics/module_soil_pre.F90 +++ b/physics/module_soil_pre.F90 @@ -5,8 +5,6 @@ module module_soil_pre !tgs Initialize RUC LSM levels, soil temp/moisture - use machine, only: kind_phys - implicit none private @@ -28,8 +26,8 @@ SUBROUTINE init_soil_depth_3 ( zs , dzs , num_soil_levels ) INTEGER, INTENT(IN) :: num_soil_levels - REAL(kind_phys), DIMENSION(1:num_soil_levels), INTENT(OUT) :: zs, dzs - REAL(kind_phys), DIMENSION(1:num_soil_levels) :: zs2 + REAL, DIMENSION(1:num_soil_levels), INTENT(OUT) :: zs, dzs + REAL, DIMENSION(1:num_soil_levels) :: zs2 INTEGER :: l @@ -92,21 +90,21 @@ SUBROUTINE init_soil_3_real ( tsk , tmn , smois , tslb , & INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input - REAL(kind_phys) , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input - REAL(kind_phys) , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input - REAL(kind_phys) , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst + REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input + REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input + REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst - REAL(kind_phys) , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn - REAL(kind_phys) , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk - REAL(kind_phys) , DIMENSION(num_soil_layers) :: zs , dzs + REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn + REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk + REAL , DIMENSION(num_soil_layers) :: zs , dzs - REAL(kind_phys) , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois + REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois - REAL(kind_phys) , ALLOCATABLE , DIMENSION(:) :: zhave + REAL , ALLOCATABLE , DIMENSION(:) :: zhave logical :: debug_print = .false. INTEGER :: i , j , l , lout , lin , lwant , lhave, k - REAL(kind_phys) :: temp + REAL :: temp ! Allocate the soil layer array used for interpolating. diff --git a/physics/radiation_gases.f b/physics/radiation_gases.f index d6f1d7259..157da8e09 100644 --- a/physics/radiation_gases.f +++ b/physics/radiation_gases.f @@ -371,7 +371,7 @@ subroutine gas_init & endif do k = 1, LOZ - pkstr(k) = fpkapx(pstr(k)*100.0_kind_phys) + pkstr(k) = fpkapx(pstr(k)*100.0) enddo endif ! end if_ioznflg_block diff --git a/physics/satmedmfvdif.F b/physics/satmedmfvdif.F index f791a2de4..c7fe1d5c0 100644 --- a/physics/satmedmfvdif.F +++ b/physics/satmedmfvdif.F @@ -4,11 +4,9 @@ !! eddy-diffusion mass-flux (TKE-EDMF) parameterization (by Jongil Han). module satmedmfvdif - use tridi_mod use mfscu_mod use mfpblt_mod - contains !> \section arg_table_satmedmfvdif_init Argument Table diff --git a/physics/surface_perturbation.F90 b/physics/surface_perturbation.F90 index 7ddbe5279..e0429a5fc 100644 --- a/physics/surface_perturbation.F90 +++ b/physics/surface_perturbation.F90 @@ -48,7 +48,7 @@ subroutine cdfnor(z,cdfz) cdfz = 0.5 else x = 0.5*z*z - call cdfgam(x,0.5_kind_phys,del,iflag, cdfx) + call cdfgam(x,0.5,del,iflag, cdfx) if (iflag.ne.0) return if (z.gt.0.0) then cdfz = 0.5+0.5*cdfx From 49c70967f78e3a5d5293bab43c81acc6076cd060 Mon Sep 17 00:00:00 2001 From: Grant Firl Date: Tue, 10 May 2022 13:49:55 -0400 Subject: [PATCH 15/20] make sure that tsfc_wat is calculated when wet = T --- physics/scm_sfc_flux_spec.F90 | 21 +++++++++++++-------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/physics/scm_sfc_flux_spec.F90 b/physics/scm_sfc_flux_spec.F90 index fc4aaf5d1..bb2c47f48 100644 --- a/physics/scm_sfc_flux_spec.F90 +++ b/physics/scm_sfc_flux_spec.F90 @@ -153,14 +153,11 @@ subroutine scm_sfc_flux_spec_run (im, u1, v1, z1, t1, q1, p1, roughness_length, frland(i) = 1.0_kind_phys cice(i) = 0.0_kind_phys icy(i) = .false. - tsfcl(i) = T_surf(i) !GJF else frland(i) = 0.0_kind_phys if (oceanfrac(i) > 0.0_kind_phys) then if (cice(i) >= min_seaice) then icy(i) = .true. - tisfc(i) = T_surf(i) !GJF - tisfc(i) = max(timin, min(tisfc(i), tgice)) ! This cplice namelist option was added to deal with the ! situation of the FV3ATM-HYCOM coupling without an active sea ! ice (e.g., CICE6) component. By default, the cplice is true @@ -186,8 +183,6 @@ subroutine scm_sfc_flux_spec_run (im, u1, v1, z1, t1, q1, p1, roughness_length, else if (cice(i) >= min_lakeice) then icy(i) = .true. - tisfc(i) = T_surf(i) !GJF - tisfc(i) = max(timin, min(tisfc(i), tgice)) islmsk(i) = 2 else cice(i) = 0.0_kind_phys @@ -198,13 +193,23 @@ subroutine scm_sfc_flux_spec_run (im, u1, v1, z1, t1, q1, p1, roughness_length, if (cice(i) < 1.0_kind_phys) then wet(i) = .true. ! some open lake endif - if (wet(i)) then ! Water - tsfc_wat(i) = T_surf(i) - endif endif endif if (nint(slmsk(i)) /= 1) slmsk(i) = islmsk(i) enddo + + do i = 1, im + if (wet(i)) then + tsfc_wat(i) = T_surf(i) + end if + if (dry(i)) then + tsfcl(i) = T_surf(i) + end if + if (icy(i)) then + tisfc(i) = T_surf(i) + tisfc(i) = max(timin, min(tisfc(i), tgice)) + end if + end do ! to prepare to separate lake from ocean under water category do i = 1, im From b9940633e43fb6b14cfdb908c398b08fb6f1291a Mon Sep 17 00:00:00 2001 From: Dustin Swales Date: Mon, 16 May 2022 09:16:30 -0600 Subject: [PATCH 16/20] Update rte-rrtmgp submodule --- physics/rte-rrtmgp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/physics/rte-rrtmgp b/physics/rte-rrtmgp index cec1e8e12..f9377e81d 160000 --- a/physics/rte-rrtmgp +++ b/physics/rte-rrtmgp @@ -1 +1 @@ -Subproject commit cec1e8e12d969c3c8c76574dbe4f40b366419cc7 +Subproject commit f9377e81d33e4f73f4433501186465b84dd1111c From 6f38cc6e83526f907ec2c34712eb230d677f4a2a Mon Sep 17 00:00:00 2001 From: Lisa Bengtsson Date: Wed, 18 May 2022 23:58:09 +0000 Subject: [PATCH 17/20] address some review comments, fix decomposition error, correct bug in initialization --- physics/GFS_suite_interstitial_3.F90 | 37 +++++++++-- physics/GFS_suite_interstitial_3.meta | 73 +++++++++++++++++++++ physics/progsigma_calc.f90 | 93 ++++++++++++--------------- physics/samfdeepcnv.f | 22 ++----- physics/samfdeepcnv.meta | 8 --- physics/samfshalcnv.f | 20 +++--- 6 files changed, 162 insertions(+), 91 deletions(-) diff --git a/physics/GFS_suite_interstitial_3.F90 b/physics/GFS_suite_interstitial_3.F90 index 79ab481ec..9fa7d69b7 100644 --- a/physics/GFS_suite_interstitial_3.F90 +++ b/physics/GFS_suite_interstitial_3.F90 @@ -9,10 +9,13 @@ module GFS_suite_interstitial_3 !! \htmlinclude GFS_suite_interstitial_3_run.html !! subroutine GFS_suite_interstitial_3_run (otsptflag, & - im, levs, nn, cscnv, & + im, levs, nn, cscnv,imfshalcnv, imfdeepcnv, & + imfshalcnv_samf, imfdeepcnv_samf,progsigma, & + first_time_step, restart, & satmedmf, trans_trac, do_shoc, ltaerosol, ntrac, ntcw, & ntiw, ntclamt, ntrw, ntsw, ntrnc, ntsnc, ntgl, ntgnc, & - xlon, xlat, gt0, gq0, imp_physics, imp_physics_mg, & + xlon, xlat, gt0, gq0, sigmain,sigmaout,qmicro, & + imp_physics, imp_physics_mg, & imp_physics_zhao_carr, imp_physics_zhao_carr_pdf, & imp_physics_gfdl, imp_physics_thompson, dtidx, ntlnc, & imp_physics_wsm6, imp_physics_fer_hires, prsi, ntinc, & @@ -33,8 +36,9 @@ subroutine GFS_suite_interstitial_3_run (otsptflag, & imp_physics_gfdl, imp_physics_thompson, imp_physics_wsm6,imp_physics_fer_hires, & imp_physics_nssl, me, index_of_process_conv_trans integer, intent(in ), dimension(:) :: islmsk, kpbl, kinver - logical, intent(in ) :: cscnv, satmedmf, trans_trac, do_shoc, ltaerosol, ras - + logical, intent(in ) :: cscnv, satmedmf, trans_trac, do_shoc, ltaerosol, ras, progsigma + logical, intent(in ) :: first_time_step, restart + integer, intent(in ) :: imfshalcnv, imfdeepcnv, imfshalcnv_samf,imfdeepcnv_samf integer, intent(in) :: ntinc, ntlnc logical, intent(in) :: ldiag3d, qdiag3d integer, dimension(:,:), intent(in) :: dtidx @@ -48,6 +52,8 @@ subroutine GFS_suite_interstitial_3_run (otsptflag, & real(kind=kind_phys), intent(in ), dimension(:,:) :: gt0 real(kind=kind_phys), intent(in ), dimension(:,:,:) :: gq0 + real(kind=kind_phys), intent(out ), dimension(:,:) :: sigmain + real(kind=kind_phys), intent(out ), dimension(:,:) :: sigmaout,qmicro real(kind=kind_phys), intent(inout), dimension(:,:) :: rhc, save_qc ! save_qi is not allocated for Zhao-Carr MP real(kind=kind_phys), intent(inout), dimension(:,:) :: save_qi @@ -73,6 +79,27 @@ subroutine GFS_suite_interstitial_3_run (otsptflag, & errmsg = '' errflg = 0 + ! In case of using prognostic updraf area fraction, initialize area fraction here + ! since progsigma_calc is called from both deep and shallow schemes. + if(((imfshalcnv == imfshalcnv_samf) .or. (imfdeepcnv == imfdeepcnv_samf)) & + .and. progsigma)then + if(first_time_step .and. .not. restart)then + do k=1,levs + do i=1,im + sigmain(i,k)=0.0 + sigmaout(i,k)=0.0 + qmicro(i,k)=0.0 + enddo + enddo + endif + do k=1,levs + do i=1,im + sigmaout(i,k)=0.0 + enddo + enddo + endif + + if (cscnv .or. satmedmf .or. trans_trac .or. ras) then tracers = 2 do n=2,ntrac @@ -192,4 +219,4 @@ subroutine GFS_suite_interstitial_3_run (otsptflag, & end subroutine GFS_suite_interstitial_3_run - end module GFS_suite_interstitial_3 \ No newline at end of file + end module GFS_suite_interstitial_3 diff --git a/physics/GFS_suite_interstitial_3.meta b/physics/GFS_suite_interstitial_3.meta index 22a11d0ea..fbeb9f03c 100644 --- a/physics/GFS_suite_interstitial_3.meta +++ b/physics/GFS_suite_interstitial_3.meta @@ -43,6 +43,55 @@ dimensions = () type = logical intent = in +[imfdeepcnv] + standard_name = control_for_deep_convection_scheme + long_name = flag for mass-flux deep convection scheme + units = flag + dimensions = () + type = integer + intent = in +[imfdeepcnv_samf] + standard_name = identifer_for_scale_aware_mass_flux_deep_convection + long_name = flag for SAMF deep convection scheme + units = flag + dimensions = () + type = integer + intent = in +[imfshalcnv] + standard_name = control_for_shallow_convection_scheme + long_name = flag for mass-flux shallow convection scheme + units = flag + dimensions = () + type = integer + intent = in +[imfshalcnv_samf] + standard_name = identifier_for_scale_aware_mass_flux_shallow_convection + long_name = flag for SAMF shallow convection scheme + units = flag + dimensions = () + type = integer + intent = in +[progsigma] + standard_name = do_prognostic_updraft_area_fraction + long_name = flag for prognostic sigma in cumuls scheme + units = flag + dimensions = () + type = logical + intent = in +[first_time_step] + standard_name = flag_for_first_timestep + long_name = flag for first time step for time integration loop (cold/warmstart) + units = flag + dimensions = () + type = logical + intent = in +[restart] + standard_name = flag_for_restart + long_name = flag for restart (warmstart) or coldstart + units = flag + dimensions = () + type = logical + intent = in [satmedmf] standard_name = flag_for_scale_aware_TKE_moist_EDMF_PBL long_name = flag for scale-aware TKE moist EDMF PBL scheme @@ -173,6 +222,30 @@ type = real kind = kind_phys intent = in +[sigmain] + standard_name = prognostic_updraft_area_fraction_in_convection + long_name = convective updraft area fraction + units = frac + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = out +[sigmaout] + standard_name = updraft_area_fraction_updated_by_physics + long_name = convective updraft area fraction updated by physics + units = frac + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = out +[qmicro] + standard_name = instantaneous_tendency_of_specific_humidity_due_to_microphysics + long_name = moisture tendency due to microphysics + units = kg kg-1 s-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = out [imp_physics] standard_name = control_for_microphysics_scheme long_name = choice of microphysics scheme diff --git a/physics/progsigma_calc.f90 b/physics/progsigma_calc.f90 index fe74dc0c1..49a5e2a4f 100644 --- a/physics/progsigma_calc.f90 +++ b/physics/progsigma_calc.f90 @@ -14,9 +14,9 @@ !> @{ subroutine progsigma_calc (im,km,flag_init,flag_restart, & - flag_shallow,del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap, & - delt,prevsq,q,kbcon1,ktcon,cnvflg,gdx, & - sigmain,sigmaout,sigmab,errmsg,errflg) + del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap, & + delt,prevsq,q,kbcon1,ktcon,cnvflg,sigmain,sigmaout, & + sigmab,errmsg,errflg) ! ! use machine, only : kind_phys @@ -29,8 +29,8 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & real(kind=kind_phys), intent(in) :: hvap,delt real(kind=kind_phys), intent(in) :: prevsq(im,km), q(im,km),del(im,km), & qmicro(im,km),tmf(im,km),dbyo1(im,km),zdqca(im,km), & - omega_u(im,km),zeta(im,km),gdx(im) - logical, intent(in) :: flag_init,flag_restart,cnvflg(im),flag_shallow + omega_u(im,km),zeta(im,km) + logical, intent(in) :: flag_init,flag_restart,cnvflg(im) real(kind=kind_phys), intent(in) :: sigmain(im,km) ! intent out @@ -41,21 +41,20 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & ! Local variables integer :: i,k,km1 - real(kind=kind_phys) :: termA(im),termB(im),termC(im),termD(im), & - mcons(im),fdqa(im),form(im,km), & - qadv(im,km),sigmamax(im),dp(im,km),inbu(im,km) + real(kind=kind_phys) :: termA(im),termB(im),termC(im),termD(im) + real(kind=kind_phys) :: mcons(im),fdqa(im),form(im,km), & + qadv(im,km),dp(im,km),inbu(im,km) real(kind=kind_phys) :: gcvalmx,epsilon,ZZ,cvg,mcon,buy2, & fdqb,dtdyn,dxlim,rmulacvg,tem, & - alpha,DEN,betascu,dp1,invdelt + DEN,betascu,dp1,invdelt !Parameters gcvalmx = 0.1 rmulacvg=10. epsilon=1.E-11 km1=km-1 - alpha=7000. betascu = 3.0 invdelt = 1./delt @@ -70,10 +69,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & !Initialization 1D do i=1,im - if(cnvflg(i))then - sigmab(i)=0. - endif - sigmamax(i)=0.95 + sigmab(i)=0. termA(i)=0. termB(i)=0. termC(i)=0. @@ -82,6 +78,21 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & mcons(i)=0. enddo + !Initial computations, dynamic q-tendency + if(flag_init .and. .not.flag_restart)then + do k = 1,km + do i = 1,im + qadv(i,k)=0. + enddo + enddo + else + do k = 1,km + do i = 1,im + qadv(i,k)=(q(i,k) - prevsq(i,k))*invdelt + enddo + enddo + endif + do k = 2,km1 do i = 1,im if(cnvflg(i))then @@ -102,33 +113,13 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & enddo do i=1,im - if(sigmab(i) < 1.E-5)then !after advection - sigmab(i)=0. + if(cnvflg(i))then + if(sigmab(i) < 1.E-5)then !after advection + sigmab(i)=0. + endif endif enddo - - !Initial computations, sigmamax - do i=1,im - sigmamax(i)=alpha/gdx(i) - sigmamax(i)=MIN(0.95,sigmamax(i)) - enddo - - !Initial computations, dynamic q-tendency - if(flag_init .and. .not.flag_restart)then - do k = 1,km - do i = 1,im - qadv(i,k)=0. - enddo - enddo - else - do k = 1,km - do i = 1,im - qadv(i,k)=(q(i,k) - prevsq(i,k))*invdelt - enddo - enddo - endif - - + !compute termD "The vertical integral of the latent heat convergence is limited to the !buoyant layers with positive moisture convergence (accumulated from the surface). !Lowest level: @@ -194,7 +185,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & endif enddo else - do i = 1,im + do i = 1,im if(cnvflg(i))then DEN=MIN(termC(i)+termB(i),1.E8) cvg=termD(i)*delt @@ -204,7 +195,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & cvg=MAX(0.0,cvg) sigmab(i)=(ZZ*(termA(i)+cvg))/(DEN+(1.0-ZZ)) if(sigmab(i)>0.)then - sigmab(i)=MIN(sigmab(i),sigmamax(i)) + sigmab(i)=MIN(sigmab(i),0.95) sigmab(i)=MAX(sigmab(i),0.01) endif endif!cnvflg @@ -218,19 +209,15 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & endif enddo enddo - - !Since updraft velocity is much lower in shallow cu region, termC becomes small in shallow cu application, thus the area fraction - !in this regime becomes too large compared with the deep cu region. To address this simply apply a scaling factor for shallow cu - !before computing the massflux to reduce the total strength of the SC MF: - - if(flag_shallow)then - do i= 1, im - if(cnvflg(i)) then - sigmab(i)=sigmab(i)/betascu - endif - enddo - endif + !Reduce area fraction before coupling back to mass-flux computation. + !This tuning could be addressed in updraft velocity equation instead. + do i= 1, im + if(cnvflg(i)) then + sigmab(i)=sigmab(i)/betascu + endif + enddo + end subroutine progsigma_calc diff --git a/physics/samfdeepcnv.f b/physics/samfdeepcnv.f index 8398af769..3968061a2 100644 --- a/physics/samfdeepcnv.f +++ b/physics/samfdeepcnv.f @@ -86,7 +86,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & & CNV_DQLDT,CLCN,CNV_FICE,CNV_NDROP,CNV_NICE,mp_phys,mp_phys_mg,& & clam,c0s,c1,betal,betas,evef,pgcon,asolfac, & & do_ca, ca_closure, ca_entr, ca_trigger, nthresh,ca_deep, & - & ca_micro, rainevap, sigmain, sigmaout, errmsg,errflg) + & rainevap,sigmain, sigmaout, errmsg,errflg) ! use machine , only : kind_phys use funcphys , only : fpvs @@ -107,7 +107,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & real(kind=kind_phys), intent(in) :: ca_deep(:) real(kind=kind_phys), intent(in) :: sigmain(:,:),qmicro(:,:), & & tmf(:,:),q(:,:), prevsq(:,:) - real(kind=kind_phys), intent(out) :: rainevap(:), ca_micro(:) + real(kind=kind_phys), intent(out) :: rainevap(:) real(kind=kind_phys), intent(out) :: sigmaout(:,:) logical, intent(in) :: do_ca,ca_closure,ca_entr,ca_trigger @@ -214,8 +214,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & ! ! parameters for prognostic sigma closure real(kind=kind_phys) omega_u(im,km),zdqca(im,km),qlks(im,km), - & omegac(im),zeta(im,km),dbyo1(im,km),sigmab(im) - logical flag_shallow + & omegac(im),zeta(im,km),dbyo1(im,km),sigmab(im) real(kind=kind_phys) gravinv c physical parameters ! parameter(grav=grav,asolfac=0.958) @@ -2884,25 +2883,14 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & !> - From Bengtsson et al. (2022) Prognostic closure scheme, equation 8, compute updraft area fraction based on a moisture budget if(progsigma)then - flag_shallow = .false. - call progsigma_calc(im,km,first_time_step,restart,flag_shallow, + call progsigma_calc(im,km,first_time_step,restart, & del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt, - & prevsq,q,kbcon1,ktcon,cnvflg,gdx, + & prevsq,q,kbcon1,ktcon,cnvflg, & sigmain,sigmaout,sigmab,errmsg,errflg) 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 for the grid sizes where the quasi-equilibrium assumption of Arakawa-Schubert is not valid any longer. !! As discussed in Han et al. (2017) \cite han_et_al_2017 , when dtconv is larger than tauadv, the convective mixing is not fully conducted before the cumulus cloud is advected out of the grid cell. In this case, therefore, the cloud base mass flux is further reduced in proportion to the ratio of tauadv to dtconv. - - do i=1,im - ca_micro(i)=0. - enddo - - do i=1,im - if(cnvflg(i))then - ca_micro(i)=sigmab(i) - endif - enddo do i= 1, im if(cnvflg(i) .and. .not.asqecflg(i)) then diff --git a/physics/samfdeepcnv.meta b/physics/samfdeepcnv.meta index 1764a74fd..3f28035b6 100644 --- a/physics/samfdeepcnv.meta +++ b/physics/samfdeepcnv.meta @@ -652,14 +652,6 @@ type = real kind = kind_phys intent = in -[ca_micro] - standard_name = output_prognostic_sigma_two - long_name = output of prognostic area fraction two - units = frac - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = out [rainevap] standard_name = physics_field_for_coupling long_name = physics_field_for_coupling diff --git a/physics/samfshalcnv.f b/physics/samfshalcnv.f index ef0366b84..b7943725b 100644 --- a/physics/samfshalcnv.f +++ b/physics/samfshalcnv.f @@ -62,7 +62,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & prslp,psp,phil,qtr,prevsq,q,q1,t1,u1,v1,fscav, & & rn,kbot,ktop,kcnv,islimsk,garea, & & dot,ncloud,hpbl,ud_mf,dt_mf,cnvw,cnvc, & - & clam,c0s,c1,evef,pgcon,asolfac,hwrf_samfshal, + & clam,c0s,c1,evef,pgcon,asolfac,hwrf_samfshal, & & sigmain,sigmaout,errmsg,errflg) ! use machine , only : kind_phys @@ -77,7 +77,9 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & real(kind=kind_phys), intent(in) :: delt real(kind=kind_phys), intent(in) :: psp(:), delp(:,:), & & prslp(:,:), garea(:), hpbl(:), dot(:,:), phil(:,:), & - & qmicro(:,:),tmf(:,:),prevsq(:,:),q(:,:),sigmain(:,:) + & qmicro(:,:),tmf(:,:),prevsq(:,:),q(:,:) + + real(kind=kind_phys), intent(in) :: sigmain(:,:) ! real(kind=kind_phys), dimension(:), intent(in) :: fscav integer, intent(inout) :: kcnv(:) @@ -164,8 +166,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & real(kind=kind_phys) omega_u(im,km),zdqca(im,km),qlks(im,km), & omegac(im),zeta(im,km),dbyo1(im,km), & sigmab(im) - logical flag_shallow - real(kind=kind_phys) gravinv + real(kind=kind_phys) gravinv,dxcrtas c physical parameters ! parameter(g=grav,asolfac=0.89) @@ -196,6 +197,8 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & parameter(bet1=1.875,cd1=.506,f1=2.0,gam1=.5) parameter(betaw=.03,dxcrt=15.e3,dxcrtc0=9.e3) parameter(h1=0.33333333) +! progsigma + parameter(dxcrtas=30.e3) 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), @@ -340,6 +343,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & enddo endif !! + !> - Return to the calling routine if deep convection is present or the surface buoyancy flux is negative. totflg = .true. do i=1,im @@ -1932,20 +1936,20 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & c c Prognostic closure if(progsigma)then - flag_shallow = .true. - call progsigma_calc(im,km,first_time_step,restart,flag_shallow, + call progsigma_calc(im,km,first_time_step,restart, & del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt, - & prevsq,q,kbcon1,ktcon,cnvflg,gdx, + & prevsq,q,kbcon1,ktcon,cnvflg, & sigmain,sigmaout,sigmab,errmsg,errflg) 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. !! As discussed in Han et al. (2017) \cite han_et_al_2017 , when dtconv is larger than tauadv, the convective mixing is not fully conducted before the cumulus cloud is advected out of the grid cell. In this case, therefore, the cloud base mass flux is further reduced in proportion to the ratio of tauadv to dtconv. + do i= 1, im if(cnvflg(i)) then k = kbcon(i) rho = po(i,k)*100. / (rd*to(i,k)) - if(progsigma)then + if(progsigma .and. gdx(i) < dxcrtas)then xmb(i) = advfac(i)*sigmab(i)*((-1.0*omegac(i))*gravinv) else xmb(i) = advfac(i)*betaw*rho*wc(i) From fc79cc39ecaef4bc567b052842dc560e3253f0d8 Mon Sep 17 00:00:00 2001 From: Lisa Bengtsson Date: Thu, 19 May 2022 19:59:43 +0000 Subject: [PATCH 18/20] Change intent to inout for conditional variables --- physics/GFS_MP_generic_post.F90 | 4 ++-- physics/GFS_MP_generic_post.meta | 4 ++-- physics/GFS_suite_interstitial_3.F90 | 4 ++-- physics/GFS_suite_interstitial_3.meta | 4 ++-- physics/satmedmfvdifq.F | 4 ++-- physics/satmedmfvdifq.meta | 2 +- 6 files changed, 11 insertions(+), 11 deletions(-) diff --git a/physics/GFS_MP_generic_post.F90 b/physics/GFS_MP_generic_post.F90 index 3cdfe91f3..992c9b969 100644 --- a/physics/GFS_MP_generic_post.F90 +++ b/physics/GFS_MP_generic_post.F90 @@ -81,8 +81,8 @@ subroutine GFS_MP_generic_post_run( real(kind=kind_phys), dimension(:), intent(inout) :: diceprv real(kind=kind_phys), dimension(:), intent(inout) :: dsnowprv real(kind=kind_phys), dimension(:), intent(inout) :: dgraupelprv - real(kind=kind_phys), dimension(:,:), intent(out) :: dqdt_qmicro - real(kind=kind_phys), dimension(:,:), intent(out) :: prevsq + real(kind=kind_phys), dimension(:,:), intent(inout) :: dqdt_qmicro + real(kind=kind_phys), dimension(:,:), intent(inout) :: prevsq real(kind=kind_phys), intent(in) :: dtp ! CCPP error handling diff --git a/physics/GFS_MP_generic_post.meta b/physics/GFS_MP_generic_post.meta index 20881fbb3..7ba09363a 100644 --- a/physics/GFS_MP_generic_post.meta +++ b/physics/GFS_MP_generic_post.meta @@ -738,7 +738,7 @@ dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys - intent = out + intent = inout [prevsq] standard_name = specific_humidity_on_previous_timestep long_name = specific_humidity_on_previous_timestep @@ -746,7 +746,7 @@ dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys - intent = out + intent = inout [lssav] standard_name = flag_for_diagnostics long_name = logical flag for storing diagnostics diff --git a/physics/GFS_suite_interstitial_3.F90 b/physics/GFS_suite_interstitial_3.F90 index 9fa7d69b7..0f666d4a5 100644 --- a/physics/GFS_suite_interstitial_3.F90 +++ b/physics/GFS_suite_interstitial_3.F90 @@ -52,8 +52,8 @@ subroutine GFS_suite_interstitial_3_run (otsptflag, & real(kind=kind_phys), intent(in ), dimension(:,:) :: gt0 real(kind=kind_phys), intent(in ), dimension(:,:,:) :: gq0 - real(kind=kind_phys), intent(out ), dimension(:,:) :: sigmain - real(kind=kind_phys), intent(out ), dimension(:,:) :: sigmaout,qmicro + real(kind=kind_phys), intent(inout ), dimension(:,:) :: sigmain + real(kind=kind_phys), intent(inout ), dimension(:,:) :: sigmaout,qmicro real(kind=kind_phys), intent(inout), dimension(:,:) :: rhc, save_qc ! save_qi is not allocated for Zhao-Carr MP real(kind=kind_phys), intent(inout), dimension(:,:) :: save_qi diff --git a/physics/GFS_suite_interstitial_3.meta b/physics/GFS_suite_interstitial_3.meta index fbeb9f03c..fe52a1adc 100644 --- a/physics/GFS_suite_interstitial_3.meta +++ b/physics/GFS_suite_interstitial_3.meta @@ -229,7 +229,7 @@ dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys - intent = out + intent = inout [sigmaout] standard_name = updraft_area_fraction_updated_by_physics long_name = convective updraft area fraction updated by physics @@ -237,7 +237,7 @@ dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys - intent = out + intent = inout [qmicro] standard_name = instantaneous_tendency_of_specific_humidity_due_to_microphysics long_name = moisture tendency due to microphysics diff --git a/physics/satmedmfvdifq.F b/physics/satmedmfvdifq.F index c7a6fadc9..865e4481c 100644 --- a/physics/satmedmfvdifq.F +++ b/physics/satmedmfvdifq.F @@ -99,7 +99,7 @@ subroutine satmedmfvdifq_run(im,km,progsigma,ntrac,ntcw,ntrw, & real(kind=kind_phys), intent(in) :: dspfac, bl_upfr, bl_dnfr real(kind=kind_phys), intent(in) :: rlmx, elmx real(kind=kind_phys), intent(inout) :: dv(:,:), du(:,:), & - & tdt(:,:), rtg(:,:,:) + & tdt(:,:), rtg(:,:,:), tmf(:,:) real(kind=kind_phys), intent(in) :: & & u1(:,:), v1(:,:), & & t1(:,:), q1(:,:,:), & @@ -123,7 +123,7 @@ subroutine satmedmfvdifq_run(im,km,progsigma,ntrac,ntcw,ntrw, & & dtsfc(:), dqsfc(:), & & hpbl(:) real(kind=kind_phys), intent(out) :: & - & dkt(:,:), dku(:,:), tmf(:,:) + & dkt(:,:), dku(:,:) ! logical, intent(in) :: dspheat character(len=*), intent(out) :: errmsg diff --git a/physics/satmedmfvdifq.meta b/physics/satmedmfvdifq.meta index 88ab676b8..8538c6aa7 100644 --- a/physics/satmedmfvdifq.meta +++ b/physics/satmedmfvdifq.meta @@ -215,7 +215,7 @@ dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys - intent = out + intent = inout [u1] standard_name = x_wind long_name = x component of layer wind From d4d0b719ab297741ea484bee00e8285cca423a64 Mon Sep 17 00:00:00 2001 From: "samuel.trahan" Date: Tue, 24 May 2022 18:18:45 +0000 Subject: [PATCH 19/20] Simplify machine.F and remove unused types. --- physics/machine.F | 38 ++++++++++---------------------------- 1 file changed, 10 insertions(+), 28 deletions(-) diff --git a/physics/machine.F b/physics/machine.F index 9b09d235c..e9c572ef2 100644 --- a/physics/machine.F +++ b/physics/machine.F @@ -6,42 +6,24 @@ module machine implicit none -#ifndef SINGLE_PREC integer, parameter :: kind_io4 = 4, kind_io8 = 8 , kind_ior = 8 & &, kind_evod = 8, kind_dbl_prec = 8 & &, kind_sngl_prec = 4 -# ifdef __PGI - &, kind_qdt_prec = 8 & -# else - &, kind_qdt_prec = 16 & -# endif - &, kind_rad = 8 & - &, kind_phys = 8 ,kind_taum=8 & - &, kind_grid = 8 & - &, kind_REAL = 8 &! used in cmp_comm - &, kind_LOGICAL = 4 & - &, kind_INTEGER = 4 ! -,,- +#ifdef SINGLE_PREC + integer, parameter :: kind_rad = kind_sngl_prec & + &, kind_phys = kind_sngl_prec & + &, kind_grid = kind_dbl_prec &! atmos_cubed_sphere requres kind_grid=8 + &, kind_REAL = kind_sngl_prec ! used in cmp_comm #else - integer, parameter :: kind_io4 = 4, kind_io8 = 8 , kind_ior = 8 & - &, kind_evod = 4, kind_dbl_prec = 8 & - &, kind_sngl_prec = 4 -# ifdef __PGI - &, kind_qdt_prec = 8 & -# else - &, kind_qdt_prec = 16 & -# endif - &, kind_rad = 4 & - &, kind_phys = 4 ,kind_taum=4 & - &, kind_grid = 8 &! atmos_cubed_sphere requres kind_grid=8 - &, kind_REAL = 4 &! used in cmp_comm - &, kind_LOGICAL = 4 & - &, kind_INTEGER = 4 ! -,,- - + integer, parameter :: kind_rad = kind_dbl_prec & + &, kind_phys = kind_dbl_prec & + &, kind_grid = kind_dbl_prec &! atmos_cubed_sphere requres kind_grid=8 + &, kind_REAL = kind_dbl_prec ! used in cmp_comm #endif #ifdef OVERLOAD_R4 - integer, parameter :: kind_dyn = 4 + integer, parameter :: kind_dyn = 4 #else integer, parameter :: kind_dyn = 8 #endif From 942f9adcef364f463158c7e7a097a97b4ddb76f7 Mon Sep 17 00:00:00 2001 From: "samuel.trahan" Date: Wed, 25 May 2022 22:08:44 +0000 Subject: [PATCH 20/20] correct bug in machine.F --- physics/machine.F | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/physics/machine.F b/physics/machine.F index e9c572ef2..eb1dcd257 100644 --- a/physics/machine.F +++ b/physics/machine.F @@ -8,7 +8,8 @@ module machine integer, parameter :: kind_io4 = 4, kind_io8 = 8 , kind_ior = 8 & &, kind_evod = 8, kind_dbl_prec = 8 & - &, kind_sngl_prec = 4 + &, kind_sngl_prec = 4, kind_INTEGER = 4 & + &, kind_LOGICAL = 4 #ifdef SINGLE_PREC integer, parameter :: kind_rad = kind_sngl_prec &