The bug of in-cloud ice water path (ciwp) in RRTMG

[lizhuo1108]

I copied the example code from:
https://github.com/climlab/climlab-rrtmg/blob/main/climlab_rrtmg/tests/test_climlab_rrtmg.py
to a jupyter notebook. I changed the value of ciwp from 0 to 0.1 (Line 67 in the example code) and the kernel died. Maybe there is a bug on ciwp? Has someone else run into similar problems before?

Below is my code:
import numpy as np
from climlab_rrtmg import rrtmg_lw, rrtmg_sw
import matplotlib.pyplot as plt

cp = 1004. # Specific heat at constant pressure

ps = 1000. # Surface pressure in hPa

ncol = 1
nlay = 30
deltap = ps / nlay # pressure interval
plev = np.linspace(ps, 0., nlay+1) # pressure bounds
plev = plev[np.newaxis, ...]
play = np.linspace(ps-deltap/2., deltap/2., nlay)
play = play[np.newaxis, ...]

tsfc = 288.
tlay = np.linspace(288.-10., 200., nlay)
tlay = tlay[np.newaxis, ...]
tlev = np.linspace(288., 200., nlay+1)
tlev = tlev[np.newaxis, ...]

specific_humidity = np.array([4.13141097e-03, 3.41509495e-03, 2.81099479e-03, 2.30359570e-03,
1.87920573e-03, 1.52578624e-03, 1.23279279e-03, 9.91026303e-04,
7.92494475e-04, 6.30283118e-04, 4.98437246e-04, 3.91851633e-04,
3.06170488e-04, 2.37695932e-04, 1.83304857e-04, 1.40373783e-04,
1.06711275e-04, 8.04974602e-05, 6.02302082e-05, 4.46774859e-05,
3.28354282e-05, 2.38916443e-05, 1.71932832e-05, 1.22193649e-05,
8.55682965e-06, 5.87957411e-06, 5.00000000e-06, 5.00000000e-06,
5.00000000e-06, 5.00000000e-06])

h2ovmr = specific_humidity * 28.97 / 18.01528
h2ovmr = h2ovmr[np.newaxis, ...]

o3vmr = np.array([2.25573888e-08, 2.38730436e-08, 2.52586476e-08, 2.66442517e-08,
2.80298557e-08, 2.97254145e-08, 3.14254923e-08, 3.31238355e-08,
3.46767916e-08, 3.62297478e-08, 3.76122833e-08, 3.86410454e-08,
3.96698075e-08, 4.08899052e-08, 4.21303310e-08, 4.39781220e-08,
4.60528063e-08, 4.87636254e-08, 5.16974065e-08, 5.57122567e-08,
6.17914190e-08, 7.15771368e-08, 9.29020109e-08, 1.29109217e-07,
1.75914529e-07, 2.45552383e-07, 3.92764464e-07, 7.61726407e-07,
2.25137178e-06, 7.27500161e-06])
o3vmr = o3vmr[np.newaxis, ...]
co2vmr = 348. / 1E6 * np.ones_like(play)
ch4vmr = 1650. / 1E9 * np.ones_like(play)
n2ovmr = 306. / 1E9 * np.ones_like(play)
o2vmr = 0.21 * np.ones_like(play)
cfc11vmr = 0. * np.ones_like(play)
cfc12vmr = 0. * np.ones_like(play)
cfc22vmr = 0. * np.ones_like(play)
ccl4vmr = 0. * np.ones_like(play)

cloud_level_index = 8
cldfrac = 0.5np.exp(-(play-play[0,cloud_level_index])**2/(225.)**2) # Layer cloud fraction: a Gaussian centered on a pressure level
clwp = 60. * np.ones_like(play) # in-cloud liquid water path (g/m2)
ciwp = 0.1 * np.ones_like(play) # in-cloud ice water path (g/m2)
relq = 14. * np.ones_like(play) # Cloud water drop effective radius (microns)
reic = 0 * np.ones_like(play) # Cloud ice particle effective size (microns)`

nbndlw = int(rrtmg_lw.parrrtm.nbndlw)
ngptlw = int(rrtmg_lw.parrrtm.ngptlw)
rrtmg_lw.climlab_rrtmg_lw_ini(cp)

