DOI: 10.5061/dryad.w9ghx3fxz
title: Cloud_ICA: A deterministic cloud-overlap algorithm for generating a complete set of independent column atmospheres

description: In calculating solar radiation, climate models make many simplifications,
 in part to reduce computational cost and enable climate modeling, and in
 part from lack of understanding of critical atmospheric information.
 Whether known errors or unknown errors, the community's concern is
 how these could impact the modeled climate. The simplifications are well
 known and most have published studies evaluating them, but with individual
 studies it is difficult to compare. Here, we collect a wide range of such
 simplifications in either radiative transfer modeling or atmospheric
 conditions and assess potential errors within a consistent framework on
 climate‐relevant scales. We build benchmarking capability around a solar
 heating code (Solar‐J) that doubles as a photolysis code for chemistry and
 can be readily adapted to consider other errors and uncertainties. The
 broad classes here include: use of broad wavelength bands to integrate
 over spectral features; scattering approximations that alter phase
 function and optical depths for clouds and gases; uncertainty in ice‐cloud
 optics; treatment of fractional cloud cover including overlap; and
 variability of ocean surface albedo. We geographically map the errors in W
 m−2 using a full climate re‐creation for January 2015 from a weather
 forecasting model. For many approximations assessed here, mean errors are
 ∼2 W m−2 with greater latitudinal biases and are likely to affect a
 model's ability to match the current climate state. Combining this
 work with previous studies, we make priority recommendations for fixing
 these simplifications based on both the magnitude of error and the ease or
 computational cost of the fix.
