IFS DOCUMENTATION – Cy48r1

cloud sub-grid structure using the Monte Carlo Independent Column Approximation (McICA, Pincus et al., 2003). To reduce the overall computational cost of the radiation scheme, full radiation calculations are performed on a reduced (coarser) radiation grid and/or on a reduced time frequency. The results are then interpolated back to the original grid. The turbulent diffusion scheme (Chapter 3) represents the vertical exchange of heat, momentum and moisture through sub-grid scale turbulence. The vertical turbulent transport is treated differently in the surface layer and above. In the surface layer, the turbulence fluxes are computed using a first order K-diffusion closure based on the Monin–Obukhov (MO) similarity theory. Above the surface layer a K- diffusion turbulence closure is used everywhere, except for unstable boundary layers where an Eddy- Diffusivity Mass-Flux (EDMF) framework is applied, to represent the non-local boundary layer eddy fluxes (K ¨ohler et al., 2011). The scheme is written in moist conserved variables (liquid static energy and total water). Convective clouds are treated separately by the shallow convection scheme. The effects of unresolved orography on the atmospheric flow are parametrized as a sink of momentum (drag). The turbulent diffusion scheme includes a parametrization in the lower atmosphere to represent the turbulent orographic form drag induced by small scale (< 5 km) orography (Beljaars et al., 2004b). In addition, in stably stratified flow, the orographic drag parametrization (Chapter 4) represents the effects of low-level blocking due to unresolved orography (blocked flow drag) and the absorption and/or reflection of vertically propagating gravity waves (gravity wave drag) on the momentum budget (Lott and Miller, 1997). The non-orographic gravity wave drag parametrization (Chapter 5) accounts for the effects of unresolved non-orographic gravity waves. These waves are generated in nature by processes like deep convection, frontal disturbances, and shear zones. Propagating upward from the troposphere the waves break in the middle atmosphere, comprising the stratosphere and the mesosphere, where they exert a strong drag on the flow. The parametrization uses a globally uniform wave spectrum, and propagates it vertically through changing horizontal winds and air density, thereby representing the wave breaking effects due to critical level filtering and non-linear dissipation. A description of the scheme and its effects on the middle atmosphere circulation can be found in Orr et al. (2010). The moist convection scheme (Chapter 6) is based on the mass-flux approach and represents deep (including congestus), shallow and mid-level (elevated moist layers) convection. The distinction between deep and shallow convection is made on the the basis of the cloud depth (≤ 200 hPa for shallow). For deep convection the mass-flux is determined by assuming that convection removes Convective Available Potential Energy (CAPE) over a given time scale. The intensity of shallow convection is based on the budget of the moist static energy, i.e. the convective flux at cloud base equals the contribution of all other physical processes when integrated over the subcloud layer. Finally, mid- level convection can occur for elevated moist layers, and its mass flux is set according to the large-scale vertical velocity. The scheme, originally described in Tiedtke (1989), has evolved over time and amongst many changes includes a modified entrainment formulation leading to an improved representation of tropical variability of convection (Bechtold et al., 2008), a modified CAPE closure leading to an improved diurnal cycle of convection (Bechtold et al., 2014) and an additional closure dependence on total moisture convergence that improves mesoscale organisation (Becker et al., 2021). Clouds and large-scale precipitation (Chapter 7) are parametrized with a number of prognostic equations for cloud liquid, cloud ice, rain and snow water contents and a sub-grid fractional cloud cover. The cloud scheme represents the sources and sinks of cloud and precipitation due to the major generation and destruction processes, including cloud formation by detrainment from cumulus convection, condensation, deposition, evaporation, collection, melting and freezing. The scheme is based on (Tiedtke, 1993) but with an enhanced representation of mixed-phase clouds and prognostic precipitation (Forbes and Tompkins, 2011; Forbes et al., 2011). Supersaturation with respect to ice is commonly observed in the upper troposphere and is also represented in the parametrization (Tompkins et al., 2007). The surface parametrization scheme (Chapter 8) represents the surface fluxes of energy and water and corresponding sub-surface quantities. The scheme is based on a tiled approach (TESSEL) representing different sub-grid surface types for vegetation, bare soil, snow and open water. Each tile has its own 6 IFS Documentation – Cy48r1 Part IV: Physical Processes properties defining separate heat and water fluxes used in an energy balance equation which is solved for the tile skin temperature. Four soil layers are represented as well as snow mass and density. The evaporative fluxes consider separately the fractional contributions from snow cover, wet and dry vegetation and bare soil. An interception layer collects water from precipitation and dew fall, and infiltration and run-off are represented depending on soil texture and subgrid orography (HTESSEL, Balsamo et al., 2009, Balsamo et al., 2011a). A carbon cycle is included and land-atmosphere exchanges of carbon dioxide are parametrized to respond to diurnal and synoptic variations in the water and energy cycles (CHTESSEL, Boussetta et al., 2013a). Chapter 9 ‘Methane oxidation’ describes a simple parametrization of the upper-stratospheric moisture source due to methane oxidation. A parametrization representing photolysis of vapour in the mesosphere is also included. Chapter 10 ‘Ozone chemistry parametrization’ gives a brief description of the Hybrid Linear Ozone (HLO) scheme, which assumes that chemical changes in ozone can be described by a linear relaxation towards a photochemical equilibrium. HLO uses the same formulation as Cariolle and D´equ´e (1986), but with the coefficients derived from both an ozone reanalysis and from a full chemistry model. Chapter 11 ‘Climatological data’ describes the distributions of climatological fields, including land sea mask, orography (mean, standard deviation, anisotropy, orientation and slope), land use, soil type, vegetation type, leaf area index, albedo, aerosols, ozone and trace gases. 1.2 OVERVIEW OF THE CODE The physics parametrization package is called by the IFS after the explicit dynamical computations. The physics computations are performed only in the vertical. The input information for the physics consists of the values of the gridbox mean prognostic variables (wind components u/v, temperature T, specific humidity q, cloud fraction a, water contents for cloud liquid ql , cloud ice qi, rain qr, and snow qs), the provisional dynamical tendencies for the same variables and various surface fields, both fixed and variable. The time integration of the physics is based on the following: (i) It has to be compatible with the adiabatic part of the IFS. (ii) The tendencies from the different physical processes are computed in separate routines. (iii) As a general approach, the value of a prognostic variable is updated with the tendency from one process and the next process starts from this updated value, in what is usually referred to as the ‘method of fractional steps’ (details are different for different processes). (iv) Implicit schemes are used when needed for stability. CALLPAR is the routine that controls the physical parametrization package with the exception of the main radiation routine RADINTG. RADINTG controls the computation of the short-wave transmissivities and the long-wave fluxes. RADINTG is called outside CALLPAR because of the need to make the radiation space interpolation compatible with the distributed memory version of the IFS. GP MODEL is a high level routine that controls all computations in grid-point space. It calls both CALLPAR via interface routines EC PHYS DRV, EC PHYSG, EC PHYS and RADINTG via driver RADDRV. In CALLPAR the physics routines are called in the following order: IFS Documentation – Cy48r1 7 Chapter 1: Overview SURFRAD LAYER Computes radiative properties of the surface. RADFLUX LAYER Calls RADHEATN to compute the temperature tendencies and the downward radiation fluxes at the surface with updated (every time-step) values for the zenith angle. GWDRAG LAYER Computes the tendencies for u, v and T due to the parametrization of subgrid- scale orographic drag. It also computes subgrid orographic coefficients for use in VDFMAIN. TURBULENCE LAYER Calls VDFOUTER/VDFMAIN in two sub time steps for numerical stability. VDFMAIN computes the vertical exchange of u, v, T, q, a, ql , qi by turbulence. CONVECTION LAYER Interface to call CUCALLN/CUMASTRN that controls the computation of the tendencies for u, v, T, q, chemical tracers and the cloud detrainment term due to the parametrization of moist convective processes. CLOUD SATADJ Performs subgrid cloud saturation adjustment due to tendencies of T, q from dynamics, radiation, turbulence and convection. Computes tendencies of T, q, a, ql and qi. CLOUD LAYER Calls CLOUDSC to compute tendencies for T, q, a, ql , qi, qr and qs due to the parametrization of cloud and precipitation processes. GWDRAGWMS LAYER Calls GWDRAG WMS to compute the tendencies for u, v and T due to the parametrization of non-orographic gravity waves. METHOX Computes tendencies for q due to methane oxidation and water vapour photolysis. SURFTSTP LAYER Calls SURFTSTP to control the soil/surface scheme. STOCHPERT LAYER Optionally add stochastic perturbations to physics tendencies. O3CHEM Computes tendencies for O3 due to ozone chemistry. SLTEND LAYER Optionally average tendencies from radiation, convection and cloud in time and space along the semi-Lagrangian trajectory. 8 IFS Documentation – Cy48r1 Part IV: Physical Processes Chapter 2 Radiation Table of contents 2.1 Introduction 2.2 The pre-Cy32r2 Short-Wave Radiation Scheme 2.2.1 Spectral integration 2.2.2 Vertical integration 2.2.3 Multiple reflections between layers 2.2.4 Pre-CY32R2 cloud short-wave optical properties 2.3 The pre-Cy22r3 Long-Wave Radiation Scheme 2.3.1 The pre-cycle Cy22r3 scheme 2.3.2 Vertical integration 2.3.3 Spectral integration 2.3.4 The incorporation of the effects of clouds 2.4 Inputs to the radiation scheme 2.4.1 Model variables 2.4.2 Clouds 2.4.3 Aerosols 2.4.4 Greenhouse and other gases 2.4.5 Surface albedo and emissivity 2.4.6 Solar irradiance and zenith angle 2.5 The ecRad radiation scheme 2.5.1 Overview 2.5.2 Gas optics 2.5.3 Aerosol optics 2.5.4 Cloud optics 2.5.5 Solver 2.6 Spatial and temporal interpolation 2.6.1 The radiation grid 2.6.2 Approximate radiation updates: short-wave 2.6.3 Approximate radiation updates: long-wave 2.7 Radiation diagnostics 2.7.1 Conventions 2.7.2 Surface fluxes 2.7.3 Top-of-atmosphere fluxes 2.7.4 Other variables 2.8 Code Appendix A. List of symbols 2.1 INTRODUCTION Table 2.1 gives the timeline of the major changes affecting the representation of the radiation transfer (RT) in the ECMWF model since the late 80s; see Table 2 of Hogan et al. (2017) for additional details. The original Morcrette (1991) radiation scheme, which was operational in the high-resolution 10-day IFS Documentation – Cy48r1 9 Chapter 2: Radiation Table 2.1 Major changes in the representation of radiation transfer in the ECMWF forecasting system. Cycle Implementation Description date Morcrette scheme SPM 32 02/05/1989 RT schemes from Univ. Lille SPM 46 01/02/1993 Optical properties for ice and mixed phase clouds IFS 14r3 13/02/1996 Revised LW and SW absorption coefficients from HITRAN’92 IFS 16r2 15/05/1997 Voigt profile in long-wave RT scheme IFS 16r4 27/08/1997 Revised ocean albedo from ERBE IFS 18r3 16/12/1997 Revised LW and SW absorption coefficients from HITRAN’96 IFS 18r5 01/04/1998 Seasonal land albedo from ERBE Morcrette scheme with RRTMGLW IFS 22r3 27/06/2000 RRTMGLW as long-wave RT scheme short-wave RT scheme with 4 spectral intervals IFS 23r4 12/06/2001 Hourly, instead of 3-hourly, calls to RT code during data assimilation cycle IFS 25r1 09/04/2002 Short-wave RT scheme with 6 spectral intervals IFS 26r3 07/10/2003 New aerosol climatology adapted from Tegen et al. (1997), new radiation grid IFS 28r3 28/09/2004 Radiation called hourly in high resolution forecasts McRad scheme IFS 32r2 05/06/2007 McICA approach to RT with RRTMGLW and RRTMGSW revised cloud optical properties, 60-km MODIS land albedo IFS 35r3 08/09/2009 GEMS climatology for ozone and other gases IFS 41r1 12/05/2015 MACC climatology for ozone and other gases, Updated RRTMG coefficients, 5-km MODIS land albedo IFS 41r2 08/03/2016 Approximate updates to mitigate coastal temperature errors IFS 43r1 22/11/2016 CAMS climatology for ozone, total solar irradiance (TSI) set to 1361 ± solar cycle ecRad scheme IFS 43r3 11/07/2017 Replace McRad with ecRad (gas optics unchanged) IFS 46r1 11/06/2020 Turned on long-wave scattering for clouds IFS 47r1 30/06/2020 CMIP6 for long-term trends in greenhouse gases and TSI, MODIS land albedo uses solar zenith angle dependence IFS 48r1 Radiatively interactive prognostic ozone from ‘HLO’ scheme forecasts from May 1989 to June 2000, is still described in this document as it forms the basis for the linearized schemes for short-wave (Section 2.2) and long-wave (Section 2.3) radiation transfer used in data assimilation. In 2007, with Cycle 32r2 of the ECMWF IFS operational libraries, a new approach to the inclusion of sub- grid cloud effects on radiation fields (the Monte-Carlo Independent Column Approximation, McICA) was introduced in the IFS (the McRad radiation scheme, Morcrette et al., 2008a). McRad also included both the long-wave and short-wave gas-optics treatment from the Rapid Radiation Transfer Model for GCMs (RRTMG), originally developed at AER, Inc. (Mlawer et al., 1997; Iacono et al., 2008). In 2017, the ecRad radiation scheme was introduced (Hogan and Bozzo, 2018), which improved computational efficiency, as well as providing much more flexibility with regard to the treatment of gas, aerosol and cloud optical properties, and the selection of different solvers for treating cloud structure. An offline version of ecRad is available from GitHub (https://github.com/ecmwf-ifs/ecrad) and has since been incorporated into other models such as ICON, AROME and MesoNH. The inputs to the radiation scheme are described in Section 2.4, followed by a description of ecRad in Section 2.5 10 IFS Documentation – Cy48r1 Part IV: Physical Processes The radiation scheme has traditionally been a relatively expensive part of the model, so the IFS reduces the overall cost by the use of a coarser-resolution radiation grid (Morcrette et al., 2008b) and by calling the full radiation scheme only every hour rather than every timestep. These reduced radiation configurations are described in Section 2.6 together with the approximate update scheme used for obtaining the radiative fluxes at every grid point and every time step for the relevant instantaneous temperature profile and solar zenith angle. Finally, Section 2.7 describes the radiation diagnostics that can be output from the IFS and Section 2.8 outlines the structure of the code. The radiation scheme is configured via the NAERAD namelist group. Throughout this section when namelist options are described they are parameters in this group. Data files used by the radiation scheme are read from the ifsdata subdirectory of the directory provided in the DATA environment variable in the IFS. At ECMWF the files for Cycle 48r1 can be found in /home/rdx/data/48/ifsdata. 2.2 THE PRE-CY32R2 SHORT-WAVE RADIATION SCHEME Please note that this section refers to a very old version of the radiation scheme that is only used in the linearized physics used in data assimilation. The symbols used here may differ from those in the other sections of this chapter. The rate of atmospheric heating by absorption and scattering of short-wave radiation is ∂T ∂t = − g cp ∂FSW ∂p (2.1) where FSW is the net total short-wave flux (the subscript SW will be omitted in the remainder of this section). F (δ) = Z ∞ 0 dν Z 2π 0 dϕ Z +1 −1 μLν(δ, μ, ϕ) dμ  (2.2) is the diffuse radiance at wavenumber ν, in a direction given by the azimuth angle, ϕ, and the zenith angle, θ, with μ = cos θ. In (2.2), we assume a plane parallel atmosphere, and the vertical coordinate is the optical depth δ, a convenient variable when the energy source is outside the medium δ(p) = Z 0 p βext v (p′) dp′ (2.3) βext ν (p) is the extinction coefficient, equal to the sum of the scattering coefficient βsca ν of the aerosol (or cloud particle absorption coefficient βabs ν ) and the purely molecular absorption coefficient kν. The diffuse radiance Lν is governed by the radiation transfer equation μ dLν(δ, μ, ϕ) dδ = Lν(δ, μ, ϕ) − ϖν(δ) 4 Pν(δ, μ, ϕ, μ0, ϕ0)E0 ν exp(−δ/μr ) − ϖν(δ) 4 Z 2π 0 dϕ′ Z +1 −1 Φν(δ, μ, ϕ, μ′, ϕ′)Lν(δ, μ′, ϕ′) dμ′  (2.4) E0 ν is the incident solar irradiance in the direction μ0 = cos θ0, ϖν, is the single scattering albedo (= βsca ν /kv) and Φ(δ, μ, ϕ, μ′, ϕ′) is the scattering phase function which defines the probability that radiation coming from direction (μ′, ϕ′) is scattered in direction (μ, ϕ). The short-wave part of the scheme, originally developed by Fouquart and Bonnel (1980) solves the radiation transfer equation and integrates the fluxes over the whole short-wave spectrum between 0.2 and 4 μm. Upward and downward fluxes are obtained from the reflectances and transmittances of the layers, and the photon- path-distribution method allows to separate the parametrization of the scattering processes from that of the molecular absorption. 2.2.1 Spectral integration Solar radiation is attenuated by absorbing gases, mainly water vapour, uniformly mixed gases (oxygen, carbon dioxide, methane, nitrous oxide) and ozone, and scattered by molecules (Rayleigh scattering), aerosols and cloud particles. Since scattering and molecular absorption occur simultaneously, the exact IFS Documentation – Cy48r1 11 Chapter 2: Radiation amount of absorber along the photon path length is unknown, and band models of the transmission function cannot be used directly as in long-wave radiation transfer (see Section 2.2). The approach of the photon path distribution method is to calculate the probability Π(U ) dU that a photon contributing to the flux Fcons in the conservative case (i.e., no absorption, ων = 1, kν = 0) has encountered an absorber amount between U and U + dU . With this distribution, the radiative flux at wavenumber v is related to Fcons by Fv = Fcons Z ∞ 0 Π(U ) exp(−kνU ) dU (2.5) and the flux averaged over the spectral interval ∆ν can then be calculated with the help of any band model of the transmission function t∆ν F = 1 ∆ν Z ∆ν Fν dν = Fcons Z ∞ 0 Π(U )t∆ν(U ) dν (2.6) To find the distribution function Π(U ), the scattering problem is solved first, by any method, for a set of arbitrarily fixed absorption coefficients k1, thus giving a set of simulated fluxes Fk1 . An inverse Laplace transform is then performed on (2.5) (Fouquart, 1974). The main advantage of the method is that the actual distribution Π(U ) is smooth enough that (2.5) gives accurate results even if Π(U ) itself is not known accurately. In fact, Π(U ) needs not be calculated explicitly as the spectrally integrated fluxes are F = Fconst∆v(⟨U ⟩) in the limiting case of weak absorption F = Fconst∆v(⟨U 1/2⟩) in the limiting case of strong absorption where ⟨U ⟩ = R ∞ 0 Π(U )U dU and ⟨U 1/2⟩ = R ∞ 0 Π(U )U 1/2 dU . The atmospheric absorption in the water vapour bands is generally strong, and the scheme determines an effective absorber amount Ue between ⟨U ⟩ and ⟨U 1/2⟩ derived from Ue = ln(Fke /Fcons)/ke (2.7) where ke is an absorption coefficient chosen to approximate the spectrally averaged transmission of the clear sky atmosphere ke = 1 Utot/μ0 ln(t∆v(Utot/μ0)) (2.8) where Utot is the total amount of absorber in a vertical column and μ0 = cos θ0. Once the effective absorber amounts of H2O and uniformly mixed gases are found, the transmission functions are computed using Pad´e approximants t∆ν(U ) = ∑N i=0 aiU i−1 ∑N j=0 bjU j−1 (2.9) Absorption by ozone is also taken into account, but since ozone is located at low pressure levels for which molecular scattering is small and Mie scattering is negligible, interactions between scattering processes and ozone absorption are neglected. Transmission through ozone is computed using (2.6) where UO3 the amount of ozone is U d O3 = M Z 0 p dUO3 for the downward transmission of the direct solar beam U u O3 = r Z 0 ps dUO3 + U d O3 (psurf) for the upward transmission of the diffuse radiation r = 1.66 is the diffusivity factor (see Section 2.2), and M is the magnification factor (Rodgers, 1967) used instead of r to account for the sphericity of the atmosphere at very small solar elevations M = 35/ q μ2 0 + 1 (2.10) To perform the spectral integration, it is convenient to discretize the solar spectrum into subintervals in which the surface reflectance, molecular absorption characteristics, and cloud optical properties can be 12 IFS Documentation – Cy48r1 Part IV: Physical Processes considered as constants. One of the main causes for such a spectral variation is the sharp increase in the reflectivity of the vegetation in the near-infrared. Also, water vapour does not absorb below 0.69 μm nor do liquid water clouds. Till June 2000, the ECMWF short-wave scheme considered only two spectral intervals, one for the visible (0.2–0.69 μm), one for the near-infrared (0.69–4.00 μm) parts of the solar spectrum. From June 2000 to April 2002, the near-infrared interval was sub-divided into three intervals (0.69–1.19–2.38–4.00 μm) to account better for the spectral variations of the cloud optical properties. Till April 2002, all the molecular absorption coefficients (for O3, H2O, uniformly mixed gases) were derived from statistical models of the transmission function using spectroscopic parameters derived from various versions of the HITRAN database (Rothman et al., 1986, 1992). In April 2002, following the recomputation of all the molecular absorption coefficients from an updated version of the short- wave line-by-line model of Dubuisson et al. (1996) using spectroscopic data from HAWKS (2000), the ultraviolet and visible part of the spectrum are now considered in three spectral intervals (0.20–0.25– 0.69 μm) making the scheme having a total of six spectral intervals over which the aerosol and cloud optical properties are also defined. The cut-off at 0.69 μm allows the scheme to be more computational efficient, in as much as the interactions between gaseous absorption (by water vapour and uniformly mixed gases) and scattering processes are accounted for only in the near-infrared interval(s). 2.2.2 Vertical integration Considering an atmosphere where a fraction Ctot cld (as seen from the surface or the top of the atmosphere) is covered by clouds (the fraction Ctot cld depends on which cloud-overlap assumption is assumed for the calculations), the final fluxes are given as a weighted average of the fluxes in the clear sky and in the cloudy fractions of the column F −(j) = Ctot cldF − cld(j) + (1 − Ctot cld)F − clr where the subscripts ‘clr’ and ‘cld’ refer to the clear-sky and cloudy fractions of the layer, respectively. In contrast to the scheme of Geleyn and Hollingsworth (1979), the fluxes are not obtained through the solution of a system of linear equations in a matrix form. Rather, assuming an atmosphere divided into homogeneous layers, the upward and downward fluxes at a given layer interface j are given by F −(j) = F0 N ∏ k=j Tbot(k) F +(j) = F −(j)Rtop(j − 1) (2.11) where Rtop(j) and Tbot(j) are the reflectance at the top and the transmittance at the bottom of the jth layer. Computation of the values of Rtop starts at the surface and works upwards, whereas determining values of Tbot starts at the top of the atmosphere and works downward. Rtop and Tbot account for the presence of cloud in the layer by using Rtop = CcldRcld + (1 − Ccld)Rclr Tbot = CcldTcld + (1 − Ccld)Tclr (2.12) where Ccld is the cloud fractional coverage of the layer within the cloudy fraction Ctot cld of the column. (a) Cloudy fraction layer Rtcdy and Rbcdy are the reflectance at the top and transmittance at the bottom of the cloudy fraction of the layer calculated with the Delta-Eddington approximation. Given δc, δa, and δg, the optical thicknesses for the cloud, the aerosol and the molecular absorption of the gases (= keU), respectively, and gc and ga the cloud and aerosol asymmetry factors, Rtcdy and Rbcdy are calculated as functions of the total optical thickness of the layer δ = δc + δa + δg (2.13) of the total single scattering albedo ϖ∗ = δc + δa δc + δa + δg (2.14) IFS Documentation – Cy48r1 13 Chapter 2: Radiation of the total asymmetry factor g∗ = δc δc + δa gc + δa δc + δa ga (2.15) of the reflectance R− of the underlying medium (surface or layers below the jth interface), and of the cosine of an effective solar zenith angle μeff(j) which accounts for the decrease of the direct solar beam and the corresponding increase of the diffuse part of the downward radiation by the upper scattering layers μeff(j) = [(1 − Ceff cld(j))/μ + rCeff cld(j)]−1 (2.16) with Ceff cld(j) the effective total cloudiness over level j Ceff cld(j) = 1 − N ∏ i=j+1 (1 − Ccld(i)E(i)) (2.17) and E(i) = 1 − exp  − (1 − ϖc(i)gc(i)2)δc(i) μ  (2.18) δc(i), ϖc(i) and gc(i) are the optical thickness, single scattering albedo and asymmetry factor of the cloud in the ith layer, and r is the diffusivity factor. The scheme follows the Eddington approximation first proposed by Shettle and Weinman (1970), then modified by Joseph et al. (1976) to account more accurately for the large fraction of radiation directly transmitted in the forward scattering peak in case of highly asymmetric phase functions. Eddington’s approximation assumes that, in a scattering medium of optical thickness δ∗, of single scattering albedo ω, and of asymmetry factor g, the radiance L can be written as L(δ, μ) = L0(δ) + μL1(δ) (2.19) In that case, when the phase function is expanded as a series of associated Legendre functions, all terms of order greater than one vanish when (2.2) is integrated over μ and ϕ. The phase function is therefore given by P(Θ) = 1 + β1(Θ)μ where Θ is the angle between incident and scattered radiances. The integral in (2.2) thus becomes Z 2π 0 dϕ′ Z +1 −1 p(μ, ϕ, μ′, ϕ′)L(μ′, ϕ′) dμ′  = 4π(L0 + πL1) (2.20) where g = β1 3 = 1 2 Z +1 −1 P(Θ)μ dμ is the asymmetry factor. Using (2.20) in (2.2) after integrating over μ and dividing by 2π, we get μ d dδ (L0 + μL1) = −(L0 + μL1) + ϖ(L0 + gμL1) + 1/4ϖF0 exp(−δ/μ0)(1 + 3gμ0μ) (2.21) We obtain a pair of equations for L0 and L1 by integrating (2.21) over μ dL0 dδ = −3(1 − ϖ)L0 + 3 4 ϖF0 exp(−δ/μ0) dL1 dδ = −(1 − ϖg)L1 + 3 4 ϖgμ0F0 exp(−δ/μ0) (2.22) For the cloudy layer assumed non-conservative (ϖ < 1), the solutions to (2.21) and (2.22), for 0 ≤ δ ≤ δ∗, are L0(δ) = C1 exp(−Kδ) + C2 exp(+Kδ) − α exp(−δ/μ0) L1(δ) = P{C1 exp(−Kδ) − C2 exp(+Kδ) − β exp(−δ/μ0)} (2.23) 14 IFS Documentation – Cy48r1 Part IV: Physical Processes where K = {3(1 − ϖ)(1 − ϖg)}1/2 P = {3(1 − ϖ)/(1 − ϖg)}1/2 α = 3ϖF0μ0{1 + 3g(1 − ϖ)}/{4(1 − K2μ2 0)} β = 3ϖF0μ0{1 + 3g(1 − ϖ)μ2 0}/(4(1 − K2μ2 0)) The two boundary conditions allow to solve the system for C1 and C2; the downward directed diffuse flux at the top of the atmosphere is zero, that is F −(0) =  L0(0) + 2 3 L1(0)  = 0 which translates into (1 + 2P/3)C1 + (1 − 2P/3)C2 = α + 2β/3 (2.24) The upward directed flux at the bottom of the layer is equal to the product of the downward directed diffuse and direct fluxes and the corresponding diffuse and direct reflectance (Rd and R−, respectively) of the underlying medium F +(δ∗) =  L0(δ∗) − 2 3 L1(δ∗)  = R−  L0(δ∗) + 2 3 L1(δ∗)  + Rdμ0F0 exp(−δ∗/μ0) which translates into {1 − R− − 2(1 + R−)P/3}C1 exp(−Kδ∗) + {1 − R− + 2(1 + R−)P/3}C2(+Kδ∗) = {(1 − R−)α − 2(1 + R−)β/3 + Rdμ0F0} exp(−δ∗/μ0) (2.25) In the Delta-Eddington approximation, the phase function is approximated by a Dirac delta function forward-scatter peak and a two-term expansion of the phase function P(θ) = 2 f (1 − μ) + (1 − f )(1 + 3g′μ) where f is the fractional scattering into the forward peak and g′ the asymmetry factor of the truncated phase function. As shown by Joseph et al. (1976), these parameters are f = g2 g′ = g/(g + 1) (2.26) The solution of the Eddington’s equations remains the same provided that the total optical thickness, single scattering albedo and asymmetry factor entering (2.21) and (2.25) take their transformed values δ′∗ = (1 + ϖ f )δ∗ ω′ = (1 − f )ϖ 1 − ϖ f (2.27) Practically, the optical thickness, single scattering albedo, asymmetry factor and solar zenith angle entering (2.21)–(2.25) are δ∗, ϖ∗, g∗ and μeff defined in (2.15) and (2.16). (b) Clear-sky fraction of the layers In the clear-sky part of the atmosphere, the short-wave scheme accounts for scattering and absorption by molecules and aerosols. The following calculations are practically done twice, once for the clear- sky fraction (1 − Ctot cld) of the atmospheric column μ with equal to μ0, simply modified for the effect of IFS Documentation – Cy48r1 15 Chapter 2: Radiation Rayleigh and aerosol scattering, the second time for the clear-sky fraction of each individual layer within the fraction Ctot cld of the atmospheric column containing clouds, with μ equal to μe. As the optical thickness for both Rayleigh and aerosol scattering is small, Rclr(j − 1) and Tclr(j), the reflectance at the top and transmittance at the bottom of the jth layer can be calculated using respectively a first- and a second-order expansion of the analytical solutions of the two-stream equations similar to that of Coakley Jr. and Chylek (1975). For Rayleigh scattering, the optical thickness, single scattering albedo and asymmetry factor are respectively δR, ϖR = 1 and gR = 0, so that RR = δR 2μ + δR TR = 2μ (2μ + δR) (2.28) The optical thickness δR of an atmospheric layer is simply δR = δ∗{p(j) − p(j − 1)}/psurf (2.29) where δ∗ R is the Rayleigh optical thickness of the whole atmosphere parametrized as a function of the solar zenith angle (Deschamps et al., 1983) For aerosol scattering and absorption, the optical thickness, single scattering albedo and asymmetry factor are respectively δa, ϖa, with 1 − ϖa ≪ 1 and ga, so that den = 1 + {1 − ϖa + back(μe)ϖa}(δa/μe) + (1 − ϖa){1 − ϖa + 2 back(μe)ϖa}(δ2 a /μ2 e ) (2.30) R(μe) = (back(μe)ϖaδa)/μa den T (μe) = 1/den (2.31) where back(μe) = (2 − 3μe ga)/4 is the backscattering factor. Practically, Rclr and Tclr are computed using (2.31) and the combined effect of aerosol and Rayleigh scattering comes from using modified parameters corresponding to the addition of the two scatterers with provision for the highly asymmetric aerosol phase function through Delta-approximation of the forward scattering peak (as in (2.22) and (2.23)). δ+ = δR + δa(1 − ϖa g2 a ) g+ = ga 1 + ga δa (δR + δa) ϖ+ = δR δa + δa ϖR + δa δR + δa ϖa(1 − g2 a) 1 − ϖa g2 a (2.32) As for their cloudy counterparts, Rclr and T˙clr must account for the multiple reflections due to the layers underneath Rclr = R(μe) + R−T (μe)(1 − R∗R−) (2.33) and R− is the reflectance of the underlying medium R− = Rt(j − 1) and r is the diffusivity factor. Since interactions between molecular absorption and Rayleigh and aerosol scattering are negligible, the radiative fluxes in a clear-sky atmosphere are simply those calculated from (2.9) and (2.27) attenuated by the gaseous transmissions (2.7). 2.2.3 Multiple reflections between layers To deal properly with the multiple reflections between the surface and the cloud layers, it should be necessary to separate the contribution of each individual reflecting surface to the layer reflectance and 16 IFS Documentation – Cy48r1 Part IV: Physical Processes transmittances in as much as each such surface gives rise to a particular distribution of absorber amount. In the case of an atmosphere including N cloud layers, the reflected light above the highest cloud consists of photons directly reflected by the highest cloud without interaction with the underlying atmosphere, and of photons that have passed through this cloud layer and undergone at least one reflection on the underlying atmosphere. In fact, (2.4) should be written F = N ∑ i=0 Fcl Z ∞ 0 P1(U )t∆v(U )dv (2.34) where Fcl and P1(U ) are the conservative fluxes and the distributions of absorber amount corresponding to the different reflecting surfaces. Fouquart and Bonnel (1980) have shown that a very good approximation to this problem is obtained by evaluating the reflectance and transmittance of each layer (using (2.21) and (2.27)) assuming successively a non-reflecting underlying medium (R− = 0), then a reflecting underlying medium (R−̸ = 0). First calculations provide the contribution to reflectance and transmittance of those photons interacting only with the layer into consideration, whereas the second ones give the contribution of the photons with interactions also outside the layer itself. From those two sets of layer reflectance and transmittances (Tt0, Tb0) and (Rt̸ =, Tb̸ =) respectively, effective absorber amounts to be applied to computing the transmission functions for upward and downward fluxes are then derived using (2.5) and starting from the surface and working the formulas upward U − e0 = ln (Tb0/Tbc)/ke U − e̸ = = ln (Tb̸ =/Tbc)/ke U + e0 = ln (Rt0/Rtc)/ke U + e̸ = = ln (Rt̸ =/Rtc)/ke (2.35) where Rtc and Tbc are the layer reflectance and transmittance corresponding to a conservative scattering medium. Finally the upward and downward fluxes are obtained as F +(j) = F0{Rt0t∆ν(U + e0) + (Rt̸ = − Rt0)t∆ν(U + e̸ =)} (2.36) F −(j) = F0{Tb0t∆ν(U + e0) + (Tb̸ = − Tb0)t∆ν(U − e̸ =)} (2.37) 2.2.4 Pre-CY32R2 cloud short-wave optical properties As seen in Subsection 2.2.2(a), the cloud radiative properties depend on three different parameters: the optical thickness δc, the asymmetry factor gc, and the single scattering albedo ϖc. In this older scheme the cloud optical properties were derived from Fouquart (1987) for the water clouds, and Ebert and Curry (1992) for the ice clouds. The optical thickness δc is related to the cloud liquid water amount ULWP by δc = 3ULWP 2re where re is the mean effective radius of the size distribution of the cloud water droplets. In an older formulation of this scheme, re was parametrized as a linear function of height from 10 μm at the surface to 45 μm at the top of the atmosphere, in an empirical attempt at dealing with the variation of water cloud type with height. Smaller water droplets are observed in low-level stratiform clouds whereas larger droplets are found in mid-level cumuliform water clouds. This parametrisation was then revised, using the formulation of Martin et al. (1994) to compute the effective radius of the liquid water cloud particles from the cloud liquid water content and specified concentrations of cloud condensation IFS Documentation – Cy48r1 17 Chapter 2: Radiation nuclei over land and ocean. For ice clouds, the effective dimension of the cloud particles was diagnosed from temperature using a revision of the formulation by Ou and Liou (1995). In the two-, four-, and six-spectral interval versions of the short-wave radiation scheme, the optical properties of liquid water clouds were defined from Fouquart (1987) and those for ice clouds from Ebert and Curry (1992). Alternative optical properties were also available for liquid water clouds (Slingo, 1989) and ice clouds (Fu, 1996). 2.3 THE PRE-CY22R3 LONG-WAVE RADIATION SCHEME Please note that this section refers to a very old version of the radiation scheme that is only used in the linearized physics used in data assimilation. The symbols used here may differ from those in the other sections of this chapter. As already noted, since cycle Cy22r3, two long-wave radiation schemes have been available in the ECMWF model, the pre-cycle Cy22r3 by Morcrette (1991), and the current long-wave radiation transfer scheme, the Rapid Radiation Transfer Model (RRTMG). The rate of atmospheric cooling by emission-absorption of long-wave radiation is ∂T ∂t = g cp ∂FLW ∂p (2.38) where FLW is the net long-wave radiation flux (the subscript ‘LW’ is omitted in the remainder of this section). Assuming a non-scattering atmosphere in local thermodynamic equilibrium, F is given by F = Z 1 −1 μ dμ Z ∞ 0 dv  Lv(psurf, μ)tv(psurf, p, μ) + Z 0 p′ =psurf Lv(p′, μ) dtv  (2.39) where Lv(p, μ) is the monochromatic radiance at wavenumber v at level p, propagating in a direction θ (the angle that this direction makes with the vertical), where μ = cos θ and tv(p, p′; r) is the monochromatic transmission through a layer whose limits are at p and p′ seen under the same angle θ, with r = sec θ. The subscript ‘surf’ refers to the earth’s surface. 2.3.1 The pre-cycle Cy22r3 scheme After separating the upward and downward components (indicated by superscripts + and −, respectively), and integrating by parts, we obtain the radiation transfer equation as it is actually estimated in the long-wave part of the radiation code F + v (p) = [Bv(Tsurf) − Bv(T0+)]tv(psurf, p; r) + Bv(T(p)) + Z p p′ =psurf tv(p, p′; r) dBv F − v (p) = [Bv(T∞) − Bv(Ttop)]tv(p, 0; r) + Bv(T(p)) + Z 0 p′ =p tv(p′, p; r) dBv (2.40) where, taking benefit of the isotropic nature of the long-wave radiation, the radiance Lv of (2.39) has been replaced by the Planck function Bv(T) in units of flux, Wm−2 (here, and elsewhere, Bv is assumed to always includes the π factor). Tsurf is the surface temperature, T0+ that of the air just above the surface, T(p) is the temperature at pressure-level p, Ttop that at the top of the atmospheric model. The transmission tv is evaluated as the radiance transmission in a direction θ to the vertical such that r = sec θ is the diffusivity factor (Elsasser, 1942). Such an approximation for the integration over the angle is usual in radiative transfer calculations, and tests on the validity of this approximation have been presented by Rodgers and Walshaw (1966) and Liu and Schmetz (1988) among others. The use of the diffusivity factor gives cooling rates within 2% of those obtained with a 4-point Gaussian quadrature. 2.3.2 Vertical integration The integrals in (2.40) are evaluated numerically, after discretization over the vertical grid, considering the atmosphere as a pile of homogeneous layers. As the cooling rate is strongly dependent on local 18 IFS Documentation – Cy48r1 Part IV: Physical Processes conditions of temperature and pressure, and energy is mainly exchanged with the layers adjacent to the level where fluxes are calculated, the contribution of the distant layers is simply computed using a trapezoidal rule integration, but the contribution of the adjacent layers is evaluated with a 2-point Gaussian quadrature, thus at the ith level Z pi p′ =psurf tv(p, p′; r) dBv = 2 ∑ l=1 dBv(l)wl tv(pi, pl ; r) + 1 2 i−2 ∑ j=1 dBv(j)[tv(pi, pj; r) + tv(pi, pj−1; r)] (2.41) where pl is the pressure corresponding to the Gaussian root and wl is the Gaussian weight. dBv(j) and dBv(l) are the Planck function gradients calculated between two interfaces, and between mid-layer and interface, respectively. 2.3.3 Spectral integration The integration over wavenumber v is performed using a band emissivity method, as first discussed by Rodgers (1967). The long-wave spectrum is divided into six spectral regions. (i) 0–350 cm−1 and 1450–1880 cm−1 (ii) 500–800 cm−1 (iii) 800–970 cm−1 and 1110–1250 cm−1 (iv) 970–1110 cm−1 (v) 350–500 cm−1 (vi) 1250–1450 cm−1 and 1880–2820 cm−1 corresponding to the centres of the rotation and vibration-rotation bands of H2O, the 15 μm band of CO2, the atmospheric window, the 9.6 μm band of O3, the 25 μm ‘window’ region, and the wings of the vibration-rotation band of H2O, respectively. Over these spectral regions, band fluxes are evaluated with the help of band transmissivities pre-calculated from the narrow-band model of Morcrette and Fouquart (1985) – See Appendix of Morcrette et al. (1986) for details. Integration of (2.40) over wavenumber ν within the kth spectral region gives the upward and downward fluxes as F + k (p) = {Bk(Tsurf) − Bk(T0+)}tBk {rU (psurf, p), TU (psurf, p)} + Bk(Tp) + Z p p′ =psurf tdBk {rU (p, p′), TU (p, p′)} dBk (2.42) F − k (p) = {Bk(T0) − Bk(T∞)}tBk {rU (p, 0), TU (p, 0)} − Bk(Tp) − Z 0 p′ =p tdBk {rU (p′, p), TU (p′, p)} dBk (2.43) The formulation accounts for the different temperature dependencies involved in atmospheric flux calculations, namely that on Tp, the temperature at the level where fluxes are calculated, and that on TU , the temperature that governs the transmission through the temperature dependence of the intensity and half-widths of the lines absorbing in the concerned spectral region. The band transmissivities are non- isothermal accounting for the temperature dependence that arises from the wavenumber integration of the product of the monochromatic absorption and the Planck function. Two normalized band transmissivities are used for each absorber in a given spectral region: the first one for calculating the first right-hand-side term in (2.40), involving the boundaries; it corresponds to the weighted average of the transmission function by the Planck function tB(U p, Tp, TU ) = R v2 v1 Bv(Tp)tv(U p, TU ) dv R v2 v1 Bv(Tp) dv (2.44) the second one for calculating the integral term in (2.40) is the weighted average of the transmission function by the derivative of the Planck function tdB(U p, Tp, TU ) = R v2 v1 {dB(Tp)/dT}tv(U p, TU ) dv R v2 v1 {dB(Tp)/dT} dv (2.45) IFS Documentation – Cy48r1 19 Chapter 2: Radiation where U p is the pressure weighted amount of absorber. The effect on absorption of the Doppler broadening of the lines (important only for pressure lower than 10 hPa) is included simply using the pressure correction method of Fels (1979). A finite line width (assumed to represent the Doppler half-width of the line) is retained under low pressure conditions where the pure Lorentz line width (proportional to pressure) would normally become negligible (Giorgetta and Morcrette, 1995). In the scheme, the actual dependence on Tp is carried out explicitly in the Planck functions integrated over the spectral regions. Although normalized relative to B(Tp) or dB(Tp)/dT, the transmissivities still depend on TU , both through Wien’s displacement of the maximum of the Planck function with temperature and through the temperature dependence of the absorption coefficients. For computational efficiency, the transmissivities have been developed into Pad´e approximants t(U p, Tu) = ∑2 i=0 ciU i/2 eff ∑2 j=0 djU i/2 eff (2.46) where Ueff = r(U p)Ψ(TU , U p) is an effective amount of absorber which incorporates the diffusivity factor r, the weighting of the absorber amount by pressure U p, and the temperature dependence of the absorption coefficients. The function Ψ(TU , U p) takes the form Ψ(TU , U p) = exp[a(U p)(TU − 250) + b(U p)(TU − 250)2] (2.47) The temperature dependence due to Wien’s law is incorporated although there is no explicit variation of the coefficients ci and dj with temperature. These coefficients have been computed for temperatures between 187.5 and 312.5 K with a 12.5 K step, and transmissivities corresponding to the reference temperature the closest to the pressure weighted temperature TU are actually used in the scheme. 2.3.4 The incorporation of the effects of clouds The incorporation of the effects of clouds on the long-wave fluxes follows the treatment discussed by Washington and Williamson (1977). Whatever the state of the cloudiness of the atmosphere, the scheme starts by calculating the fluxes corresponding to a clear-sky atmosphere and stores the terms of the energy exchange between the different levels (the integrals in (2.40)) Let F + 0 (i) and F − 0 (i) be the upward and downward clear-sky fluxes. For any cloud layer actually present in the atmosphere, the scheme then evaluates the fluxes assuming a unique overcast cloud of emissivity unity. Let F + n (i) and F + n (i) the upward and downward fluxes when such a cloud is present in the nth layer of the atmosphere. Downward fluxes above the cloud, and upward fluxes below the cloud, are assumed to be given by the clear-sky values F + n (i) = F + 0 (i) for i ≤ n F − n (i) = F − 0 (i) for i > n (2.48) Upward fluxes above the cloud (F + n (k) for k ≤ n + 1) and downward fluxes below it (F − n (k) for k > n) can be expressed with expressions similar to (2.41) provided the boundary terms are now replaced by terms corresponding to possible temperature discontinuities between the cloud and the surrounding air F + n (k) = {F + cld − B(n + 1)}t(pk, pn+1; r) + B(k) + Z pk p′ =pn−1 t(pk, p′; r) dB F − n (k) = {F − cld − B(n)}t(pk, pn; r) + B(k) + Z pn p′ =pk t(pk, p′; r) dB (2.49) where B(i) is now the total Planck function (integrated over the whole long-wave spectrum) at level i, and F + cld and F − cld are the long-wave fluxes at the upper and lower boundaries of the cloud. Terms under the integrals correspond to exchange of energy between layers in clear-sky atmosphere and have already been computed in the first step of the calculations. This step is repeated for all cloudy layers. The fluxes for the actual atmosphere (with semi-transparent, fractional and/or multi-layered clouds) are 20 IFS Documentation – Cy48r1 Part IV: Physical Processes derived from a linear combination of the fluxes calculated in previous steps with some cloud overlap assumption in the case of clouds present in several layers. Let N be the index of the layer containing the highest cloud, Ccld(i)) the fractional cloud cover in layer i, with Ccld(0) = 1 for the upward flux at the surface, and with Ccld(N + 1) = 1 and F − N+1 = F − 0 to have the right boundary condition for downward fluxes above the highest cloud. Whereas the maximum and random overlap assumptions are also available in the code (Morcrette and Fouquart, 1986), the maximum-random overlap assumption is used in this version of the radiation scheme, and the cloudy upward F + and downward F − fluxes are obtained as F +(i) = F + 0 (i) for i = 1 F −(i) = Ccld(i − 1)F + i−1(i−) + i−2 ∑ n=0 Ccld(n)F + n (i) i−1 ∏ l=n+1 {1 − Ccld(l)} for 2 ≤ i ≤ N + 1 F +(i) = Ccld(N)F + N (i) + N−1 ∑ n=0 Ccld(n)F + n (i) N ∏ l=n+1 {1 − Ccld(l)} for i ≥ N + 2 (2.50) In the case of semi-transparent clouds, the fractional cloudiness entering the calculations is an effective cloud cover equal to the product of the emissivity due to the condensed water and the gases in the layer by the horizontal coverage of the cloud layer, with the emissivity, εcld, related to the condensed water amount by εcld = 1 − exp(−kabsULWP) (2.51) where kabs is the condensed water mass absorption coefficient (in m2kg−1) following Smith and Shi (1992). 2.4 INPUTS TO THE RADIATION SCHEME 2.4.1 Model variables Temperature values are needed at the boundaries of the layers (half-levels), where the fluxes are computed. This is because the long-wave equations are solved in each model layer assuming the Planck function varies linearly with optical depth through the layer. They are derived from the model’s full- level temperatures with a pressure-weighted interpolation: Tk+1/2 = Tk pk(pk+1 − pk+1/2) pk+1/2(pk+1 − pk) + Tk+1 pk+1/2(pk+1/2 − pk) pk+1/2(pk+1 − pk) . (2.52) Hogan and Bozzo (2016) investigated several options for how to treat the temperature of the base of lowest atmospheric layer, Tn+1/2 (where n is the number of model layers): (i) Extrapolate linearly from the temperature at the half-level at the top of this layer, Tn−1/2, and the full-level temperature, Tn, which has the advantage that it respects the mean temperature of the layer. This can lead to a large discontinuity with the skin temperature, Ts, but it should be noted that in reality the temperature is far from linear in the lowest 20 m of the atmosphere represented by the lowest model layer. (ii) Set it equal to the skin temperature Ts, thereby avoiding the discontinuity but not conserving the layer-mean temperature. (iii) Set it equal to half-way between options (i) and (ii). The impact on 2-m temperature forecasts between options (i) and (iii) can be as much as 4 K (Hogan and Bozzo, 2016), primarily due to the very strong 15-μm CO2 band which makes the downwelling long-wave radiation sensitive to atmospheric temperature very close to the surface. The IFS currently uses option (iii), which lies between options (i) and (ii) in terms of impact. A potential avenue for future model improvement would be to investigate how to represent a more realistic vertical profile of temperature than linear in the lowest model layer. IFS Documentation – Cy48r1 21 Chapter 2: Radiation 2.4.2 Clouds Cloud fraction and the mass mixing ratios of liquid, ice, rain and snow (ql , qi, qr and qs), are provided in all layers by the cloud scheme. These are used to calculate the effective radius for cloud water drops and ice particles required for the calculation of cloud-affected radiative fluxes. Note that ‘in-cloud’ water contents are used in all the effective radius equations (i.e. gridbox mean divided by the cloud fraction). Internally, the radiation scheme only knows about two cloud types, ‘liquid’ and ‘ice’. The liquid cloud mixing ratio fed to radiation comes only from the model’s ql variable, but rain mixing is used in the computation of liquid effective radius. The ice cloud mixing ratio fed to radiation is the sum of the model’s ice and snow, i.e. ql + qs. For both phases, effective radius is defined as re = 3V 4A , (2.53) where V is the total volume of liquid or ice in a given volume of air, and A is the total projected area of the particles in the same volume of air. Thus, re relates what the model provides (the mass of liquid or ice) to what radiation is primarily sensitive to when particles are much larger than the wavelength (the projected area). For liquid clouds this may be written as re = 3 4 LWC ρl A , (2.54) where LWC is the liquid water content and ρl is the density of liquid water. A similar expression may be used for ice clouds using IWC and the density of solid ice. (a) Cloud droplet effective radius The formulation of the liquid cloud drop effective radius follows Martin et al. (1994) and Wood (2000). The liquid water content, LWC (equal to ql ρair), of a population of spherical cloud droplets can be written as an integral of the volume of the water drops over the droplet size distribution: LWC = 4 3 πρw Z ∞ 0 nd(r)r3dr, where ρw is the density of water and nd(r)dr is the number density of droplets of radius between r and r + dr. The average droplet radius (or volume radius) rv in a unit volume of air with a droplet number density Nd can then be written as rv =  3LWC 4πρw Nd 1/3 . The effective radius, re, of a cloud droplet population can be expressed as a function of its volume radius r by a factor k dependent on the shape of the cloud droplet size spectrum: r3 e = r3 v k , so that re =  3LWC 4πρwkNd 1/3 . (2.55) The constant k can be formulated in terms of the relative dispersion D of the droplet size distribution. This is a normalised measure of the width of the distribution, equal to the standard deviation divided by the mean. The relationship is defined as k = (1 + D2)3 (1 + 3D2)2 . In the IFS, the spectral dispersion is given a value of D = 0.33 over ocean corresponding to k = 0.77. Over land, D = 0.43 and therefore k = 0.69. The values of the constant k are both close to the values given in Martin et al. (1994) from maritime and continental observations of stratocumulus. 22 IFS Documentation – Cy48r1 Part IV: Physical Processes Wood (2000) proposed an adjustment to the parametrization of effective radius in Martin et al. (1994) due to observations of an increased relative dispersion in the drop size spectra (a broader spectrum) when drizzle drops are present in the cloud. This formulation was introduced to the IFS in Cycle 36r4. The effective radius enhancement factor accounting for this increased dispersion follows Eq. 19 of Wood (2000): Ed = 1 + ϕ2/3 1 + 0.2(kr /k)1/3ϕ , where ϕ = RWC/LWC is the ratio of rain water content to liquid water content, kr represents the large drops and is set to a constant 0.222, and k for the small droplets is defined from the dispersion D above. The enhancement factor, Ed, multiplies the effective radius in Eq. 2.55 and the in-cloud liquid water content LWC includes both cloud liquid and rain drops to give: re =  Ed3(LWC + RWC) 4πρwkNd 1/3 . (2.56) The effective radius is then limited to be in the range 4 μm to 30 μm. The remaining unknown is the number concentration of cloud droplets, Nd, described below. Note that even though RWC appears in the definition of re, the rain water content is not used elsewhere in the radiation scheme (in contrast to snow, which is added to the ice water content). The cloud droplet number concentration (cm−3) (due to activated CCN) is parametrized in the IFS following Martin et al. (1994). The formulation is based on extensive aircraft observations of stratocumulus clouds from both maritime and continental regions of the world. Over the sea the parametrization is: Nd = −1.15 × 10−3 Na2 + 0.963Na + 5.30, (2.57) and over land Nd = −2.10 × 10−4 Na2 + 0.568Na − 27.9, (2.58) where Nd is the cloud droplet number concentration (assumed to be equal to the number of activated cloud condensation nuclei) and Na is the aerosol number concentration in the size range 0.05 to 1.5 μm just below cloud base. Before IFS Cycle 31r2, the aerosol concentration was fixed in time with Na = 50 cm−3 over sea and Na = 900 cm−3 over land. From Cycle 31r2 onwards, a more complex parametrization of the aerosol number concentration was introduced to better represent the spatial and temporal variability of aerosol. The near surface aerosol mass concentration (qa) is calculated with a parametrization for the injection of aerosol from the surface into the atmosphere based on the 10-m wind speed. The formulation follows a parametrization derived for sea salt aerosol concentrations (Erickson et al., 1986; Genthon, 1992): qa = eaWs +b, (2.59) where Ws is the 10-m wind speed (m s−1) predicted by the model with coefficients depending on the wind speed; a = 0.16 and b = 1.45 for wind speeds less than 15 m s−1, and a = 0.13 and b = 1.89 otherwise. This is applied over both ocean and land, with qa limited to a value of 327 μg m−3 for wind speeds greater than 30 m s−1 over the ocean. The conversion of aerosol mass concentration qa to a number concentration Na then follows an empirical formulation as described by Boucher and Lohmann (1995), Lowenthal et al. (2004) and others: Na = 10c+d log10(qa ). (2.60) In the IFS, the coefficients are c = 1.2, d = 0.5 over ocean and c = 2.21, d = 0.3 over land. This Na is then used in Eqns. 2.57 and 2.58 to calculate the final cloud droplet number concentration. For very low wind speeds Nd is equal to 120 cm−3 over land and 40 cm−3 over sea, with increasing values as wind speed increases. Note that the above formula in Boucher and Lohmann (1995) and Lowenthal et al. (2004) is derived to give the number of activated cloud condensation nuclei rather than the aerosol number concentration, so there is some inconsistency in the use of both Eq. 2.60 and Eqns. 2.57 and 2.58. A major revision of the effective radius calculation is planned for a future update. IFS Documentation – Cy48r1 23 Chapter 2: Radiation (b) Ice cloud effective radius For ice clouds, the effective dimension is a function of temperature, T (K), and in-cloud ice water content, IWC (g m−3), based on Sun and Rikus (1999) with a revision from Sun (2001). Wherever there is ice cloud present, the snow water content is added to the ice to form a total ice water content as input to the scheme and therefore the same optical properties are currently assumed for both ice and snow. The ice particle effective diameter, Dice e (μm), is defined as Dice e = 1.2351 + 0.0105(T − T0) [a(IWC) + b(IWC)(T − 83.15)] , (2.61) where T0 = 273.15 K, a(IWC) = 45.8966IWC0.2214 and b(IWC) = 0.7957IWC0.2535. The effective diameter is then limited between a maximum of 155 μm and a latitude dependent minimum varying between 20 μm at the poles to 60 μm at the Equator, defined by 20 + 40 cos(latitude). The ice particle effective radius (μm) is then defined as rice e = 3√3 8 Dice e = 0.64952Dice e , (2.62) as suggested by Sun and Rikus (1999). (c) Cloud overlap The radiation scheme requires as input a profile of ‘overlap parameter’ α, on half-levels, defining the overlap of clouds in adjacent layers between the limits of ‘maximum’ (α = 1) and ‘random’ overlap (α = 0) overlap. Following the formulation of Hogan and Illingworth (2000), we compute this from decorrelation length, z0, such that α = exp(∆z/z0), where ∆z is the vertical distance between the midpoints of the adjacent layers. This is often referred to as ‘exponential-random’ overlap, since clouds in adjacent layers have an overlap parameter with an exponential dependence on separation, but if clouds in two layers are separated by a clear layer between them then their overlap is random. The IFS offers three means to specify ∆z. If namelist parameter NDECOLAT=0 then a fixed value of z0 = 2 km is used. Shonk et al. (2010) analysed data from cloud radar sites at a range of latitudes and proposed the following dependence on latitude (ϕ), in degrees, which is obtained with NDECOLAT=1: z0[km] = 2.899 − 0.0259|ϕ|. The larger decorrelation lengths in the tropics are consistent with the greater prevalence of convective clouds and the tendency for reduced wind shear compared to the mid-latitudes. The IFS default (NDECOLAT=2) uses a formulation that is consistent with the same radar observations, but has a smooth dependence with latitude across the equator: z0[km] = 0.75 + 2.149 cos2 ϕ. The radiation also needs to know how cloud sub-grid heterogeneities are correlated in the vertical. Hogan and Illingworth (2003) presented radar observations of ice clouds that suggested the vertical decorrelation length for cloud heterogeneities was half that for cloud boundaries, z0, so in the absence of more specific information for different cloud types, this value of one half is used for all cloud types. 2.4.3 Aerosols In the ECMWF forecast system from Cycles 26r3 to 43r1 the aerosol fields were provided by the Tegen et al. (1997) climatology. This is still used in the linearized version of the Morcrette (1991) radiation scheme used in data assimilation, but in other configurations of the IFS, the radiation scheme takes as input the mass mixing ratios of the 11 aerosol species from the Copernicus Atmospheric Monitoring System (CAMS). In CAMS air-quality forecasts, these values come directly from the 3D prognostic aerosol species after passing through the interpolation to the radiation grid. In other IFS configurations, we use a monthly mean 3D climatology derived from a CAMS aerosol reanalysis. As with trace gases 24 IFS Documentation – Cy48r1 Part IV: Physical Processes Table 2.2 Summary of the properties of the aerosol species used in the CAMS system, either in prognostic or climatological configuration, taken from Bozzo et al. (2020) except where updates have been made. The size distribution is described by the modal radius rmod and the geometric standard deviation σ; values are for dry aerosol apart from sea salt which is given for a relative humidity of 80%. Hydro- Radius Optical Density ρ rmod Aerosol type (code) Bin philic? range (μm) model (g cm−3) (μm) σ 1 Y 0.03–0.5 0.1992, 1.992† 1.9, 2.0† Sea salt (SU) 2 Y 0.5–5.0 OPAC 1.183 3 Y 5.0–20 1 N 0.03–0.55 Woodward or Composite‡ Desert dust (DU) 2 N 0.55–0.9 2.61 0.29 2.0 3 N 0.9–20 Organic matter (OM) 1 Y 0.005–20 OPAC See Bozzo et al. (2020) 2 N Black carbon (BC) 1 N* 0.005–0.5 OPAC 1.0 0.0118 2.0 2 N Sulfate (SU) – Y 0.005–20 GACP 1.76 0.0355 2.0 †Sea salt is described by a bimodal lognormal distribution with a 23.33 times larger number concentration in the small mode than the large. ‡Desert dust is treated using Woodward’s refractive indices when used in a climatological configuration, and ‘Composite’ refractive indices when used in a prognostic configuration. *Black carbon bin 1 is treated as hydrophilic in the CAMS aerosol scheme but as hydrophobic in the radiation scheme. (discussed in Section 2.4.4), the time interpolation is done linearly assuming the climatology data correspond to the 15th of each month. The details of the implementation together with the comparison against the Tegen et al. (1997) climatology and the impact on the model performance were discussed by Bozzo et al. (2020); a summary of the properties of the aerosol species is given in Table 2.2. See Section 11.7 for a visualization of the seasonal aerosol optical depth (AOD) of the CAMS climatology. 2.4.4 Greenhouse and other gases Molecular oxygen is assigned a global-mean volume mixing ratio of 0.20944. Naturally, water vapour is taken directly from the model’s prognostic variable. From Cycle 48r1, the model’s prognostic ozone variable is used in the radiation scheme. The formulation of the linear ozone chemistry scheme follows the description in Chapter 10, but from the same cycle has changed from using the Cariolle coefficients to using coefficients derived to best fit the CAMS ozone analyses; this is known as the Hybrid Linear Ozone (HLO) scheme. The five well-mixed greenhouse gases treated by the radiation scheme (CO2, CH4, N2O, CFC-11 and CFC-12) are derived from a monthly/zonal-mean climatology defined in the file greenhouse gas climatology 48r1.nc. The file name may be overridden by setting namelist parameter CGHGCLIMFILE. This file also contains a climatology of ozone, which may be used instead of the prognostic ozone by setting namelist parameter LEPO3RA to false (default true). The long-term trends in the five well-mixed greenhouse gases are treated by another file, greenhouse gas timeseries CMIP6 SSP370 CFC11equiv 47r1.nc, which implements CMIP6 Shared Socioeconomic Pathway SSP3-7.0, covering observed and projected surface concentrations for the years 0–2500 CE, and including a scaling-up of the concentration of CFC-11 to approximately represent the radiative forcing of all other minor greenhouse gases. There are a number of namelist parameters that allow the forcing to be changed in the IFS: IFS Documentation – Cy48r1 25 Chapter 2: Radiation • CGHGTIMESERIESFILE: specify an alternative file, in which case NGHGCMIP and NCP below are ignored. • NGHGCMIP: select the version of CMIP from 0 (CMIP3), 5 (CMIP5) or 6 (CMIP6, default). • NRCP: select the scenario from 1 (SSP1-2.6), 2 (SSP2-4.5), 3 (SSP3-7.0, the default) or 4 (SSP5-8.5). If NGHGCMIP=4 then the approximately equivalent scenario from CMIP5 is selected. Note that the values of NGHGCMIP and NRCP are used to select the timeseries file to read. • NCMIPFIXYR: specify a fixed year for the forcing (default −1, which means use the current year). 2.4.5 Surface albedo and emissivity The surface module of the IFS prepares the albedo and emissivity to be provided to the radiation scheme. Six bands are used to store albedo (λ < 0.25 μm, 0.25–0.44 μm, 0.44–0.69 μm, 0.69–1.19 μm, 1.19–2.38 μm and λ > 2.38 μm), although only sea ice uses different values for all six. Separate albedos are provided for diffuse and direct solar radiation. The albedos are computed separately for each tile, then averaged in proportion to the tile fraction. The snow-free land-surface albedo is taken from a monthly-mean albedo climatology derived from the MODIS dataset described by Schaaf et al. (2002). It is pre-averaged to the model grid, and stored in two spectral bands (split at 0.7 μm) with three components for each band referred to as the isotropic, volumetric and geometric albedos. These are used to represent the solar-zenith-angle dependence of the albedo to direct solar radiation as described by Schaaf et al. (2002). The 2-band direct and diffuse albedos computed from this climatology and the current solar zenith angle are then modified to add the contribution from snow according to the snow cover and the prognostic snow albedo (see Section 8.4.4). Since the latter is a broadband value in the absence of vegetation, it is split into the same two bands and modified according to the presence of high or low vegetation using the data of Moody et al. (2007). More information on the original data and plots of the monthly mean albedo are shown in Chapter 11. The ice sheets of Antarctica, Greenland and some glaciers are are designated as ‘permanent snow’; rather than using the prognostic snow albedo, they are assigned a fixed broadband value of 0.82. This is divided into two bands in the same way as the prognostic snow. The value may be overriden using the namelist parameter RALFMINPSN in namelist group NAEPHY. Over open water, the surface albedo for diffuse radiation is set as a wavelength-independent value of 0.06, while the albedo for direct parallel radiation is a fit to low-flying aircraft measurements over the ocean given by Taylor et al. (1996): αsp = 0.037 1.1μ1.4 0 + 0.15 . (2.63) For sea ice, monthly values based on Ebert and Curry (1993) albedos for the Arctic Ocean, for each of the six bands, are interpolated to the forecast time. The bare sea-ice albedo value in Ebert and Curry is taken as a representative value for summer, and the dry-snow albedo value is used for the winter months. Values for the Antarctic are shifted by six months. The thermal emissivity of the surface outside the 800–1250 cm−1 spectral region is assumed to be 0.99 everywhere. In the window region, the spectral emissivity is constant at 0.99 for open water, 0.98 for sea ice and 0.98 for snow tiles. For low and high vegetation and for shaded snow the emissivity depends on the water content in the top soil layer. Emissivity decreases linearly from 0.96 for soils at or above field capacity to 0.93 for soils at or below permanent wilting point. The same formulation is used for bare ground, except for desert areas (αsb > 0.3), where a value of 0.93 is used independently of the soil water content. 2.4.6 Solar irradiance and zenith angle The solar irradiance, S0, is the instantaneous flux of solar radiation incident at the Earth, and is computed from Total Solar Irradiance (TSI) and the Sun-Earth distance. TSI is the annual- mean solar irradiance and varies primarily with the 11-year solar cycle as described in the file total solar irradiance CMIP6 47r1.nc which includes an observational reconstruction from 1850 to 2014 CE and the CMIP6 recommended scenario from 2015 to 2300. This file may be overridden with 26 IFS Documentation – Cy48r1 Part IV: Physical Processes Table 2.3 Spectral distribution of the absorption by atmospheric gases in RRTMGLW. Note: CCl4 and CFC-22 (CHClF2) are presently not accounted for in the ECMWF model. Gases included Spectral intervals cm−1 Number of g-points Troposphere Stratosphere 10–350 8 H2O H2O 350–500 14 H2O H2O 500–630 16 H2O, CO2,N2O H2O, CO2,N2O 630–700 14 H2O, CO2 O3, CO2 700–820 16 H2O, CO2, O3, CCl4 O3, CO2, CCl4 820–980 8 H2O, CO2 CFC-11, CFC-12 980–1080 12 H2O, O3, CO2 O3, CO2 1080–1180 8 H2O, CFC-12, CFC22, O3, CFC-12, CFC-22, CO2,O3,N2O CO2,N2O 1180–1390 12 H2O, CH4, N2O CH4, N2O 1390–1480 6 H2O H2O 1480–1800 8 H2O, O2 H2O, O2 1800–2080 8 H2O 2080–2250 4 H2O, N2O O3 2250–2380 2 CO2 CO2 2380–2600 2 N2O, CO2 2600–3250 2 H2O, CH4 CH4 namelist parameter CSOLARIRRADIANCEFILE. Solar irradiance is then computed from TSI correcting for the changing Sun-Earth distance through the year. Equations to compute the annual variation of the solar declination δs and the difference between solar time and official time (the Equation of Time) can be found in Paltridge and Platt (1976). These equations are used to give the cosine of the solar angle at the ground. Because of the curvature of the Earth, the zenith angle is not quite constant along the path of a sun ray. Hence the correction applied to μa 0 to give an average μ0 for the atmosphere is μ0 = H a r (μa 0)2 + H a  2 + H a  − μa 0 . (2.64) where a is the earth radius and H is the atmospheric equivalent height. H/a is fixed at 0.001277. The cosine of solar zenith angle, μ0, is actually used twice in the model: (1) once per radiation timestep by the radiation scheme itself (including computing sea albedo in Eq. 2.63), and (2) at every intervening model timestep to scale the fluxes according to the instantaneous incoming top-of-atmosphere radiation (Eq. 2.108). In case (2), the value is computed as the average across the model timestep (using μ0 = 0 for parts of the model timestep when the sun is below the horizon), in order that when evolving temperature from timestep n to n + 1, the most accurate mean fluxes are used. In case (1), Hogan and Hirahara (2016) showed that for the most accurate fluxes and stratospheric temperatures, the value of μ0 used in the radiation scheme itself should be the average value only over the sunlit part of the radiation timestep. They provided the formula for computing the average μ0 over a particular time interval, and this is implemented in the IFS. 2.5 THE ECRAD RADIATION SCHEME 2.5.1 Overview A top-level scientific description of ecRad and the impact it had when introduced in the IFS was provided by Hogan and Bozzo (2018). The earlier memo by Hogan and Bozzo (2016) provided more of an architectural overview. Here we provide an overview of the mathematical formulation of ecRad, IFS Documentation – Cy48r1 27 Chapter 2: Radiation Table 2.4 Spectral distribution of the absorption by atmospheric gases in RRTMGSW. Gases included Spectral intervals cm−1 Number of g-points Troposphere Stratosphere 800–2600 12 H2O CO2 2600–3250 6 H2O, CH4 3250–4000 12 H2O, CO2 H2O, CO2 4000–4650 8 H2O, CH4 CH4 4650–5150 8 H2O, CO2 CO2 5150–6150 10 H2O, CH4 H2O, CH4 6150–7700 10 H2O, CO2 H2O, CO2 7700–8050 2 H2O, O2 O2 8050–12850 10 H2O 12850–16000 8 H2O, O2 O2 16000–22650 6 H2O 22650–29000 6 29000–38000 8 O3 O3 38000–50000 6 O3, O2 O3, O2 cloud properties optical properties clear−sky properties cloud optical Solver fluxes Cloud optics Gas optics Interpolate to model grid Interpolate to radiation grid ratios aerosol mixing ratios gas mixing surface optical properties Aerosol optics Figure 2.1 Schematic of the four main components of ecRad and the flow of data between them. and details of how to configure it from the IFS. The code for the ecRad radiation scheme within the IFS resides in radiation/module, while RRTMG gas-optics routines and the code for preparing inputs to ecRad resides in arpifs/phys radi. Hogan (2022) describes the full set of options for ecRad that may be configured via the RADIATION namelist group, while here we describe the more limited set of options intended for use in the IFS via the NAERAD namelist group; namelist options mentioned in the text refer to the NAERAD group. Fig. 2.1 depicts the four main components of ecRad, which are described in sections 2.5.2–2.5.5. 2.5.2 Gas optics The gas-optics module determines the spectral discretization used by the entire radiation scheme. Up to and including Cycle 48r1, the only gas-optics model available is RRTMG (Mlawer et al., 1997; Iacono et al., 2008), which has been extensively validated against observations, principally at the ARM Southern Great Plains site. It uses 16 spectral bands in the long-wave and 14 in the short-wave. It is a correlated 28 IFS Documentation – Cy48r1 Part IV: Physical Processes k-distribution scheme, which means that within each band the cumulative probability distribution of molar absorption coefficient, k, is discretized into g-points, each of which represents a similar range of absorptions but which may not be contiguous in wavelength space. Since an independent radiative transfer calculation is required for each g-point, the computational cost of the entire radiation scheme scales with the total number of g-points. The original RRTM scheme sub-divided each band into 16 ‘g-points’, but it was found that a more efficient scheme could be used by reducing the number of g-points in bands that either accounted for less radiative energy, or those with a weaker variation in gas absorption. This led to the development of RRTMG (the GCM-optimized version of RRTM), with the bands and number of g-points shown in Table 2.3 and Table 2.4; the total number of g-points in the long-wave and short-wave is 140 and 112, respectively. Note that there is some spectral overlap between the long-wave and the short-wave; this is because rather than defining these two regions by some arbitrary wavelength, we define the former as all radiation originating from emission by the surface or atmosphere, and the latter as all radiation originating from the sun. This separation makes the subsequent solvers more efficient since only the long-wave solver needs to treat atmospheric emission and only the short-wave solver needs to treat an additional ‘direct’ path of radiation through the atmosphere at the angle from the sun. The gas-optics module takes as input the mass-mixing ratios of H2O, O2, O3, CO2, CH4, N2O, CFC-11 and CFC-12. Some gas-optics models support additionally N2, CCl4 and/or HCFC-22, but these are not currently used in the IFS version of ecRad. The source of these mixing ratios, whether from prognostic variables, a climatology or a global constant, is described in Section 2.4.4. The module is also fed the temperature and pressure at half-levels (i.e. at the interfaces between layers), Ti+1/2 and pi+1, and the skin temperature Ts. Additionally it is provided with the solar irradiance, S0, which varies through the year according to the varying sun-earth distance, and through the 11-year solar cycle. From these it computes the following variables: • Short-wave absorption optical depth τa and Rayleigh scattering optical depth τs of each layer and g-point. These are converted to the extinction optical depth τ = τa + τs and the single-scattering albedo ω = τs/τ. • Long-wave absorption optical depth τ of each layer at each g-point, making use of the mid-point pressure and the pressure-weighted mean temperature of the layer. Rayleigh scattering is so weak in the long-wave that it may be neglected. • The solar irradiance associated with each g-point, calculated from the assumed solar spectrum which in the case of RRTMG is the Kurucz spectrum. • The Planck function at each half-level and long-wave g-point, Bi+1/2, i.e. the power emitted by a black body at temperature Ti+1/2 that is attributable the wavelengths associated with that g-point. Thus, if all g-points were summed we would expect to recover the total power emitted by a black body (σTi+1/2, where σ is the Stefan-Boltzmann constant). • The Planck function at at the surface in each long-wave g-point, Bs. Since version 1.5, the offline version of ecRad also supports the ECMWF Correlated K-Distribution (ecCKD) gas-optics scheme of Hogan and Matricardi (2022), in which the number of g-points is more flexible. This will be introduced into the IFS in a future cycle. 2.5.3 Aerosol optics Whether produced from an aerosol climatology or the CAMS prognostic aerosol package, the radiation scheme is passed profiles of aerosol mass mixing ratios for a number of different aerosol species. As part of the setup phase of ecRad, each aerosol species is mapped to pre-computed aerosol optical properties stored in a netCDF file1 These properties are computed using Mie theory as described by Bozzo et al. (2020) and pre-averaged to the RRTMG bands (16 in the long-wave and 14 in the short-wave); for aerosol species j they are (in both the long-wave and short-wave) the mass-extinction coefficient βj (in m2 kg−1), the single scattering albedo ωj and the asymmetry factor gj. Some of the species are hydrophilic, which means that their optical properties are stored in the file as a function of relative humidity in bins of width 1The latest aerosol optical property file aerosol ifs rrtm 48r1 v2.nc is available from the ecRad web site: https:// confluence.ecmwf.int/display/ECRAD/. IFS Documentation – Cy48r1 29 Chapter 2: Radiation 10% (except between 80 and 100% where they are stored in bins of width 5%). When optical properties are computed for these species, the nearest bin is selected according to the relative humidity. The aerosol optical depth in a layer is computed as the sum of the contribution from each species τ = ∆p g0 naer ∑ j qj βj, (2.65) where qj is the mass mixing ratio of aerosol species j, ∆p is the pressure difference between the top and base of the layer, and g0 is the acceleration due to gravity. The mean single scattering albedo and asymmetry factor are the weighted average of the contributions of the naer species, where ω is weighted by extinction coefficient and g is weighted by scattering coefficient: ω = ∑naer j qj βjωj ∑naer j qj βj ; (2.66) g = ∑naer j qj βjωj gj ∑naer j qj βjωj . (2.67) The scattering patterns of both aerosol and cloud particles can exhibit strong forward scattering, which is not well captured in the two-stream paradigm used in ecRad. Therefore, Joseph et al. (1976) proposed ‘delta scaling’ of the optical properties: a fraction f of the scattering distribution is approximated as an infinitesimally narrow delta function at a scattering angle of 0◦, essentially treating it as if it had not been scattered at all. They showed that a reasonable fit to more accurate calculations could be obtained with the assumption that f = g2, leading to the following transformation: τ′ = τ(1 − ωg2); (2.68) ω′ = ω 1 − g2 1 − ωg2 ; (2.69) g′ = g 1 + g , (2.70) where the delta-scaled quantities τ′, ω′ and g′ are used in all subsequent calculations. If long-wave aerosol scattering is turned off then only the optical depth of aerosols is required: the single- scattering albedo is implicitly set to zero. For the emission rate to be correct we require the absorption optical depth rather than extinction optical depth, and hence replace (2.65) with τ = ∆p g0 naer ∑ j qj βj(1 − ωj). (2.71) Once the combined, delta-scaled optical properties of aerosols are computed in each band, they are mapped to each g-point and combined with the existing optical properties of gases in the same way, i.e. optical depth is summed, single-scattering albedo is weighted by extinction coefficient (or equivalently by optical depth) and asymmetry factor is weighted by scattering coefficient.. Note that gases only absorb in the long-wave (so the ω of gases is zero) and their scattering is isotropic in the short-wave (so the g of gases is zero). 2.5.4 Cloud optics The optical properties of clouds in each layer are also described by τ, ω and g in the cloudy part of each layer and band. For liquid clouds we use the ‘SOCRATES’ method to parametrize these optical properties in each band as a function of effective radius using Eqs. 6–8 of Hogan and Bozzo (2016). This approach originates from the Met Office radiation scheme. For ice clouds we use the ‘Fu’ scheme, which again parametrizes the properties versus effective radius, combining Fu (1996) in the short-wave and Fu et al. (1998) in the long-wave. 30 IFS Documentation – Cy48r1 Part IV: Physical Processes As with aerosol, delta scaling (Eqs. 2.68–2.70) is applied, but the resulting optical properties are not immediately added to the clear-sky properties on g-points; rather they are passed to the solver which deals with sub-grid cloud heterogeneity. 2.5.5 Solver (a) Overview The solver takes as input the clear-sky (gas plus aerosol) optical properties of each layer at each g-point, and the cloud optical properties of the cloudy layers in each band, and computes profiles of upwelling and downwelling short-wave and long-wave fluxes at each half-level. Several solvers are available in ecRad, which differ in the way that sub-grid cloud is treated; they are controlled by the namelist parameters NLWSOLVER and NSWSOLVERS. A value of 0 indicates the Monte Carlo Independent Column Approximation (McICA) solver, the default in the IFS, while 3 indicates ‘Tripleclouds’. A value of 2 yields the ‘SPARTACUS’ solver, which can approximately represent the effects of 3D radiative transfer, but this is still an experimental code and is not perfectly stable in single precision. In addition to the choice of overarching solver, a detail of the long-wave solver that can be specified is how scattering is treated: the namelist option NLWSCATTERING may be 0 to indicate no scattering is represented, 2 to indicate that both cloud and aerosol scattering is represented, or 1 (the default) to indicate that cloud scattering is represented but not aerosol. As described by Hogan and Bozzo (2018), the latter is almost as fast as if scattering is not represented, yet achieves most of the accuracy of the full-scattering treatment. (b) McICA cloud generator Originally proposed by Pincus et al. (2003) and used by a large fraction of weather and climate models worldwide, the idea of McICA is to feed each g-point with a stochastically generated (i.e. Monte Carlo) profile of cloud properties, with each sub-column treated independently by an efficient radiative transfer algorithm that assumes the cloud to be horizontally uniform (i.e. the Independent Column Approximation). The profiles are generated such that for a large number of samples, the generated clouds will respect the input variables controlling sub-grid cloud geometry, specifically the profiles of cloud fraction, a, the fractional standard deviation of in-cloud water content, FSD, and the overlap parameter between each layer, α. By default FSD is set to 1 everywhere, but this value can be overridden with namelist parameter RCLOUD FRAC STD. The IFS originally used the cloud generator proposed by R¨ais¨anen et al. (2004), but with the introduction of ecRad this was changed to the scheme described by Hogan and Bozzo (2018), which is more efficient via the use of fewer random numbers, and introduces less stochastic noise in the fluxes via the exact enforcement of the total cloud cover, which is computed from the profiles of cloud fraction and overlap parameter as described in section 3.2 of Hogan and Bozzo (2018). If the gas-optics scheme specifies ng g-points are needed to represent either the long-wave or short-wave spectrum, then the cloud generator produces ng profiles of the cloud scaling factor O. This contains zeros for layers to be treated as cloud- free, while for cloudy layers the in-cloud-mean optical depth is multiplied by O before being merged with the clear-sky optical properties. Over a sufficiently large number of samples, the non-zero values of O in a layer have a mean of unity and a standard deviation of FSD. Note that inhomogeneity is enacted by scaling both cloud water content and optical depth up and down by the same factor, while leaving cloud effective radius unchanged. (c) Tripleclouds The concept of the Tripleclouds solver was originally described by Shonk and Hogan (2008): each cloudy layer is divided horizontally into three regions, one clear-sky, one optically thin cloudy region and one optically thick cloudy region. The optical depths and area fractions of the two cloudy regions are chosen such that the cloud has the same mean in-cloud optical depth as specified by the input data, and the same FSD. The fractional overlap of each region in one layer with each region in the layer below is computed in a way that is explicitly consistent with the specified overlap parameter. Shonk and Hogan IFS Documentation – Cy48r1 31 Chapter 2: Radiation (2008) described in detail how the two-stream radiative transfer equations could be solved subject to this geometry. (d) Long-wave layer properties The long-wave two-stream equations express the relationship between upwelling and downwelling fluxes (F↑ and F↓, in W m−2) for a particular g-point: − dF↑ dτ = −γ1 F↑ + γ2 F↓ + B; (2.72) dF↓ dτ = −γ1 F↓ + γ2 F↑ + B, (2.73) where τ is the optical depth of the atmosphere measured down from top-of-atmosphere for that g-point, treated as a vertical coordinate, and B is the black-body emission for the same g-point (in W m−2). The terms on the right-hand-side of each expression represent, respectively, loss of radiation by scattering into the other stream and absorption, gain of radiation by scattering from the other stream, and emission. The exchange terms are given by (e.g. Fu et al., 1997) γ1 = D [1 − (1 + g)ω/2] ; (2.74) γ2 = D(1 − g)ω/2, (2.75) where ω and g are the combined single-scattering albedo and asymmetry factor of the gases, aerosol and cloud in the volume, and D = 1.66 is the ‘diffusivity’ of the atmosphere (Elsasser, 1942), i.e. the approximate factor by which the atmosphere is more optically thick to diffuse long-wave radiation than it is to radiation propagating directly up or down. These equations are solved exactly in each layer under the assumption that the scattering properties are constant vertically through the layer, but that B varies linearly with τ across the layer. These solutions are expressed in terms of four properties. The reflectance and transmittance of the layer to diffuse incident radiation are given by Eqs. 25 and 26 of Meador and Weaver (1980): R = γ2G [1 − exp(−2kτ1)] ; (2.76) T = 2kG exp(−kτ1), (2.77) where τ1 is the optical depth of the layer and k = [(γ1 − γ2)(γ1 + γ2)]1/2 ; (2.78) G = [k + γ1 + (k − γ1) exp(−2kτ1)]−1 . (2.79) The other two layer properties are the rates of emission from the layer (in W m−2) in an upward and downward direction (e.g. Eqs. 5 and 12 of Stackhouse and Stephens, 1991): S↑ = (Bi−1/2 + C) − (Bi−1/2 − C)R − (Bi+1/2 + C)T; (2.80) S↓ = (Bi+1/2 − C) − (Bi+1/2 + C)R − (Bi−1/2 − C)T, (2.81) where Bi−1/2 and Bi+1/2 are the black-body emission rates at the top and bottom of the layer, and C = Bi+1/2 − Bi−1/2 τ(γ1 + γ2) . (2.82) (e) Long-wave adding method Given the four variables, Ri, Ti, S↑ i and S↓ i for each model layer i, we compute fluxes using the ‘Adding Method’, which successively adds the contribution of a new layer to the previous ones, accounting for multiple reflections between the layers being added. The calculations described here are repeated for each g-point. The surface which has an albedo of An+1/2 = 1 − ϵs (where ϵs is the surface emissivity and 32 IFS Documentation – Cy48r1 Part IV: Physical Processes n is the number of layers) and an upward emitted flux of En+1/2 = ϵs Bs. We then loop up through the layers of the atmosphere computing Ai−1/2, which is the albedo of the entire surface-atmosphere system below half-level i − 1/2, and Ei−1/2, which is the upwelling flux at half-level i − 1/2 originating from emission below that half-level (e.g. Eqs. 9 and 11 of Shonk and Hogan, 2008): Ai−1/2 = Ri + T2 i Ai+1/2 1 − Ri Ai+1/2 ; (2.83) Ei−1/2 = S↑ i + Ti Ei+1/2 + Ai+1/2S↓ i 1 − Ri Ai+1/2 ! . (2.84) This is repeated up to top-of-atmosphere (half-level 1/2). The denominator in each equation is the result of a geometric series representing multiple reflections between the new layer (with reflectance Ri) and the layers below (with reflectance Ai+1/2). At top-of-atmosphere we have no downwelling long-wave radiation (F↓ 1/2 = 0) and the upwelling is equal to the total emission (F↑ 1/2 = E1/2). We now pass back down through the atmosphere computing up- and down-welling fluxes at each half-level (e.g. Eqs. 12 and 14 of Shonk and Hogan, 2008): F↓ i+1/2 = Ti + F↓ i−1/2 + Ri Ei+1/2 + S↓ i 1 − Ri Ai+1/2 ; (2.85) F↑ i+1/2 = Ai+1/2 F↓ i+1 + Ei+1/2. (2.86) This is repeated down to the surface, after which the fluxes from each g-point are summed to obtain the broadband fluxes. The upward-then-downward passes through the atmosphere are mathematically equivalent to solving a tridiagonal system of equations by Gaussian elimination followed by back- substitution, but with the benefit that the intermediate terms always have a physical meaning. The method described here allows for the possibility of scattering in all layers. If long-wave aerosol scattering is turned off then the reflectance Ri in clear-sky becomes zero and the computation may be simplified, as described by Hogan and Bozzo (2018). (f ) Short-wave layer properties The short-wave two-stream equations are as in the long-wave except that they include an additional ‘direct’ or ‘unscattered’ solar flux, F0 (defined as the power into a plane perpendicular to the solar direction), which acts as a source term for the diffuse fluxes: dF0 dτ = F0 μ0 ; (2.87) − dF↑ dτ = −γ1 F↑ + γ2 Fdiff↓ + ωγ3 F0; (2.88) dFdiff↓ dτ = −γ1 Fdiff↓ + γ2 F↑ + ωγ4 F0, (2.89) where μ0 is the cosine of the solar zenith angle, and we use expressions for the γ coefficients originating from Zdunkowski et al. (1980): γ1 = 2 − ω 5 + 3g 4 ; (2.90) γ2 = 3 4 ω(1 − g); (2.91) γ3 = 1 2 − 3 4 μ0 g; (2.92) γ4 = 1 − γ3. (2.93) In (2.88) and (2.89) we use Fdiff↓ to indicate the diffuse component of the downwelling, to distinguish from the total downwelling F↓ = Fdiff↓ + μ0 F0 (whereas in the long-wave there is no direct beam so IFS Documentation – Cy48r1 33 Chapter 2: Radiation F↓ = Fdiff↓). As in the long-wave case we seek exact solutions to these equations for individual layers under the assumption that the scattering properties are constant vertically through the layer. The diffuse reflectance, R, and transmittance, T, are already given by (2.76) and (2.77). The unscattered transmittance to direct solar radiance is simply the solution to (2.87): T0 = exp(−τ1/μ0). (2.94) We also require the reflectance to direct incoming radiation, Rdir, and the equivalent transmittance Tdir, but note that the latter is the fraction of direct flux at the top of the layer that emerges as diffuse radiation at the base of the layer (adapted from Meador and Weaver, 1980): Rdir = Gdirh (1 − kμ0)(α2 + kγ3) − (1 + kμ0)(α2 − kγ3) exp(−2kτ1) − 2k exp(−kτ1)(γ4 − α2μ0)T0i ; (2.95) Tdir = Gdirn 2k exp(−kτ1)(γ4 + α1μ0) − T0 [(1 + kμ0)(α1 + kγ4) − (1 − kμ0)(α1 − kγ4) exp(−2kτ1)] o , (2.96) where Gdir = μ0ω 1 − k2μ2 0 G; (2.97) α1 = γ1γ4 + γ2γ3; (2.98) α2 = γ1γ3 + γ2γ4. (2.99) The denominator of (2.97) can in principle go to zero, and inaccurate results can occur in the vicinity of zero. This is avoided by detecting when kμ0 ≃ 1 and then perturbing μ0 very slightly to avoid the singularity. (g) Short-wave adding method Given the five variables Ri, Ti, Rdir i , Tdir i and T0 i for each model layer i, we compute fluxes using the short-wave version of the Adding Method. As in the long-wave, we propagate the diffuse albedo Ai−1/2 from half-level to half-level up through the atmosphere, but rather than also propagating an emitted flux Ei−1/2 which would be zero at solar wavelengths, we propagate the direct albedo Di−1/2, which is the albedo of the entire surface-atmosphere below half-level i − 1/2 to incoming direct solar radiation. The radiation scheme is provided with diffuse surface albedo α and direct surface albedo αdir. The surface values of the propagated variables are therefore simply An+1/2 = α and Dn+1/2 = μ0αdir, where the factor of μ0 is necessary because the direct flux F0 is defined to be into a plane perpendicular to the solar direction whereas the reflected flux F↑ is for a horizontal plane. Diffuse albedo is propagated exactly as in (2.83), whereas direct albedo uses Di−1/2 = Rdir i + Ti+1/2 T0 i Di+1/2 + Tdir i Ai+1/2 1 − Ri Ai+1/2 . (2.100) At top-of-atmosphere the downwelling direct solar radiation is the solar irradiance for the specific g- point (F0 1/2 = S0) and there is no downwelling diffuse solar radiation (Fdiff↓ 1/2 = 0). We pass back down through the atmosphere computing fluxes using F0 i+1/2 = T0 i F0 i−1/2; (2.101) Fdiff↓ i+1/2 = Ti Fdiff↓ i−1/2 + F0 i−1/2(T0 i Di+1/2 Ri + Tdir i ) 1 − Ri Ai+1/2 ; (2.102) F↑ i+1/2 = F0 i+1/2 Di+1/2 + Fdiff↓ i+1/2 Ai+1/2. (2.103) As in the long-wave, the fluxes from each g-point are summed to obtain the broadband values. 34 IFS Documentation – Cy48r1 Part IV: Physical Processes 2.6 SPATIAL AND TEMPORAL INTERPOLATION As stated in the introduction, the radiation scheme is called intermittently in time and on a lower- resolution spatial grid. The grid is described in Section 2.6.1. The fluxes are interpolated back to the full model grid and are then available at each intervening model timestep between calls to the radiation scheme. The time between calls to the radiation scheme may be set by namelist parameter NRADFR: a positive number indicates the number of model timesteps between calls, whereas a negative number indicates the number of hours. The default is hourly. At each model gridbox and timestep, surface fluxes are fed into the surface scheme to evolve surface temperature, and atmospheric radiative heating rates are computed as the divergence of net radiation fluxes Fnet so that  ∂T ∂t  rad = − g0 cp ∂Fnet ∂p (2.104) where cp is set to the specific heat at constant pressure of dry air cpdry . Ho