icld = 1 # Cloud overlap method, 0: Clear only, 1: Random, 2, Maximum/random] 3: Maximum
irng = 1 # more monte carlo stuff
idrv = 0 # whether to also calculate the derivative of flux with respect to surface temp
permuteseed = 300
inflglw = 2
iceflglw = 1
liqflglw = 1

tauc = 0. * np.ones_like(play)
tauc = tauc * np.ones([nbndlw,ncol,nlay])
tauaer = 0. * np.ones_like(play)
tauaer = np.transpose(tauaer * np.ones([nbndlw,ncol,nlay]), (1,2,0))

emis = 1. * np.ones((ncol,nbndlw))

(cldfmcl, ciwpmcl, clwpmcl, reicmcl, relqmcl, taucmcl) =
rrtmg_lw.climlab_mcica_subcol_lw(
ncol, nlay, icld,
permuteseed, irng, play,
cldfrac, ciwp, clwp, reic, relq, tauc)

ispec=0

(olr_sr, uflx, dflx, hr, uflxc, dflxc, hrc, duflx_dt, duflxc_dt) =
rrtmg_lw.climlab_rrtmg_lw(ncol, nlay, icld, ispec, idrv,
play, plev, tlay, tlev, tsfc,
h2ovmr, o3vmr, co2vmr, ch4vmr, n2ovmr, o2vmr,
cfc11vmr, cfc12vmr, cfc22vmr, ccl4vmr, emis,
inflglw, iceflglw, liqflglw, cldfmcl,
taucmcl, ciwpmcl, clwpmcl, reicmcl, relqmcl,
tauaer)

[brian-rose]

Hi @lizhuo1108, I think your bug results from setting zero effective ice crystal radius reic.

Try changing reic to a realistic value in your code. Or, here is a series of working steps that elaborate on that example:

Set up a single column

import numpy as np
from climlab_rrtmg import rrtmg_lw, rrtmg_sw

# Specific heat at constant pressure
cp = 1004.

# Set up the pressure domain
ps = 1000. #  Surface pressure in hPa
# RRTM code expects arrays with (ncol, nlay)
# and with pressure decreasing from surface at element 0
ncol = 1
nlay = 30
deltap = ps / nlay  # pressure interval
plev = np.linspace(ps, 0., nlay+1)  # pressure bounds
plev = plev[np.newaxis, ...]
play = np.linspace(ps-deltap/2., deltap/2., nlay)
play = play[np.newaxis, ...]
# Set the temperatures
#  Using a linearly decreasing temperature from surface to TOA
tsfc = 288.
tlay = np.linspace(288.-10., 200., nlay)
tlay = tlay[np.newaxis, ...]
tlev = np.linspace(288., 200., nlay+1)
tlev = tlev[np.newaxis, ...]

# atmospheric composition
# specific humidity profile computed from those temperatures using Manabe's
# fixed relative humidity profile
specific_humidity = np.array([4.13141097e-03, 3.41509495e-03, 2.81099479e-03, 2.30359570e-03,
       1.87920573e-03, 1.52578624e-03, 1.23279279e-03, 9.91026303e-04,
       7.92494475e-04, 6.30283118e-04, 4.98437246e-04, 3.91851633e-04,
       3.06170488e-04, 2.37695932e-04, 1.83304857e-04, 1.40373783e-04,
       1.06711275e-04, 8.04974602e-05, 6.02302082e-05, 4.46774859e-05,
       3.28354282e-05, 2.38916443e-05, 1.71932832e-05, 1.22193649e-05,
       8.55682965e-06, 5.87957411e-06, 5.00000000e-06, 5.00000000e-06,
       5.00000000e-06, 5.00000000e-06])
# Convert to volume mixing ratio from mass mixing ratio
#  just multiplying by ratio of molecular weights of dry air and H2O
h2ovmr = specific_humidity * 28.97 / 18.01528
h2ovmr = h2ovmr[np.newaxis, ...]
#  A global-mean ozone climatology
o3vmr = np.array([2.25573888e-08, 2.38730436e-08, 2.52586476e-08, 2.66442517e-08,
       2.80298557e-08, 2.97254145e-08, 3.14254923e-08, 3.31238355e-08,
       3.46767916e-08, 3.62297478e-08, 3.76122833e-08, 3.86410454e-08,
       3.96698075e-08, 4.08899052e-08, 4.21303310e-08, 4.39781220e-08,
       4.60528063e-08, 4.87636254e-08, 5.16974065e-08, 5.57122567e-08,
       6.17914190e-08, 7.15771368e-08, 9.29020109e-08, 1.29109217e-07,
       1.75914529e-07, 2.45552383e-07, 3.92764464e-07, 7.61726407e-07,
       2.25137178e-06, 7.27500161e-06])