description: A unique and useful advantage of Cloud-ICA is that with low
 computational cost it readily generates a complete set of ICAs such that
 the sum of wtICAs equals one.  Cloud-ICA can then calculate a set of
 up-to-four quadrature column atmospheres (QCAs) that can be used to
 approximate the integral over the full set of ICAs.  The method of parsing
 partially overlapping clouds within a single column atmosphere (SCA) into
 a large number of ICAs and then into a set of only four QCAs originated in
 Neu, Prather & Penner (2007) and was fully developed with
 observation-based overlap models in the chemistry photolysis code Cloud-J
 v7.3c in Prather (2015).  Cloud-J has continued to develop (Prather
 & Hsu, 2019;  2019d; Hsu & Prather, 2021) and the current
 Cloud-ICA code is taken from Cloud-J v8.0c (Prather, 2023). 
 <strong>Note added
 (2024/05/20): </strong> The program as written limits the number of
 ICAs to the dimension size ICA_ = 20,000 = 2 x
 10<sup>4</sup>.  If using the preferred 6 max-overlap vertical
 blocks plus cirrus shield (LNRG=6) and the quantization of cloud fraction
 in units of 1/10, then the maximum number of ICAs is limited to about 5 x
 10<sup>6</sup>.   While this maximum number has not been found
 by us with realistic atmospheres, it is possible. We did not want to
 dimension and scan/sort such large numbers of ICAs when running multi-year
 CTM simulations and thus truncated it at 2 x 10<sup>4</sup>. 
 This truncation is a soft landing and just combines upper level
 decorrelated blocks into a single max-overlap block until the ICA_ limit
 is met.  The invocation of this limit does invoke a 'write'
 statement that may disrupt parallel optimization, and it can be commented
 away. NICA greater
 than ICA_       22680       20000        7560           7          
 6 Beware that if you run LNRG=0, then alternating
 cloudy and clear layers can generate 2<sup>CLD_/2</sup>
 ICAs. The primary Cloud-ICA references
 are:   Neu, J.L., M.J.
 Prather, J.E. Penner (2007) Global atmospheric chemistry:  integrating
 over fractional cloud cover, J. Geophys. Res., 112, D11306,
 doi:10.1029/2006JD008007  
 Prather, M.J. (2015) Photolysis rates in correlated overlapping
 cloud fields: Cloud-J 7.3c, Geosci. Model Dev., 8, 2587-2595,
 doi:10.5194/gmd-8-2587-2015  
 Prather, M.J. and J.C. Hsu (2019) A round Earth for climate
 models, Proc Natl Acad Sci, 116 (39) 19330-19335;
 https://doi.org/10.1073/pnas.1908198116.
   Prather, Michael and Hsu, Juno
 (2019d), Solar-J and Cloud-J models version 7.6c, v2, UC Irvine Dash,
 Dataset, doi.org/10.7280/D1096P  
 Hsu, J.C. and M.J. Prather (2021), Assessing uncertainties and
 approximations in solar heating of the climate system. Journal of Advances
 in Modeling Earth Systems, 13, e2020MS002131. doi:
 10.1029/2020MS002131  
 Prather, Michael (2023d). An updated cloud-overlap photolysis
 module for atmospheric chemistry models, UCI Cloud-J v8.0, with near-UV
 H2O absorption [Dataset]. Dryad.
 https://doi.org/10.7280/D1Q398  
 A description of the cloud overlap algorithm and it basis in
 observations is take directly from Prather (2015), see references
 within.   “<em>Our recommended
 cloud-overlap model uses the information on vertical correlations (Pincus
 et al., 2005; Naud and DelGenio, 2006; Kato et al., 2010; Oreopoulis et
 al., 2012), which shows cloud decorrelation lengths on the order of 1.5 km
 in the lower atmosphere increasing to 3 km or more in the upper
 troposphere. Since a true COR model scales as 2**NL and becomes rapidly
 impractical for high-resolution models, we define vertical groups of cloud
 layers globally according to the decorrelation lengths: 0–1.5 km altitude,
 1.5–3.5, 3.5–6, 6–9, 9–13, and &gt;13 km. We assume that the cloud
 layers within a decorrelation length are highly correlated with one
 another and thus form a MAX group. When such MAX groups are adjacent they
 have a mean separation of one decorrelation length, and we choose a cloud
 correlation factor of cc=0.33, similar to 1 e-fold. When there is a
 clear-sky gap between a pair of G6 layers, the MAX groups are separated by
 more than one decorrelation length; thus, we reduce the factor cc with
 successive multiples (i.e., with two missing G6 MAX groups between two
 cloudy layers,the effective cc=0.333**2 =0.036). This model is denoted
 G6/.33. Two other G6 models were tested: cc=0.00 corresponds to randomly
 overlapped adjacent groups (MAXRAN, G6/.00); and cc=0.99 is almost
 maximally overlapped (MAX, G6/.99).</em>   I<em>n looking at how this
 model aligned the clouds for realistic FCAs, we found that extensive
 cirrus fractions in the uppermost layers prevented the expected overlap of
 small fraction cumulus below. Thus, a seventh MAX group is added if there
 was a cirrus shield (defined from top down as adjacent ice-only clouds
 with f &gt;0.5). Because of the cloud fraction binning into 10%
 intervals, the number of ICAs is bounded by 5x106 (including the cirrus
 shield). This limit is resolution independent and was never reached in any
 FCAs examined here (highest number of ICAs for one FCA was 3500). The
 major computational cost comes with the Fast-J computation, and the
 methods for approximating the average of J values over all ICAs (Sect. 3)
 use at most four Fast-J calculations no matter how many ICAs.  Two other
 cloud-overlap models tested here are the MAXRAN groupings G0 and G3 (Feng
 et al., 2004; Neu et al., 2007). Model G0 assumes that all vertically
 adjacent cloudy layers are a MAX group (maximally overlapped), and all
 such groups separated by a clear layer are RAN overlapped.  This model
 seems logical but has difficulty finding a clear layer when the FCA has
 been averaged over several hours or taken from a parameterized
 cloud-resolving model.  It our tests, using meteorological data with
 NLD36, the maximum number of G0 ICAs was 375. Model G3 has at most three
 MAX-RAN groups demarcated by atmospheric regimes: a fixed altitude (1.5
 km, stratus top) and temperature (the liquid-to-ice cloud transition). The
 maximum possible number of ICAs per FCA for G3 is 103, and in our tests we
 found 288.</em>   <em>Our recommended
 cloud-overlap model is G6/.33 since it is based on the observed–modeled
 cloud decorrelation lengths. For a given FCA, we treat the J values
 calculated by summing Fast-J over all the ICAs generated by G6/.33 as the
 correct value.  We calculate errors for the other cloud-overlap models
 (here) or various ICA-approximation models using the G6/.33 model (Sect.
 3).</em>   <em>We use a high-resolution
 snapshot from the European Center for Medium-rangeWeather Forecasts,
 similar to what is used (at lower resolution) in the UC Irvine and
 University of Oslo chemistry-transport models (Søvde et al., 2012; Hsu and
 Prather, 2014). The 640 FCAs are a 3 h average of a single longitudinal
 belt just above the Equator (T319L60 Cycle 36) and have clouds only in the
 lowermost 36 layers. Profiles of temperature and ozone are taken from
 tropical mean observations; the Rayleigh-scattering optical depth at 600
 nm is about 0.12, and a mix of aerosol layers has a total optical depth of
 0.23. J value errors are calculated separately for each FCA and then
 averaged. The number of ICAs per FCA averages 169 for model G6, 21 for
 model G3, and 19 for model G0; see Fig. 2 for the probability distribution
 of ICA numbers.</em>   Further
 discussion about the deterministic cloud cover generator is found in Hsu
 and Prather (2021)   <em>“In a manner similar to
 Hogan and Bozzo's (2018) deterministic cloud-cover generator that
 goes from MAXRAN to EXP-RAN, Cloud-J developed a deterministic ICA
 generator for MAX-RAN and then adapted it to use vertical decorrelation
 lengths in its MAX-COR algorithm (M. J. Prather, 2015). Chemistry models
 need the selection of ICAs for any overlap method to be deterministic
 because many critical applications require perturbation-control pairs
 without stochastic noise (e.g., M. J. Prather & Hsu, 2010). Thus
 Solar-J cannot use a stochastic cloud generator (e.g., Räisänen et al.,
 2004), and this drove the structure of our cloud overlap algorithm.
 MAX-COR was designed to be (i) deterministic, (ii) linear in cost with
 increasing numbers of layers, and (iii) robust when cloud data are
 averaged in time or space, because such averaging tends to eliminate
 cloud-free layers and revert to MAX overlap. Based on observations of
 decorrelation length (Kato et al., 2010; Naud et al., 2008; Oreopoulos et
 al., 2012; Pincus et al., 2005), MAX-COR defines 6-layer groupings by
 altitude range. Because decorrelation is small across the vertical range
 of each group, we assume MAX overlap within each group and a decorrelation
 of the overlap of each MAX group with its neighbor. Adopting terminology
 of climate community, MAX-COR is effectively a MAX-EXP algorithm. By
 quantizing the cloud fraction to the nearest 10% and allowing an
 independent cirrus shield at the top, the absolute maximum number of ICAs
 under MAX-EXP is &lt;5 × 106 and thus ICAs can be rapidly defined and
 binned with low computational overhead. Deterministic EXP-EXP or EXP-RAN
 models in our code would have to enumerate up to 233 ICAs for our model
 that has potentially 33 cloudy layers, which is truly prohibitive and not
 linearly scalable with resolution. We believe that a MAX-COR or MAX-EXP
 algorithm is likely the most stable and scalable deterministic ICA
 generator for vertical cloud decorrelation algorithms. The RRTMG v4.0 code
 available at the time of this study uses primarily MAX-RAN cloud overlap,
 but the new v5.0 code includes an EXP-RAN option. Thus, our comparisons of
 cloud-overlap results with the RRTMG code are limited to MAX-RAN. Within
 Solar-J we can run both MAX-RAN (SJ/RAN) and the standard MAX-COR (SJ) and
 thus compare with J. K. P. Shonk and Hogan (2010), as discussed
 below.</em>   <em>Let us accept that ICAs
 generated by cloud overlap algorithms can be solved with 1D RT as
 horizontally homogeneous plane parallel layers, then the next step is how
 to solve the RT problem for all ICAs and average the results. The number
 of ICAs are often numerous enough that no practical climate RT code can
 solve them all, and most codes do not even count them all (Räisänen et
 al., 2004). RRTMG randomly selects an ICA for each wavelength bin in the
 RT solution, a method designated Monte Carlo ICA McICA, Pincus et al.,
 2003). McICA has errors at each time step by mixing ICAs across
 wavelengths and by not accurately sampling the average of ICAs (e.g.,
 average cloud optical depth) in that time step. McICA is intended to
 deliver the correct mean when averaged enough times over the same cloud
 system, but it has hourly grid-cell rms errors of 40 W m−2 (H. W. Barker
 et al., 2008; Pincus et al., 2003). A key underlying premise is that solar
 heating errors propagate symmetrically and linearly in the climate system
 and average out, as was found for simple forecast models. Assessing net
 bias errors caused by noisy heating rates would need to examine nonlinear
 processes in hydrology, cloud systems, ecosystem productivity, and air
 quality in Earth system models (e.g., Pincus & Stevens,
 2013).</em>   <em>With a deterministic ICA
 generator, we can calculate an "exact" non-stochastic answer as
 was done for limited test cases in M. J. Prather (2015), but we could not
 afford to do this for our January climate metric.  Solar-J identifies and
 sorts all ICAs by cloud optical depth and then selects up to four
 representative quadrature column atmospheres (QCAs) each with a fractional
 area to represent the distribution of ICAs. The full-wavelength RT
 solutions are completed for each QCA (Neu et al., 2007). See Figure S3 for
 a global picture of the average frequency of occurrence of the 4 QCA bins
 for January 2015. Cloud quadrature does a very good job of averaging over
 the ICAs with net bias errors of ∼1% in solar intensity and rms errors of
 2%–4%. To reach equivalent accuracy for a single time step using random
 selection would require about 50 ICAs each with full wavelength
 calculation (not as in McICA) versus an average of 2.8 QCAs (many grid
 cells have less than 4 QCAs).</em>
