From 7f31e4939fd64abb8b38e3d16a3a9ab9e0e1537b Mon Sep 17 00:00:00 2001 From: brianhenn Date: Fri, 16 Jun 2023 22:24:08 +0000 Subject: [PATCH 1/7] update fv3gfs-fortran to rad-progcld6 branch --- external/fv3gfs-fortran | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/external/fv3gfs-fortran b/external/fv3gfs-fortran index cb106f2eb8..e7695a697b 160000 --- a/external/fv3gfs-fortran +++ b/external/fv3gfs-fortran @@ -1 +1 @@ -Subproject commit cb106f2eb806e8c635d28d8b76ee8e80a0e20bc3 +Subproject commit e7695a697b6083d8d0595e88fa115e21a25edffc From 60be91cb08380249f09c76660d62c5594478a97f Mon Sep 17 00:00:00 2001 From: brianhenn Date: Mon, 19 Jun 2023 23:05:47 +0000 Subject: [PATCH 2/7] port progcld6 to python radiation and use namelist flag --- external/radiation/radiation/config.py | 2 + .../gfdl_cloud_microphys/__init__.py | 1 + .../gfdl_cloud_microphys/_cloud_eff_rad.py | 291 ++++++++++++++++++ .../gfdl_cloud_microphys/_constants.py | 210 +++++++++++++ .../radiation/radiation/radiation_clouds.py | 187 ++++++++++- .../radiation/radiation/radiation_driver.py | 70 +++-- external/radiation/radiation/wrapper_api.py | 1 + .../tests/test_radiation.py | 27 +- 8 files changed, 763 insertions(+), 26 deletions(-) create mode 100644 external/radiation/radiation/gfdl_cloud_microphys/__init__.py create mode 100644 external/radiation/radiation/gfdl_cloud_microphys/_cloud_eff_rad.py create mode 100644 external/radiation/radiation/gfdl_cloud_microphys/_constants.py diff --git a/external/radiation/radiation/config.py b/external/radiation/radiation/config.py index ef1d3afc2c..857e2640c4 100644 --- a/external/radiation/radiation/config.py +++ b/external/radiation/radiation/config.py @@ -66,6 +66,7 @@ class GFSPhysicsControlConfig: default_factory=lambda: [[-999.0], [-999.0], [-999.0], [-999.0], [-999.0]] ) do_only_clearsky_rad: bool = False + do_progcld6: bool = False swhtr: bool = True lwhtr: bool = True lprnt: bool = False @@ -93,6 +94,7 @@ def from_namelist(cls, namelist: Mapping[Hashable, Any]): "swhtr": "swhtr", "lwhtr": "lwhtr", "levr": "levr", + "rad_progcld6": "do_progcld6", } CORE_NAMELIST_TO_GFS_CONTROL = {"npz": "levs"} diff --git a/external/radiation/radiation/gfdl_cloud_microphys/__init__.py b/external/radiation/radiation/gfdl_cloud_microphys/__init__.py new file mode 100644 index 0000000000..9c1e007b4c --- /dev/null +++ b/external/radiation/radiation/gfdl_cloud_microphys/__init__.py @@ -0,0 +1 @@ +from ._cloud_eff_rad import cld_eff_rad diff --git a/external/radiation/radiation/gfdl_cloud_microphys/_cloud_eff_rad.py b/external/radiation/radiation/gfdl_cloud_microphys/_cloud_eff_rad.py new file mode 100644 index 0000000000..ee78f7ec8e --- /dev/null +++ b/external/radiation/radiation/gfdl_cloud_microphys/_cloud_eff_rad.py @@ -0,0 +1,291 @@ +from typing import Optional, Tuple +import numpy as np +import numba +from ._constants import ( + retab, + ccn_l, + ccn_o, + qcmin, + grav, + zvir, + rdgas, + pi, + rhow, + rewmin, + rewmax, + reimin, + reimax, + rermin, + rermax, + resmin, + resmax, + regmin, + regmax, + mur, + mus, + mug, + edar, + edas, + edag, + edbr, + edbs, + edbg, +) + +liq_ice_combine = False # combine all liquid water, combine all solid water +snow_graupel_combine = True # combine snow and graupel +prog_ccn = False # do prognostic ccn (yi ming's method) + +rewflag = 1 # cloud water effective radius scheme (only ported option is 1) +# 1: Martin et al. (1994) +# 2: Martin et al. (1994), GFDL revision +# 3: Kiehl et al. (1994) +# 4: effective radius from PSD + +reiflag = 4 # cloud ice effective radius scheme (only ported option is 4) +# 1: Heymsfield and Mcfarquhar (1996) +# 2: Donner et al. (1997) +# 3: Fu (2007) +# 4: Kristjansson et al. (2000) +# 5: Wyser (1998) +# 6: Sun and Rikus (1999), Sun (2001) +# 7: effective radius from PSD + +rerflag = 1 # rain effective radius scheme (only ported option is 1) +# 1: effective radius from PSD + +resflag = 1 # snow effective radius scheme (only ported option is 1) +# 1: effective radius from PSD + +regflag = 1 # graupel effective radius scheme (only ported option is 1) +# 1: effective radius from PSD + + +@numba.jit +def cld_eff_rad( + _is: int, + ie: int, + ks: int, + ke: int, + lsm: np.ndarray, + p: np.ndarray, + delp: np.ndarray, + t: np.ndarray, + qv: np.ndarray, + qw: np.ndarray, + qi: np.ndarray, + qr: np.ndarray, + qs: np.ndarray, + qg: np.ndarray, + qa: np.ndarray, + cloud: np.ndarray, + cnvw: Optional[np.ndarray] = None, + cnvi: Optional[np.ndarray] = None, + cnvc: Optional[np.ndarray] = None, +) -> Tuple[ + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, + np.ndarray, +]: + + qmw = qw + qmi = qi + qmr = qr + qms = qs + qmg = qg + + # outputs + qcw = np.zeros((ie, ke)) + qci = np.zeros((ie, ke)) + qcr = np.zeros((ie, ke)) + qcs = np.zeros((ie, ke)) + qcg = np.zeros((ie, ke)) + rew = np.zeros((ie, ke)) + rei = np.zeros((ie, ke)) + rer = np.zeros((ie, ke)) + res = np.zeros((ie, ke)) + reg = np.zeros((ie, ke)) + cld = cloud + + # ----------------------------------------------------------------------- + # merge convective cloud to total cloud + # ----------------------------------------------------------------------- + if cnvw is not None: + qmw = qmw + cnvw + if cnvi is not None: + qmi = qmi + cnvi + if cnvc is not None: + cld = cnvc + (1 - cnvc) * cld + + for i in range(_is, ie): + for k in range(ks, ke): + # ----------------------------------------------------------------------- + # combine liquid and solid phases + # ----------------------------------------------------------------------- + if liq_ice_combine: + qmw[i, k] = qmw[i, k] + qmr[i, k] + qmr[i, k] = 0.0 + qmi[i, k] = qmi[i, k] + qms[i, k] + qmg[i, k] + qms[i, k] = 0.0 + qmg[i, k] = 0.0 + + # ----------------------------------------------------------------------- + # combine snow and graupel + # ----------------------------------------------------------------------- + if snow_graupel_combine: + qms[i, k] = qms[i, k] + qmg[i, k] + qmg[i, k] = 0.0 + + qmw[i, k] = max(qmw[i, k], qcmin) + qmi[i, k] = max(qmi[i, k], qcmin) + qmr[i, k] = max(qmr[i, k], qcmin) + qms[i, k] = max(qms[i, k], qcmin) + qmg[i, k] = max(qmg[i, k], qcmin) + + cld[i, k] = min(max(cld[i, k], 0.0), 1.0) + + mask = min(max(lsm[i], 0.0), 2.0) + + dpg = np.abs(delp[i, k]) / grav + rho = p[i, k] / (rdgas * t[i, k] * (1.0 + zvir * qv[i, k])) + + if rewflag == 1: + + # ----------------------------------------------------------------------- + # cloud water (Martin et al. 1994) + # ----------------------------------------------------------------------- + + if prog_ccn: + # boucher and lohmann (1995) + ccnw = (1.0 - abs(mask - 1.0)) * ( + 10.0 ** 2.24 * (qa[i, k] * rho * 1.0e9) ** 0.257 + ) + abs(mask - 1.0) * ( + 10.0 ** 2.06 * (qa[i, k] * rho * 1.0e9) ** 0.48 + ) + else: + ccnw = ccn_o * abs(mask - 1.0) + ccn_l * (1.0 - abs(mask - 1.0)) + + if qmw[i, k] > qcmin: + qcw[i, k] = dpg * qmw[i, k] * 1.0e3 + rew[i, k] = ( + np.exp( + 1.0 + / 3.0 + * np.log((3.0 * qmw[i, k] * rho) / (4.0 * pi * rhow * ccnw)) + ) + * 1.0e4 + ) + rew[i, k] = max(rewmin, min(rewmax, rew[i, k])) + else: + qcw[i, k] = 0.0 + rew[i, k] = rewmin + + else: + raise ValueError( + f"Only ported option for `rewflag` is 1, got {rewflag}." + ) + + if reiflag == 4: + + # ----------------------------------------------------------------------- + # cloud ice (Kristjansson et al. 2000) + # ----------------------------------------------------------------------- + + if qmi[i, k] > qcmin: + qci[i, k] = dpg * qmi[i, k] * 1.0e3 + ind = min(max(int(t[i, k] - 136.0), 44), 138 - 1) + cor = t[i, k] - int(t[i, k]) + rei[i, k] = retab[ind] * (1.0 - cor) + retab[ind] * cor + rei[i, k] = max(reimin, min(reimax, rei[i, k])) + else: + qci[i, k] = 0.0 + rei[i, k] = reimin + else: + raise ValueError( + f"Only ported option for `reiflag` is 4, got {reiflag}." + ) + + if rerflag == 1: + + # ----------------------------------------------------------------------- + # rain derived from PSD + # ----------------------------------------------------------------------- + + if qmr[i, k] > qcmin: + qcr[i, k] = dpg * qmr[i, k] * 1.0e3 + der = calc_ed(qmr[i, k], rho, mur, eda=edar, edb=edbr) + rer[i, k] = 0.5 * der * 1.0e6 + rer[i, k] = max(rermin, min(rermax, rer[i, k])) + else: + qcr[i, k] = 0.0 + rer[i, k] = rermin + else: + raise ValueError( + f"Only ported option for `rerflag` is 4, got {rerflag}." + ) + + if resflag == 1: + + # ----------------------------------------------------------------------- + # snow derived from PSD + # ----------------------------------------------------------------------- + + if qms[i, k] > qcmin: + qcs[i, k] = dpg * qms[i, k] * 1.0e3 + des = calc_ed(qms[i, k], rho, mus, eda=edas, edb=edbs) + res[i, k] = 0.5 * des * 1.0e6 + res[i, k] = max(resmin, min(resmax, res[i, k])) + else: + qcs[i, k] = 0.0 + res[i, k] = resmin + else: + raise ValueError( + f"Only ported option for `resflag` is 1, got {resflag}." + ) + + if regflag == 1: + + # ----------------------------------------------------------------------- + # graupel derived from PSD + # ----------------------------------------------------------------------- + + if qmg[i, k] > qcmin: + qcg[i, k] = dpg * qmg[i, k] * 1.0e3 + deg = calc_ed(qmg[i, k], rho, mug, eda=edag, edb=edbg) + reg[i, k] = 0.5 * deg * 1.0e6 + reg[i, k] = max(regmin, min(regmax, rer[i, k])) + else: + qcg[i, k] = 0.0 + reg[i, k] = regmin + else: + raise ValueError( + f"Only ported option for `regflag` is 4, got {regflag}." + ) + + return ( + qcw, + qci, + qcr, + qcs, + qcg, + rew, + rei, + rer, + res, + reg, + cld, + ) + + +def calc_ed(q: float, den: float, mu: float, eda: float, edb: float) -> float: + """calculation of effective diameter (ed)""" + ed = eda / edb * np.exp(1.0 / (mu + 3) * np.log(6 * den * q)) + return ed diff --git a/external/radiation/radiation/gfdl_cloud_microphys/_constants.py b/external/radiation/radiation/gfdl_cloud_microphys/_constants.py new file mode 100644 index 0000000000..487548f926 --- /dev/null +++ b/external/radiation/radiation/gfdl_cloud_microphys/_constants.py @@ -0,0 +1,210 @@ +import numpy as np +from scipy.special import gamma + +grav = 9.80665 # gfs: acceleration due to gravity +rhow = 1.0e3 # density of cloud water (kg/m^3) +rhoi = 9.17e2 # density of cloud ice (kg/m^3) +rhor = 1.0e3 # density of rain (kg/m^3) +rhos = 0.1e3 # density of cloud snow (kg/m^3) +rhog = 0.4e3 # density of cloud graupel (kg/m^3) +rdgas = 287.05 # gfs: gas constant for dry air +rvgas = 461.50 # gfs: gas constant for water vapor +zvir = rvgas / rdgas - 1.0 # 0.6077338443 +pi = 3.1415926535897931 # gfs: ratio of circle circumference to diameter + +ccn_o = 90.0 # ccn over ocean (cm^ - 3) +ccn_l = 270.0 # ccn over land (cm^ - 3) + +qcmin = 1.0e-12 # min value for cloud condensates + +rewmin = 5.0 +rewmax = 15.0 +reimin = 10.0 +reimax = 150.0 +rermin = 15.0 +rermax = 10000.0 +resmin = 150.0 +resmax = 10000.0 +regmin = 150.0 +regmax = 10000.0 + +n0r_sig = 8.0 # intercept parameter (significand) of rain (Lin et al. 1983) (1/m^4) +# (Marshall and Palmer 1948) +n0s_sig = 3.0 # intercept parameter (significand) of snow (Lin et al. 1983) (1/m^4) +# (Gunn and Marshall 1958) +n0g_sig = 4.0 # intercept parameter (significand) of graupel (Rutledge and Hobbs 1984) +# (1/m^4) (Houze et al. 1979) + +n0r_exp = 6 # intercept parameter (exponent) of rain (Lin et al. 1983) (1/m^4) +# (Marshall and Palmer 1948) +n0s_exp = 6 # intercept parameter (exponent) of snow (Lin et al. 1983) (1/m^4) +# (Gunn and Marshall 1958) +n0g_exp = ( + 6 # intercept parameter (exponent) of graupel (Rutledge and Hobbs 1984) (1/m^4) +) +# (Houze et al. 1979) + +mur = 1.0 # shape parameter of rain in Gamma distribution (Marshall and Palmer 1948) +mus = 1.0 # shape parameter of snow in Gamma distribution (Gunn and Marshall 1958) +mug = 1.0 # shape parameter of graupel in Gamma distribution (Houze et al. 1979) + +edar = ( + np.exp(-1.0 / (mur + 3) * np.log(n0r_sig)) + * (mur + 2) + * np.exp(-n0r_exp / (mur + 3) * np.log(10.0)) +) +edas = ( + np.exp(-1.0 / (mus + 3) * np.log(n0s_sig)) + * (mus + 2) + * np.exp(-n0s_exp / (mus + 3) * np.log(10.0)) +) +edag = ( + np.exp(-1.0 / (mug + 3) * np.log(n0g_sig)) + * (mug + 2) + * np.exp(-n0g_exp / (mug + 3) * np.log(10.0)) +) + +edbr = np.exp(1.0 / (mur + 3) * np.log(pi * rhor * gamma(mur + 3))) +edbs = np.exp(1.0 / (mus + 3) * np.log(pi * rhos * gamma(mus + 3))) +edbg = np.exp(1.0 / (mug + 3) * np.log(pi * rhog * gamma(mug + 3))) + +retab = [ + 0.05000, + 0.05000, + 0.05000, + 0.05000, + 0.05000, + 0.05000, + 0.05500, + 0.06000, + 0.07000, + 0.08000, + 0.09000, + 0.10000, + 0.20000, + 0.30000, + 0.40000, + 0.50000, + 0.60000, + 0.70000, + 0.80000, + 0.90000, + 1.00000, + 1.10000, + 1.20000, + 1.30000, + 1.40000, + 1.50000, + 1.60000, + 1.80000, + 2.00000, + 2.20000, + 2.40000, + 2.60000, + 2.80000, + 3.00000, + 3.20000, + 3.50000, + 3.80000, + 4.10000, + 4.40000, + 4.70000, + 5.00000, + 5.30000, + 5.60000, + 5.92779, + 6.26422, + 6.61973, + 6.99539, + 7.39234, + 7.81177, + 8.25496, + 8.72323, + 9.21800, + 9.74075, + 10.2930, + 10.8765, + 11.4929, + 12.1440, + 12.8317, + 13.5581, + 14.2319, + 15.0351, + 15.8799, + 16.7674, + 17.6986, + 18.6744, + 19.6955, + 20.7623, + 21.8757, + 23.0364, + 24.2452, + 25.5034, + 26.8125, + 27.7895, + 28.6450, + 29.4167, + 30.1088, + 30.7306, + 31.2943, + 31.8151, + 32.3077, + 32.7870, + 33.2657, + 33.7540, + 34.2601, + 34.7892, + 35.3442, + 35.9255, + 36.5316, + 37.1602, + 37.8078, + 38.4720, + 39.1508, + 39.8442, + 40.5552, + 41.2912, + 42.0635, + 42.8876, + 43.7863, + 44.7853, + 45.9170, + 47.2165, + 48.7221, + 50.4710, + 52.4980, + 54.8315, + 57.4898, + 60.4785, + 63.7898, + 65.5604, + 71.2885, + 75.4113, + 79.7368, + 84.2351, + 88.8833, + 93.6658, + 98.5739, + 103.603, + 108.752, + 114.025, + 119.424, + 124.954, + 130.630, + 136.457, + 142.446, + 148.608, + 154.956, + 161.503, + 168.262, + 175.248, + 182.473, + 189.952, + 197.699, + 205.728, + 214.055, + 222.694, + 231.661, + 240.971, + 250.639, +] diff --git a/external/radiation/radiation/radiation_clouds.py b/external/radiation/radiation/radiation_clouds.py index ad3b63f555..17d6f5d298 100644 --- a/external/radiation/radiation/radiation_clouds.py +++ b/external/radiation/radiation/radiation_clouds.py @@ -1,4 +1,5 @@ import numpy as np +from .gfdl_cloud_microphys import cld_eff_rad from .phys_const import con_ttp, con_pi, con_g, con_rd, con_thgni from .radphysparam import lcrick @@ -7,7 +8,7 @@ class CloudClass: VTAGCLD = "NCEP-Radiation_clouds v5.1 Nov 2012 " gfac = 1.0e5 / con_g gord = con_g / con_rd - NF_CLDS = 9 + NF_CLDS = 11 NK_CLDS = 3 ptopc = np.array([[1050.0, 650.0, 400.0, 0.0], [1050.0, 750.0, 500.0, 0.0]]).T climit = 0.001 @@ -1460,6 +1461,190 @@ def progclduni( return clouds, clds, mtop, mbot, de_lgth + def progcld6( + self, + plyr, + tlyr, + tvly, + qlyr, + qstl, + rhly, + cnvw, + cnvc, + xlat, + qw, + qi, + qr, + qs, + qg, + qa, + slmsk, + cldtot, + dz, + delp, + IX, + NLAY, + NLP1, + ): + # ================= subprogram documentation block ================ ! + # ! + # subprogram: progcld6 computes cloud related quantities using ! + # GFDL Lin MP prognostic cloud microphysics scheme, as run in ! + # SHiELD physics ! + # ! + # abstract: this program computes cloud fractions from cloud ! + # condensates, calculates liquid/ice cloud droplet effective radius, ! + # and computes the low, mid, high, total and boundary layer cloud ! + # fractions and the vertical indices of low, mid, and high cloud ! + # top and base. the three vertical cloud domains are set up in the ! + # initial subroutine "cld_init". ! + # ! + # usage: call progcld6 ! + # ! + # subprograms called: cld_eff_rad ! + # ! + # attributes: ! + # language: fortran 90 ! + # machine: ibm-sp, sgi ! + # ! + # ! + # ==================== definition of variables ==================== ! + # ! + # input variables: ! + # plyr (IX,NLAY) : model layer mean pressure in mb (100Pa) ! + # plvl (IX,NLP1) : model level pressure in mb (100Pa) ! + # tlyr (IX,NLAY) : model layer mean temperature in k ! + # tvly (IX,NLAY) : model layer virtual temperature in k ! + # qlyr (IX,NLAY) : layer specific humidity in gm/gm ! + # qstl (IX,NLAY) : layer saturate humidity in gm/gm ! + # rhly (IX,NLAY) : layer relative humidity (=qlyr/qstl) ! + # cnvw (IX,NLAY) : layer convective cloud condensate ! + # cnvc (IX,NLAY) : layer convective cloud cover ! + # xlat (IX) : grid latitude in radians, default to pi/2 -> -pi/2! + # range, otherwise see in-line comment ! + # qw (IX,NLAY) : liquid cloud mixing ratio in kg/kg ! + # qi (IX,NLAY) : ice cloud mixing ratio in kg/kg ! + # qr (IX,NLAY) : rain mixing ratio in kg/kg ! + # qs (IX,NLAY) : snow mixing ratio in kg/kg ! + # qg (IX,NLAY) : graupel mixing ratio in kg/kg ! + # qa (IX,NLAY) : aerosol mixing ratio ! + # slmsk (IX) : sea/land mask array (sea:0,land:1,sea-ice:2) ! + # cldtot (ITX) : cloud area fraction ! + # dz (ix,nlay) : layer thickness (km) ! + # delp (ix,nlay) : model layer pressure thickness in mb (100Pa) ! + # IX : horizontal dimention ! + # NLAY,NLP1 : vertical layer/level dimensions ! + # ! + # output variables: ! + # clouds(IX,NLAY,NF_CLDS) : cloud profiles ! + # clouds(:,:,1) - layer total cloud fraction ! + # clouds(:,:,2) - layer cloud liq water path (g/m**2) ! + # clouds(:,:,3) - mean eff radius for liq cloud (micron) ! + # clouds(:,:,4) - layer cloud ice water path (g/m**2) ! + # clouds(:,:,5) - mean eff radius for ice cloud (micron) ! + # clouds(:,:,6) - layer rain drop water path (g/m**2) ! + # clouds(:,:,7) - mean eff radius for rain drop (micron) ! + # *** clouds(:,:,8) - layer snow flake water path (g/m**2) ! + # clouds(:,:,9) - mean eff radius for snow flake (micron) ! + # clouds(:,:,10)- layer graupel water path (g/m**2) ! + # clouds(:,:,11)- mean eff radius for graupel (micron) ! + # *** fu's scheme need to be normalized by snow density (g/m**3/1.0e6) ! + # clds (IX,5) : fraction of clouds for low, mid, hi, tot, bl ! + # mtop (IX,3) : vertical indices for low, mid, hi cloud tops ! + # mbot (IX,3) : vertical indices for low, mid, hi cloud bases ! + # de_lgth(ix) : clouds decorrelation length (km) ! + # ! + # module variables: ! + # ivflip : control flag of vertical index direction ! + # =0: index from toa to surface ! + # =1: index from surface to toa ! + # lsashal : control flag for shallow convection ! + # lcrick : control flag for eliminating CRICK ! + # =t: apply layer smoothing to eliminate CRICK ! + # =f: do not apply layer smoothing ! + # lcnorm : control flag for in-cld condensate ! + # =t: normalize cloud condensate ! + # =f: not normalize cloud condensate ! + # ! + # ==================== end of description ===================== ! + # + + cldcnv = np.zeros((IX, NLAY)) + ptop1 = np.zeros((IX, self.NK_CLDS + 1)) + rxlat = np.zeros(IX) + clouds = np.zeros((IX, NLAY, self.NF_CLDS)) + de_lgth = np.zeros(IX) + + (cwp, cip, crp, csp, cgp, rew, rei, rer, res, reg, cldtot) = cld_eff_rad( + 0, + IX, + 0, + NLAY, + slmsk, + plyr, + delp, + tlyr, + qlyr, + qw, + qi, + qr, + qs, + qg, + qa, + cldtot, + cnvw=cnvw, + cnvc=cnvc, + ) + cldcnv = 0.0 + + for i in range(IX): + rxlat[i] = np.abs(xlat[i] / con_pi) # if xlat in pi/2 -> -pi/2 range + + # --- find top pressure for each cloud domain for given latitude + # ptopc(k,i): top presure of each cld domain (k=1-4 are sfc,L,m,h; + # --- i=1,2 are low-lat (<45 degree) and pole regions) + + for id in range(4): + tem1 = self.ptopc(id, 1) - self.ptopc(id, 0) + for i in range(IX): + tem2 = xlat[i] / con_pi # if xlat in pi/2 -> -pi/2 range + ptop1[i, id] = self.ptopc(id, 1) + tem1 * max( + 0.0, 4.0 * abs(tem2) - 1.0 + ) + + for k in range(NLAY): + for i in range(IX): + clouds[i, k, 0] = cldtot[i, k] + clouds[i, k, 1] = cwp[i, k] + clouds[i, k, 2] = rew[i, k] + clouds[i, k, 3] = cip[i, k] + clouds[i, k, 4] = rei[i, k] + clouds[i, k, 5] = crp[i, k] + clouds[i, k, 6] = rer[i, k] + clouds[i, k, 7] = csp[i, k] + clouds[i, k, 8] = res[i, k] + clouds[i, k, 9] = cgp[i, k] + clouds[i, k, 10] = reg[i, k] + + # --- ... estimate clouds decorrelation length in km + # this is only a tentative test, need to consider change later + + if self.iovr == 3: + for i in range(IX): + de_lgth[i] = max(0.6, 2.78 - 4.6 * rxlat[i]) + + # --- compute low, mid, high, total, and boundary layer cloud fractions + # and clouds top/bottom layer indices for low, mid, and high clouds. + # The three cloud domain boundaries are defined by ptopc. The cloud + # overlapping method is defined by control flag 'iovr', which may + # be different for lw and sw radiation programs. + + clds, mtop, mbot = self.gethml( + plyr, ptop1, cldtot, cldcnv, dz, de_lgth, IX, NLAY + ) + + return clouds, clds, mtop, mbot, de_lgth + def gethml(self, plyr, ptop1, cldtot, cldcnv, dz, de_lgth, IX, NLAY): # =================================================================== ! # ! diff --git a/external/radiation/radiation/radiation_driver.py b/external/radiation/radiation/radiation_driver.py index ff2bb85289..c72b742a4f 100644 --- a/external/radiation/radiation/radiation_driver.py +++ b/external/radiation/radiation/radiation_driver.py @@ -825,27 +825,55 @@ def _GFS_radiation_driver( if gfs_physics_control.config.imp_physics == 99: ccnd[:IM, :LMK, 0] = ccnd[:IM, :LMK, 0] + cnvw[:IM, :LMK] - clouds, cldsa, mtopa, mbota, de_lgth = self.cld.progcld4( - plyr, - plvl, - tlyr, - tvly, - qlyr, - qstl, - rhly, - ccnd[:IM, :LMK, 0], - cnvw, - cnvc, - Grid["xlat"], - Grid["xlon"], - Sfcprop["slmsk"], - cldcov, - dz, - delp, - IM, - LMK, - LMP, - ) + if not gfs_physics_control.do_progcld6: + clouds, cldsa, mtopa, mbota, de_lgth = self.cld.progcld4( + plyr, + plvl, + tlyr, + tvly, + qlyr, + qstl, + rhly, + ccnd[:IM, :LMK, 0], + cnvw, + cnvc, + Grid["xlat"], + Grid["xlon"], + Sfcprop["slmsk"], + cldcov, + dz, + delp, + IM, + LMK, + LMP, + ) + else: + q_aerosol = np.zeros((IM, LMK)) + clouds, cldsa, mtopa, mbota, de_lgth = self.cld.progcld6( + plyr, + tlyr, + tvly, + qlyr, + qstl, + rhly, + ccnd[:IM, :LMK, 0], + cnvw, + cnvc, + Grid["xlat"], + tracer1[:IM, :LMK, ntcw - 1], + tracer1[:IM, :LMK, ntiw - 1], + tracer1[:IM, :LMK, ntrw - 1], + tracer1[:IM, :LMK, ntsw - 1], + tracer1[:IM, :LMK, ntgl - 1], + q_aerosol, + Sfcprop["slmsk"], + cldcov, + dz, + delp, + IM, + LMK, + LMP, + ) # --- ... start radiation calculations # remember to set heating rate unit to k/sec! diff --git a/external/radiation/radiation/wrapper_api.py b/external/radiation/radiation/wrapper_api.py index 5baaabfae0..052d03915c 100644 --- a/external/radiation/radiation/wrapper_api.py +++ b/external/radiation/radiation/wrapper_api.py @@ -73,6 +73,7 @@ class GFSPhysicsControl: nslwr: int lsswr: bool = True lslwr: bool = True + do_progcld6: bool = False @classmethod def from_config( diff --git a/workflows/prognostic_c48_run/tests/test_radiation.py b/workflows/prognostic_c48_run/tests/test_radiation.py index cdc265c93e..c73279f519 100644 --- a/workflows/prognostic_c48_run/tests/test_radiation.py +++ b/workflows/prognostic_c48_run/tests/test_radiation.py @@ -1,3 +1,4 @@ +from collections.abc import Hashable, Mapping from runtime.segmented_run.prepare_config import HighLevelConfig from runtime.segmented_run.api import create, append from runtime.steppers.radiation import RadiationStepper @@ -64,6 +65,7 @@ def python_name(self): fhswr: 1800.0 hybedmf: true satmedmf: false + rad_progcld6: false fv_core_nml: npx: 13 npy: 13 @@ -105,8 +107,25 @@ def get_fv3config(): return fv3config_dict -def radiation_scheme_config(): +NAMELIST_UPDATES = { + "base": {}, + "progcld6": {"namelist": {"gfs_physics_nml": {"rad_progcld6": True}}}, +} + + +def _update_nested(config, updates): + for key, update in updates.items(): + if isinstance(config[key], Mapping) and isinstance(update, Mapping): + _update_nested(config[key], update) + elif isinstance(config[key], Hashable) and isinstance(update, Hashable): + config[key] = update + else: + raise ValueError("Invalid update to configuration dict structure.") + + +def radiation_scheme_config(version): config = get_fv3config() + _update_nested(config, NAMELIST_UPDATES[version]) config["radiation_scheme"] = {"kind": "python"} diagnostics = [ { @@ -139,9 +158,9 @@ def run_model(config, rundir): append(rundir) -@pytest.fixture(scope="module") -def completed_rundir(tmpdir_factory): - config = radiation_scheme_config() +@pytest.fixture(scope="module", params=["base", "progcld6"]) +def completed_rundir(tmpdir_factory, request): + config = radiation_scheme_config(request.param) rundir = tmpdir_factory.mktemp("rundir").join("subdir") run_model(config, str(rundir)) return rundir From 415d96a111ae385c15eb56ebc6a276e0aa310ba5 Mon Sep 17 00:00:00 2001 From: brianhenn Date: Mon, 19 Jun 2023 23:06:24 +0000 Subject: [PATCH 3/7] update to fortran branch with slight fixes --- external/fv3gfs-fortran | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/external/fv3gfs-fortran b/external/fv3gfs-fortran index e7695a697b..ad97722e44 160000 --- a/external/fv3gfs-fortran +++ b/external/fv3gfs-fortran @@ -1 +1 @@ -Subproject commit e7695a697b6083d8d0595e88fa115e21a25edffc +Subproject commit ad97722e446c19794adb76109387632b330ea868 From a0627039e4132bec2615336b3381d01675876796 Mon Sep 17 00:00:00 2001 From: brianhenn Date: Mon, 19 Jun 2023 23:22:13 +0000 Subject: [PATCH 4/7] fix radiation port validation test --- external/radiation/tests/test_driver/test_driver.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/external/radiation/tests/test_driver/test_driver.py b/external/radiation/tests/test_driver/test_driver.py index af4ea88b9b..3ec2407b44 100644 --- a/external/radiation/tests/test_driver/test_driver.py +++ b/external/radiation/tests/test_driver/test_driver.py @@ -33,13 +33,13 @@ def getscalars(indict): def get_gfs_physics_control(indict): config_names = [field.name for field in dataclasses.fields(GFSPhysicsControlConfig)] - config_kwargs = {name: indict[name] for name in config_names} + config_kwargs = {name: indict[name] for name in config_names if name in indict} config = GFSPhysicsControlConfig(**config_kwargs) indict_ = indict.copy() indict_["config"] = config indict_["nsswr"], indict_["nslwr"] = 1, 1 control_names = [field.name for field in dataclasses.fields(GFSPhysicsControl)] - kwargs = {name: indict_[name] for name in control_names} + kwargs = {name: indict_[name] for name in control_names if name in indict_} return GFSPhysicsControl(**kwargs) From 4f5438c29236818ba1f3a7027af84bc4f36ef538 Mon Sep 17 00:00:00 2001 From: brianhenn Date: Tue, 20 Jun 2023 06:54:46 +0000 Subject: [PATCH 5/7] debug progcld6 --- .../radiation/gfdl_cloud_microphys/_cloud_eff_rad.py | 2 -- external/radiation/radiation/radiation_clouds.py | 9 ++------- external/radiation/radiation/radiation_driver.py | 7 +------ external/radiation/radiation/wrapper_api.py | 1 - 4 files changed, 3 insertions(+), 16 deletions(-) diff --git a/external/radiation/radiation/gfdl_cloud_microphys/_cloud_eff_rad.py b/external/radiation/radiation/gfdl_cloud_microphys/_cloud_eff_rad.py index ee78f7ec8e..d18de52470 100644 --- a/external/radiation/radiation/gfdl_cloud_microphys/_cloud_eff_rad.py +++ b/external/radiation/radiation/gfdl_cloud_microphys/_cloud_eff_rad.py @@ -1,6 +1,5 @@ from typing import Optional, Tuple import numpy as np -import numba from ._constants import ( retab, ccn_l, @@ -61,7 +60,6 @@ # 1: effective radius from PSD -@numba.jit def cld_eff_rad( _is: int, ie: int, diff --git a/external/radiation/radiation/radiation_clouds.py b/external/radiation/radiation/radiation_clouds.py index 17d6f5d298..cf0e8e1f58 100644 --- a/external/radiation/radiation/radiation_clouds.py +++ b/external/radiation/radiation/radiation_clouds.py @@ -1465,10 +1465,7 @@ def progcld6( self, plyr, tlyr, - tvly, qlyr, - qstl, - rhly, cnvw, cnvc, xlat, @@ -1484,7 +1481,6 @@ def progcld6( delp, IX, NLAY, - NLP1, ): # ================= subprogram documentation block ================ ! # ! @@ -1595,7 +1591,6 @@ def progcld6( cnvw=cnvw, cnvc=cnvc, ) - cldcnv = 0.0 for i in range(IX): rxlat[i] = np.abs(xlat[i] / con_pi) # if xlat in pi/2 -> -pi/2 range @@ -1605,10 +1600,10 @@ def progcld6( # --- i=1,2 are low-lat (<45 degree) and pole regions) for id in range(4): - tem1 = self.ptopc(id, 1) - self.ptopc(id, 0) + tem1 = self.ptopc[id, 1] - self.ptopc[id, 0] for i in range(IX): tem2 = xlat[i] / con_pi # if xlat in pi/2 -> -pi/2 range - ptop1[i, id] = self.ptopc(id, 1) + tem1 * max( + ptop1[i, id] = self.ptopc[id, 0] + tem1 * max( 0.0, 4.0 * abs(tem2) - 1.0 ) diff --git a/external/radiation/radiation/radiation_driver.py b/external/radiation/radiation/radiation_driver.py index c72b742a4f..e34348ff04 100644 --- a/external/radiation/radiation/radiation_driver.py +++ b/external/radiation/radiation/radiation_driver.py @@ -825,7 +825,7 @@ def _GFS_radiation_driver( if gfs_physics_control.config.imp_physics == 99: ccnd[:IM, :LMK, 0] = ccnd[:IM, :LMK, 0] + cnvw[:IM, :LMK] - if not gfs_physics_control.do_progcld6: + if not gfs_physics_control.config.do_progcld6: clouds, cldsa, mtopa, mbota, de_lgth = self.cld.progcld4( plyr, plvl, @@ -852,11 +852,7 @@ def _GFS_radiation_driver( clouds, cldsa, mtopa, mbota, de_lgth = self.cld.progcld6( plyr, tlyr, - tvly, qlyr, - qstl, - rhly, - ccnd[:IM, :LMK, 0], cnvw, cnvc, Grid["xlat"], @@ -872,7 +868,6 @@ def _GFS_radiation_driver( delp, IM, LMK, - LMP, ) # --- ... start radiation calculations diff --git a/external/radiation/radiation/wrapper_api.py b/external/radiation/radiation/wrapper_api.py index 052d03915c..5baaabfae0 100644 --- a/external/radiation/radiation/wrapper_api.py +++ b/external/radiation/radiation/wrapper_api.py @@ -73,7 +73,6 @@ class GFSPhysicsControl: nslwr: int lsswr: bool = True lslwr: bool = True - do_progcld6: bool = False @classmethod def from_config( From 933bc3bf7a22a1e582783b5199e4da75b0c72826 Mon Sep 17 00:00:00 2001 From: brianhenn Date: Wed, 21 Jun 2023 00:14:37 +0000 Subject: [PATCH 6/7] add ccnorm in port progcld6 and bump fortran to same --- external/fv3gfs-fortran | 2 +- external/radiation/radiation/radiation_clouds.py | 15 +++++++++++++-- 2 files changed, 14 insertions(+), 3 deletions(-) diff --git a/external/fv3gfs-fortran b/external/fv3gfs-fortran index ad97722e44..a34355e834 160000 --- a/external/fv3gfs-fortran +++ b/external/fv3gfs-fortran @@ -1 +1 @@ -Subproject commit ad97722e446c19794adb76109387632b330ea868 +Subproject commit a34355e83437dd54062b5b9a7880736015b137a7 diff --git a/external/radiation/radiation/radiation_clouds.py b/external/radiation/radiation/radiation_clouds.py index cf0e8e1f58..423053cbc0 100644 --- a/external/radiation/radiation/radiation_clouds.py +++ b/external/radiation/radiation/radiation_clouds.py @@ -1571,14 +1571,25 @@ def progcld6( clouds = np.zeros((IX, NLAY, self.NF_CLDS)) de_lgth = np.zeros(IX) + if self.lcnorm is True: + for k in range(NLAY): + for i in range(IX): + if cldtot[i, k] >= self.climit: + tem1 = 1.0 / max(self.climit2, cldtot[i, k]) + qw[i, k] = qw[i, k] * tem1 + qi[i, k] = qi[i, k] * tem1 + qr[i, k] = qr[i, k] * tem1 + qs[i, k] = qs[i, k] * tem1 + qg[i, k] = qg[i, k] * tem1 + (cwp, cip, crp, csp, cgp, rew, rei, rer, res, reg, cldtot) = cld_eff_rad( 0, IX, 0, NLAY, slmsk, - plyr, - delp, + 100 * plyr, + 100 * delp, tlyr, qlyr, qw, From 7c65b67ad44aab363144dec5a3dcd5ded6d17c09 Mon Sep 17 00:00:00 2001 From: brianhenn Date: Wed, 21 Jun 2023 00:21:43 +0000 Subject: [PATCH 7/7] use correct fv3gfs-fortran --- external/fv3gfs-fortran | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/external/fv3gfs-fortran b/external/fv3gfs-fortran index a34355e834..fe4749cbd9 160000 --- a/external/fv3gfs-fortran +++ b/external/fv3gfs-fortran @@ -1 +1 @@ -Subproject commit a34355e83437dd54062b5b9a7880736015b137a7 +Subproject commit fe4749cbd9c8fe07ef4cb484c4f901a628129ab1