o3vmr = o3vmr[np.newaxis, ...]
# Other values taken from the AquaPlanet Experiment protocols,
# except for O2 which is set the realistic value 0.21
co2vmr = 348. / 1E6 * np.ones_like(play)
ch4vmr = 1650. / 1E9 * np.ones_like(play)
n2ovmr = 306. / 1E9 * np.ones_like(play)
o2vmr = 0.21 * np.ones_like(play)
cfc11vmr = 0. * np.ones_like(play)
cfc12vmr = 0. * np.ones_like(play)
cfc22vmr = 0. * np.ones_like(play)
ccl4vmr = 0. * np.ones_like(play)

nbndlw = int(rrtmg_lw.parrrtm.nbndlw)
ngptlw = int(rrtmg_lw.parrrtm.ngptlw)
#  Initialize absorption data
rrtmg_lw.climlab_rrtmg_lw_ini(cp)

# Lots of RRTMG parameters
icld = 0    # Cloud overlap method, 0: Clear only, 1: Random, 2,  Maximum/random] 3: Maximum
irng = 1  # more monte carlo stuff
idrv = 0  # whether to also calculate the derivative of flux with respect to surface temp
permuteseed = 300
inflglw  = 2
iceflglw = 1
liqflglw = 1

#  These arrays have an extra dimension for number of bands
# in-cloud optical depth [nbndlw,ncol,nlay]
tauc = 0. * np.ones_like(play)
#  broadcast to get [nbndlw,ncol,nlay]
tauc = tauc * np.ones([nbndlw,ncol,nlay])
# Aerosol optical depth at mid-point of LW spectral bands [ncol,nlay,nbndlw]
tauaer = 0. * np.ones_like(play)
#  broadcast and transpose to get [ncol,nlay,nbndlw]
tauaer = np.transpose(tauaer * np.ones([nbndlw,ncol,nlay]), (1,2,0))

# surface emissivity
emis = 1. * np.ones((ncol,nbndlw))

# Clear-sky only
cldfmcl = np.zeros((ngptlw,ncol,nlay))
ciwpmcl = np.zeros((ngptlw,ncol,nlay))
clwpmcl = np.zeros((ngptlw,ncol,nlay))
reicmcl = np.zeros((ncol,nlay))
relqmcl = np.zeros((ncol,nlay))
taucmcl = np.zeros((ngptlw,ncol,nlay))

ispec = 0

Compute clear sky LW flux

# Call the RRTMG_LW driver and output the upwelling LW flux
(olr_sr, uflx, dflx, hr, uflxc, dflxc, hrc, duflx_dt, duflxc_dt) = \
    rrtmg_lw.climlab_rrtmg_lw(ncol, nlay, icld, ispec, idrv,
                 play, plev, tlay, tlev, tsfc,
                 h2ovmr, o3vmr, co2vmr, ch4vmr, n2ovmr, o2vmr,
                 cfc11vmr, cfc12vmr, cfc22vmr, ccl4vmr, emis,
                 inflglw, iceflglw, liqflglw, cldfmcl,
                 taucmcl, ciwpmcl, clwpmcl, reicmcl, relqmcl,
                 tauaer)
print(uflx)

produces

[[390.0993433  375.88912721 365.99472965 357.08147831 348.80377104
  341.07153179 333.97172557 319.82028378 297.58354407 275.39583782
  265.17421804 258.84306736 254.28963431 250.11499435 246.20455494
  242.55139767 239.17104936 236.05032037 233.17039297 230.52254908
  228.0927947  225.87083026 223.827967   221.94192292 220.17919277
  218.53032198 216.9706734  215.40915488 213.85137695 212.02525401
  210.32048901]]

Put in a liquid cloud