description: # Assessing Uncertainties and Approximations in Solar Heating of the
 Climate System
 [https://doi.org/10.5061/dryad.w9ghx3fxz](https://doi.org/10.5061/dryad.w9ghx3fxz) The primary goal of the computer code published here is to provide an algorithm for sorting an atmosphere of overlapping clouds into a complete and finite set of independent column atmospheres (ICAs), wherein each layer is horizontally homogenous (i.e., plane parallel) and either uniformly cloudy or clear.  Many current solar heating codes can only estimate the distribution of ICAs with Monte Carlo sampling.  This algorithm gets around that problem.  It has been in use with UCI photochemistry codes, but it was so deeply embedded as to make it unusable by others.  Here we extract that algorithm to enable ready use in other Earth system modeling codes.  This work presents a stand-alone FORTRAN 90 computer code Cloud-ICA that reads in a single column atmosphere (p at edges, T in layer) where each layer includes a cloud fraction (CldF = 0-1) plus the liquid- and/or ice-water path of the cloud in the layer (not the average over the clear and cloudy layer).  The cloud data are used to calculate the visible optical depth of the cloud itself (CldTau), not the average optical depth of the layer.  With these two quantities and an overlap algorithm, the code calculates a subset of independent column atmospheres (ICAs), each occupying a fraction of the original cell (wtICA) and each with layers that are either 100% cloudy or clear, thus enabling plane-parallel radiative transfer solutions.  ## Description of the data and file structure There are two sample output files: **ICA_print1.out** = detailed printout for the first cloudy atmosphere that has only 36 ICAs.  This gives the values for the grouping indices and the profiles of the 36 ICAs in terms of cloud visible optical depth in leach model layer. See the rest of this readme file for a description of the variables. **ICA_nica640.out** = a quick scan through 640 different tropical atmospherrs listing the number of ICAs in each (NICA) the resulting QCAs in terms of weights and column optical depth. **atmos_PTClds.dat** = gives a standard atmospheric profile (p, T) and then 640 different profiles of cloud systems, sampled from ECMWF fields in the tropics. Variables in the printout examples **&gt;&gt;&gt;&gt;&gt;atmospheric data** L = model layer, must be from surface(L=1) to top-of-atmosphere (L=58)in this case.  The inverse layering notation (L=1=toa)is not supported and may not work.  For physics, in our model L=1:LATM\_ For clouds, it is limited to L=1:CLD\_ (&gt;&gt;&gt;&gt;cloud overlap groups and quantization** FBIN (integer) number of quantized bins allowed for CLDF.  Example, if FBIN = 10, then CLDF is quantized into fixed fractions: 0.0, 0.1, 0.2, … 0.9, 1.0.  Note that if CLF &lt; CLF_MIN = 0.02) then CLDF = 0.0 and the layer is cleared of clouds!  The optical depth of the cloud (CLTAU) is scaled to conserve CLF*CLTAU.  Quantized CLDF means any max-overlap group has at most 10 possible combinations. NCLDF(1:LCLD_) = cloud fraction quantum number (1:FBIN) for layer. NCLDF = 0 means NO clouds LNRG (= 0, 3, 6) index for Cloud-J to establish the number of semi-independent maximally overlapped cloud groups in an atmospheric column. LNRG=0.  Defines a break in the max-overlap cloud groups by the cloud optical depth quantum. NG_BRK=0 is default = clear layer, see FBIN and NCLDF below. LNRG=3:  Defines 3 separate max-overlap groups strictly by altitude, that is presently locked to our ECMWF L57 grid:  GRP=1 is L=1:8 (surf to +1 km); GRP=2 is L=9 to uppermost liquid water cloud; GRP=3 is for all-ice clouds.  Note that mixed-phase clouds count as liquid water since they dominate the liquid water cloud opacity.  These groups may be merged if they are without clouds. LNRG=6:  Current recommendation is for correlated groups of max-overlap clouds.  Uses altitude (actually calculated for the true model grid, not presumed to be ECMWF L57) to create 6 max-overlap groups, based on observed correlation lengths (see Cloud-J papers).  GRP#1 is 0-1.5km, GRP#2 is 1.5-3.5km,  GRP#3 is 3.5-6km, GRP#4 is 6-9km, GRP#5 is 9-13km, and GRP#6 is 13km+. LCIRRUS = layer number of the top cirrus shield group (if there is one).  New Group:  if there is a CLF~1 cirrus shield, pull it off the top of the maximum GRP, and create a separate GRP#7. (Could be lower GRP# is some of GRP#1-6 have no clouds. NRG (integer) final number of maximally overlapped groups, should be 1 to 7 in current code depending on numer of groups w/o clouds.  Code is dimensioned to accommodate 9 total groups. GMIN(1:NRG), GMAX(1:NRG) = layer number of the bottom and top of each max-ran group.  All possible cloud layers (1:LCLD_) must be in a group.  The numbers start with GMIN(1)=1 and end with GMAX(NRG) = LCLD_.  GNR(1:NRG) = no. of unique quantized cloud fractions in group, .le.FBIN+1 GBOT(1:NRG) = lowest model layer of max-overlap group (similar to GMIN) GTOP(1:NRG) = uppermost model layer of max-overlap group (similar to GMAX) GCMX(1:NRG) = maximum cloud fraction (quantum # 1:FBIN) of the max-overlap group GLVL(1:NRG) = G6 group # (1:7,for LNRG=6 only) of the NRG group #.  This keeps track of a missing (uncloudy) G6 group so that the decorrelation factor is calculated correctly, with extra factor for missing groups.  Note that there can be 7 G6 groups if the upper cirrus shield is separated so as to be decorrelated with clouds below. GFNR(1:NRG,1:GNR) = GFNR(N,G) = cloud fraction quantum no. (1 to CBIN_) for max-overlap group N’s cloud fraction # G. **&gt;&gt;&gt;&gt;&gt;independent column atmospheres** ICA\_ = dimensioned size for ICAs, currently 20,000.  Can be larger.  With LNRG=6 and FBIN=10, the maximum possible is about 5e6, but for LNRG=0, the maximum could be 2**(LCLD_/2) (i.e., alternating cloudy and clear layers). NICA = total number of ICAs WCOL(1:NICA) = fractional area weight of each ICA OCOL(1:NICA) = column visible optical depth of each ICA OPRF(1:LCLD_,1:NICA) = visible optical depth in each layer of each ICA ISORT(1:NICA) = number of ICA (original sequence generated by groups) when sorted in order from smallest OCOL to largest OCOL. **&gt;&gt;&gt;&gt;&gt;quadrature column atmospheres** N.B. These QCA variables are not returned in the main subroutine CLOUD-ICA because QCAs may not necessarily be invoked (need CLDFLAG = 7 or 8). In Cloud_j, the QCA calculation is called internally.  These can be extracted if the QCAs are desired. WTQCA(1:NQD_) = fractional area weight of QCA = 1:4, may be 0. ODQCA(1:NQD_) = column cloud visible optical depth of QCA QCOLN(1:LCLD_,1:NDQ_) = visible optical depth in each layer of each ICA NQ1(1:NQD_) = starting ICA number in the QCA group (based on ISORT, sorted ICAs!) NQ2(1:NQD_) = ending ICA number in the QCA group (from ISORT) ## Code/Software The core computer program and standalone driver is archived in the accompanying zenodo files.  These are FORTRAN 90 files, written on a pc (cr/lf).