icld = 1    # Cloud overlap method, 0: Clear only, 1: Random, 2,  Maximum/random] 3: Maximum

#  Cloud parameters
cloud_level_index = 8
cldfrac = 0.5*np.exp(-(play-play[0,cloud_level_index])**2/(2*25.)**2)  # Layer cloud fraction: a Gaussian centered on a pressure level
clwp = 60. * np.ones_like(play)  # in-cloud liquid water path (g/m2)
ciwp = 0. * np.ones_like(play)   # in-cloud ice water path (g/m2)
relq = 14. * np.ones_like(play) # Cloud water drop effective radius (microns)
reic = 0. * np.ones_like(play)  # Cloud ice particle effective size (microns)


#  Call the Monte Carlo Independent Column Approximation (McICA, Pincus et al., JC, 2003)
(cldfmcl, ciwpmcl, clwpmcl, reicmcl, relqmcl, taucmcl) = \
    rrtmg_lw.climlab_mcica_subcol_lw(
                ncol, nlay, icld,
                permuteseed, irng, play,
                cldfrac, ciwp, clwp, reic, relq, tauc)

# Call the RRTMG_LW driver and output the upwelling LW flux
(olr_sr, uflx, dflx, hr, uflxc, dflxc, hrc, duflx_dt, duflxc_dt) = \
            rrtmg_lw.climlab_rrtmg_lw(ncol, nlay, icld, ispec, idrv,
                 play, plev, tlay, tlev, tsfc,
                 h2ovmr, o3vmr, co2vmr, ch4vmr, n2ovmr, o2vmr,
                 cfc11vmr, cfc12vmr, cfc22vmr, ccl4vmr, emis,
                 inflglw, iceflglw, liqflglw, cldfmcl,
                 taucmcl, ciwpmcl, clwpmcl, reicmcl, relqmcl,
                 tauaer)
print(uflx)

produces

[[390.0993433  375.88912721 365.99472965 357.08147831 348.80377104
  341.07153179 333.9471018  319.68522828 297.36048042 275.07844086
  264.85645432 258.48281265 253.93329287 249.76003965 245.85126144
  242.19996138 238.8215618  235.70281816 232.82486442 230.17895328
  227.75109081 225.53098228 223.48999093 221.60596093 219.84553278
  218.1992751  216.64261065 215.08475758 213.53154191 211.71233183
  210.01478479]]

Change the cloud from water to ice

clwp = 0. * np.ones_like(play)  # in-cloud liquid water path (g/m2)
ciwp = 50 * np.ones_like(play)   # in-cloud ice water path (g/m2)
reic = 15. * np.ones_like(play)  # Cloud ice particle effective size (microns)

#  Call the Monte Carlo Independent Column Approximation (McICA, Pincus et al., JC, 2003)
(cldfmcl, ciwpmcl, clwpmcl, reicmcl, relqmcl, taucmcl) = \
    rrtmg_lw.climlab_mcica_subcol_lw(
                ncol, nlay, icld,
                permuteseed, irng, play,
                cldfrac, ciwp, clwp, reic, relq, tauc)

# Call the RRTMG_LW driver and output the upwelling LW flux
(olr_sr, uflx, dflx, hr, uflxc, dflxc, hrc, duflx_dt, duflxc_dt) = \
            rrtmg_lw.climlab_rrtmg_lw(ncol, nlay, icld, ispec, idrv,
                 play, plev, tlay, tlev, tsfc,
                 h2ovmr, o3vmr, co2vmr, ch4vmr, n2ovmr, o2vmr,
                 cfc11vmr, cfc12vmr, cfc22vmr, ccl4vmr, emis,
                 inflglw, iceflglw, liqflglw, cldfmcl,
                 taucmcl, ciwpmcl, clwpmcl, reicmcl, relqmcl,
                 tauaer)
print(uflx)

produces

[[390.0993433  375.88912721 365.99472965 357.08147831 348.80377104
  341.07153179 333.97172557 319.82028378 297.58354407 275.39583782
  265.17421804 258.84306736 254.28963431 250.11499435 246.20455494
  242.55139767 239.17104936 236.05032037 233.17039297 230.52254908
  228.0927947  225.87083026 223.827967   221.94192292 220.17919277
  218.53032198 216.9706734  215.40915488 213.85137695 212.02525401
  210.32048901]]

[lizhuo1108]

Hi @lizhuo1108, I think your bug results from setting zero effective ice crystal radius reic.

Try changing reic to a realistic value in your code. Or, here is a series of working steps that elaborate on that example:

Set up a single column

import numpy as np
from climlab_rrtmg import rrtmg_lw, rrtmg_sw

# Specific heat at constant pressure
cp = 1004.

# Set up the pressure domain
ps = 1000. #  Surface pressure in hPa
# RRTM code expects arrays with (ncol, nlay)
# and with pressure decreasing from surface at element 0
ncol = 1
nlay = 30
deltap = ps / nlay  # pressure interval
plev = np.linspace(ps, 0., nlay+1)  # pressure bounds
plev = plev[np.newaxis, ...]
play = np.linspace(ps-deltap/2., deltap/2., nlay)
play = play[np.newaxis, ...]
# Set the temperatures
#  Using a linearly decreasing temperature from surface to TOA
tsfc = 288.
tlay = np.linspace(288.-10., 200., nlay)
tlay = tlay[np.newaxis, ...]
tlev = np.linspace(288., 200., nlay+1)
tlev = tlev[np.newaxis, ...]

# atmospheric composition
# specific humidity profile computed from those temperatures using Manabe's
# fixed relative humidity profile
specific_humidity = np.array([4.13141097e-03, 3.41509495e-03, 2.81099479e-03, 2.30359570e-03,
       1.87920573e-03, 1.52578624e-03, 1.23279279e-03, 9.91026303e-04,
       7.92494475e-04, 6.30283118e-04, 4.98437246e-04, 3.91851633e-04,
       3.06170488e-04, 2.37695932e-04, 1.83304857e-04, 1.40373783e-04,
       1.06711275e-04, 8.04974602e-05, 6.02302082e-05, 4.46774859e-05,
       3.28354282e-05, 2.38916443e-05, 1.71932832e-05, 1.22193649e-05,
       8.55682965e-06, 5.87957411e-06, 5.00000000e-06, 5.00000000e-06,
       5.00000000e-06, 5.00000000e-06])
# Convert to volume mixing ratio from mass mixing ratio
#  just multiplying by ratio of molecular weights of dry air and H2O
h2ovmr = specific_humidity * 28.97 / 18.01528
h2ovmr = h2ovmr[np.newaxis, ...]
#  A global-mean ozone climatology
o3vmr = np.array([2.25573888e-08, 2.38730436e-08, 2.52586476e-08, 2.66442517e-08,
       2.80298557e-08, 2.97254145e-08, 3.14254923e-08, 3.31238355e-08,
       3.46767916e-08, 3.62297478e-08, 3.76122833e-08, 3.86410454e-08,
       3.96698075e-08, 4.08899052e-08, 4.21303310e-08, 4.39781220e-08,
       4.60528063e-08, 4.87636254e-08, 5.16974065e-08, 5.57122567e-08,
       6.17914190e-08, 7.15771368e-08, 9.29020109e-08, 1.29109217e-07,
       1.75914529e-07, 2.45552383e-07, 3.92764464e-07, 7.61726407e-07,
       2.25137178e-06, 7.27500161e-06])
o3vmr = o3vmr[np.newaxis, ...]
# Other values taken from the AquaPlanet Experiment protocols,
# except for O2 which is set the realistic value 0.21
co2vmr = 348. / 1E6 * np.ones_like(play)
ch4vmr = 1650. / 1E9 * np.ones_like(play)
n2ovmr = 306. / 1E9 * np.ones_like(play)
o2vmr = 0.21 * np.ones_like(play)
cfc11vmr = 0. * np.ones_like(play)
cfc12vmr = 0. * np.ones_like(play)
cfc22vmr = 0. * np.ones_like(play)
ccl4vmr = 0. * np.ones_like(play)

nbndlw = int(rrtmg_lw.parrrtm.nbndlw)
ngptlw = int(rrtmg_lw.parrrtm.ngptlw)
#  Initialize absorption data
rrtmg_lw.climlab_rrtmg_lw_ini(cp)

# Lots of RRTMG parameters
icld = 0    # Cloud overlap method, 0: Clear only, 1: Random, 2,  Maximum/random] 3: Maximum
irng = 1  # more monte carlo stuff
idrv = 0  # whether to also calculate the derivative of flux with respect to surface temp
permuteseed = 300
inflglw  = 2
iceflglw = 1
liqflglw = 1

#  These arrays have an extra dimension for number of bands
# in-cloud optical depth [nbndlw,ncol,nlay]
tauc = 0. * np.ones_like(play)
#  broadcast to get [nbndlw,ncol,nlay]
tauc = tauc * np.ones([nbndlw,ncol,nlay])
# Aerosol optical depth at mid-point of LW spectral bands [ncol,nlay,nbndlw]
tauaer = 0. * np.ones_like(play)
#  broadcast and transpose to get [ncol,nlay,nbndlw]
tauaer = np.transpose(tauaer * np.ones([nbndlw,ncol,nlay]), (1,2,0))

# surface emissivity
emis = 1. * np.ones((ncol,nbndlw))

# Clear-sky only
cldfmcl = np.zeros((ngptlw,ncol,nlay))
ciwpmcl = np.zeros((ngptlw,ncol,nlay))
clwpmcl = np.zeros((ngptlw,ncol,nlay))
reicmcl = np.zeros((ncol,nlay))
relqmcl = np.zeros((ncol,nlay))
taucmcl = np.zeros((ngptlw,ncol,nlay))

ispec = 0

Compute clear sky LW flux

# Call the RRTMG_LW driver and output the upwelling LW flux
(olr_sr, uflx, dflx, hr, uflxc, dflxc, hrc, duflx_dt, duflxc_dt) = \
    rrtmg_lw.climlab_rrtmg_lw(ncol, nlay, icld, ispec, idrv,
                 play, plev, tlay, tlev, tsfc,
                 h2ovmr, o3vmr, co2vmr, ch4vmr, n2ovmr, o2vmr,
                 cfc11vmr, cfc12vmr, cfc22vmr, ccl4vmr, emis,
                 inflglw, iceflglw, liqflglw, cldfmcl,
                 taucmcl, ciwpmcl, clwpmcl, reicmcl, relqmcl,
                 tauaer)
print(uflx)

produces

[[390.0993433  375.88912721 365.99472965 357.08147831 348.80377104
  341.07153179 333.97172557 319.82028378 297.58354407 275.39583782
  265.17421804 258.84306736 254.28963431 250.11499435 246.20455494
  242.55139767 239.17104936 236.05032037 233.17039297 230.52254908
  228.0927947  225.87083026 223.827967   221.94192292 220.17919277
  218.53032198 216.9706734  215.40915488 213.85137695 212.02525401
  210.32048901]]

Put in a liquid cloud

icld = 1    # Cloud overlap method, 0: Clear only, 1: Random, 2,  Maximum/random] 3: Maximum

#  Cloud parameters
cloud_level_index = 8
cldfrac = 0.5*np.exp(-(play-play[0,cloud_level_index])**2/(2*25.)**2)  # Layer cloud fraction: a Gaussian centered on a pressure level
clwp = 60. * np.ones_like(play)  # in-cloud liquid water path (g/m2)
ciwp = 0. * np.ones_like(play)   # in-cloud ice water path (g/m2)
relq = 14. * np.ones_like(play) # Cloud water drop effective radius (microns)
reic = 0. * np.ones_like(play)  # Cloud ice particle effective size (microns)


#  Call the Monte Carlo Independent Column Approximation (McICA, Pincus et al., JC, 2003)
(cldfmcl, ciwpmcl, clwpmcl, reicmcl, relqmcl, taucmcl) = \
    rrtmg_lw.climlab_mcica_subcol_lw(
                ncol, nlay, icld,
                permuteseed, irng, play,
                cldfrac, ciwp, clwp, reic, relq, tauc)

# Call the RRTMG_LW driver and output the upwelling LW flux
(olr_sr, uflx, dflx, hr, uflxc, dflxc, hrc, duflx_dt, duflxc_dt) = \
            rrtmg_lw.climlab_rrtmg_lw(ncol, nlay, icld, ispec, idrv,
                 play, plev, tlay, tlev, tsfc,
                 h2ovmr, o3vmr, co2vmr, ch4vmr, n2ovmr, o2vmr,
                 cfc11vmr, cfc12vmr, cfc22vmr, ccl4vmr, emis,
                 inflglw, iceflglw, liqflglw, cldfmcl,
                 taucmcl, ciwpmcl, clwpmcl, reicmcl, relqmcl,
                 tauaer)
print(uflx)

produces

[[390.0993433  375.88912721 365.99472965 357.08147831 348.80377104
  341.07153179 333.9471018  319.68522828 297.36048042 275.07844086
  264.85645432 258.48281265 253.93329287 249.76003965 245.85126144
  242.19996138 238.8215618  235.70281816 232.82486442 230.17895328
  227.75109081 225.53098228 223.48999093 221.60596093 219.84553278
  218.1992751  216.64261065 215.08475758 213.53154191 211.71233183
  210.01478479]]

Change the cloud from water to ice

clwp = 0. * np.ones_like(play)  # in-cloud liquid water path (g/m2)
ciwp = 50 * np.ones_like(play)   # in-cloud ice water path (g/m2)
reic = 15. * np.ones_like(play)  # Cloud ice particle effective size (microns)

#  Call the Monte Carlo Independent Column Approximation (McICA, Pincus et al., JC, 2003)
(cldfmcl, ciwpmcl, clwpmcl, reicmcl, relqmcl, taucmcl) = \
    rrtmg_lw.climlab_mcica_subcol_lw(
                ncol, nlay, icld,
                permuteseed, irng, play,
                cldfrac, ciwp, clwp, reic, relq, tauc)

# Call the RRTMG_LW driver and output the upwelling LW flux
(olr_sr, uflx, dflx, hr, uflxc, dflxc, hrc, duflx_dt, duflxc_dt) = \
            rrtmg_lw.climlab_rrtmg_lw(ncol, nlay, icld, ispec, idrv,
                 play, plev, tlay, tlev, tsfc,
                 h2ovmr, o3vmr, co2vmr, ch4vmr, n2ovmr, o2vmr,
                 cfc11vmr, cfc12vmr, cfc22vmr, ccl4vmr, emis,
                 inflglw, iceflglw, liqflglw, cldfmcl,
                 taucmcl, ciwpmcl, clwpmcl, reicmcl, relqmcl,
                 tauaer)
print(uflx)

produces

[[390.0993433  375.88912721 365.99472965 357.08147831 348.80377104
  341.07153179 333.97172557 319.82028378 297.58354407 275.39583782
  265.17421804 258.84306736 254.28963431 250.11499435 246.20455494
  242.55139767 239.17104936 236.05032037 233.17039297 230.52254908
  228.0927947  225.87083026 223.827967   221.94192292 220.17919277
  218.53032198 216.9706734  215.40915488 213.85137695 212.02525401
  210.32048901]]

Hi Brian,
Thanks for the info! I output ciwp and reic from CAM3. Maybe there exist some difference (maybe unit) in the cloud information needed for CAM3 radiation scheme (CAMRT) and RRTMG. I will take a further look at the cloud information.

[lizhuo1108]

Hi All,

I checked the code of RRTMG and I found the reason why it crashed when the effective radius becomes very large. For the source code:
climlab_rrtmg/rrtmg_sw/rrtmg_sw_v4.0/gcm_model/src/rrtmg_sw_cldprmc.f90

Line 182-183:
elseif (iceflag .eq. 1) then
if (radice .lt. 13.0_rb .or. radice .gt. 130._rb) stop &

When the iceflag (iceflglw) is set to 1 and the radius is unreasonably large (> 130) or small (<13), the RRTMG will be forced to stop and crash. A remind for other users: when you give the information of the effective ice radius (radice) to RRTMG, please check your input radice and make sure that the radice is within this range, otherwise RRTMG will crash! Same is for other two options (iceflag=2 and 3).

Li

[brian-rose]

When the iceflag (iceflglw) is set to 1 and the radius is unreasonably large (> 130) or small (<13), the RRTMG will be forced to stop and crash. A remind for other users: when you give the information of the effective ice radius (radice) to RRTMG, please check your input radice and make sure that the radice is within this range, otherwise RRTMG will crash! Same is for other two options (iceflag=2 and 3).

This is useful info! This package doesn't have any actual "documentation" but it might be helpful at least to put this into comments in the Python wrapper code.