Radiation sensitivity tests of the HARMONIE 37h1 NWP model

. When shortwave (SW) radiation ﬂuxes modelled with a numerical weather prediction (NWP) model or climate model do not match observed SW ﬂuxes it can be challenging to ﬁnd the cause of the blue differences. Several elements in the model affect SW ﬂuxes. This necessitates individual testing of each of the physical processes in the NWP model. Here we present a focussed study of the SW radiation schemes in the H IRLAM (HIgh Resolution Limited Area Model) A ladin R egional M esoscale O perational N WP I n E urope ( HARMONIE ) model, which is the primary NWP model used and developed by several National Weather Services in Europe. Detailed calculations have been made with the DISORT model run in the libRadtran framework, which is a collection of state-of-the-art radiative transfer software and datasets. These are used to test the NWP radiation calculations. Both models are given the same atmospheric properties as input. We also perform a separate test of cloud liquid optical property parametrizations with Mie calculations. This leads us to introduce a new parametrization for calculating these properties. In addition, we show that the results of a simpler radiation scheme, introduced into HAR-MONIE, compare well with those of the comprehensive default parametrizations. The methodology applied here may be used for testing radiation schemes in other NWP or climate models.


Introduction
One of the main problems in comparing numerical weather prediction (NWP) models to observations is that there is generally a long chain of elements in the models that leads to the given results.This is also true for modelled solar radiation fluxes.There can be multiple causes of biases in SW Correspondence to: Kristian Pagh Nielsen (kpn@dmi.dk)fluxes: The release of precipitation could be parametrized incorrectly; the liquid or ice water load in the clouds or atmospheric aerosols could be wrong; the surface radiative properties could be described incorrectly; the amount of soil water available for evaporation could be wrong; the fluxes of turbulent energy could by parametrized incorrectly; the radiative transfer calculations could be wrong; etc.Some of these causes are not directly related to the SW fluxes but affect the fluxes indirectly through, for example, the formation and lifetime of clouds and fog.Very often multiple solutions exist for correcting an observed SW bias, where several of the mentioned parametrizations or assumptions can be tuned to correct a given bias.However, this equifinality needs to be resolved.One way of doing this is by testing each of physical processes in the NWP model separately.Here we present a focused study of one component of the HIRLAM (HIgh Resolution Limited Area Model) Aladin Regional Mesoscale Operational NWP In Europe (HARMONIE) modelling system, i.e. the SW radiation scheme.Our particular focus on this component is because it has not previously been studied in the HARMONIE modelling system.
Correctly modelling the SW radiative transfer in an NWP model depends firstly on correctly describing the state of the atmosphere and the surface.In particular cloud liquid water and cloud ice concentrations are important, but cloud droplet number concentration, cloud ice particle concentration, cloud ice particle shape and size distributions, water vapour, aerosols, ozone, surface reflectance, Sun-Earth distance, air pressure, altitude and position of the Sun relative to the given topography can also significantly affect the net and downward global radiation (Sagan and Pollack, 1967;Hansen and Travis, 1974;Pierluissi and Peng, 1985;Shettle, 1989;Hu and Stamnes, 1993;Fu, 1996;Hess et al., 1998;Thomas and Stamnes, 2002;Reddy et al., 2005;Forster et al., 2007;Myhre et al., 2007;Senkova et al., 2007;Kahnert et al., 2008;L'Ecuyer et al., 2008;Manners et al., 2012).Secondly, when the atmospheric state is known, the optical properties of particles and gases in the atmosphere need to be correctly calculated.Thirdly, a sufficiently accurate radiative transfer scheme is needed to estimate the solar radiation fluxes.
Here we focus on the second and third steps described above, by testing various radiation parametrizations, available in the HARMONIE NWP model, against the accurate one-dimensional radiative transfer model, DISORT (Stamnes et al., 1988(Stamnes et al., , 2000) ) in prescribed atmospheric conditions.Since not all of the radiation schemes available in HAR-MONIE perform the second and third steps separately, we primarily present tests of these two steps together.DISORT is run within the framework of the libRadtran library (Mayer and Kylling, 2005). 1 It solves the radiative transfer (the third step), while various parametrizations in libRadatran are used to calculate the optical properties (the second step).In addition to using libRadtran/DISORT we also perform separate Mie calculations with the solver of Wiscombe (1980) in order to perform direct tests of the liquid cloud optical property calcuations in HARMONIE.Details on these tests can be found in: "Supplement 1: Mie calculations" (reference to the supplements).
The HARMONIE NWP system combines elements from the global IFS/Arpege model (Déqué et al., 1994) with the ALADIN nonhydrostatic dynamics (Bénard et al., 2010) and physical parametrizations of AROME (Seity et al., 2011) and HIRLAM (Undén et al., 2002).Here IFS is an abbreviation of Integrated Forecast System, ALADIN is an abbreviation of Aire Limitée Adaptation dynamique Développement In-terNational and AROME is an abbreviation of Applications of Research to Operations at MEsoscale.In the current version (37h1), the default SW radiation scheme is based on the radiation transfer code in the Integrated Forecast System (IFS cycle 25R1, European Centre for Medium-Range Weather Forecast implementation in 2002), see ECMWF (2012) and Mascart and Bougeault (2011).This is a comprehensive, but also computationally demanding, scheme with several spectral bands both in the solar and terrestrial radiation ranges.In the following, we refer to this as the IFS radiation scheme.Here we test this scheme and also the alternative HIRLAM radiation scheme, further denoted as hlradia, which is based on Savijärvi (1990).Hlradia is a simple and computationally fast radiation parametrization, where the cloud transmittance and absorptance are calculated using the parametrization of Wyser et al. (1999).Oreopoulos et al. (2012) recently compared a large set of SW and longwave radiation schemes that are currently used in atmospheric models.In their tests, realistic atmospheric profiles from the Atmospheric Radiation Measurement programme were used.Here we use idealised cases which are not necessarily realistic.The advantage of using idealised cases is that the full range of sensitivity can be tested for each of the relevant variables individually.Fixed atmospheric 1 www.libradtran.orgIn practice, a simple radiation scheme applied at each time-step (currently 1 minute) of the model integration requires as much computation time as a complicated scheme called a few times an hour.The accuracy requirements for radiation parametrizations in long-term climate simulations on relatively coarse spatial and temporal scales differ from those of the kilometre-scale short-term NWP models run in continuous data assimilation forecast cycles.For example, the importance of quickly changing cloud-radiation interactions may be greater in the latter while the former requires maximum accuracy in handling the details of gaseous transmission.It is thus of interest to understand and compare the accuracy of the simple and complex schemes against a good reference.Such a reference for testing the HAR-MONIE SW radiation parametrizations was obtained by applying the DISORT model.The computations done by the detailed radiative transfer model also provide valuable information about the sensitivity of SW radiation fluxes to a wide range of various atmospheric properties.

Methods
Different SW radiation schemes currently available in the HARMONIE framework were compared with DISORT results for 11 experiments.The experiments can be split into three groups: clear sky experiments 1-3, liquid cloud experiments 4-7, and ice cloud experiments 8-11.In each of the experiments the sensitivity to a different variable was studied as outlined in Table 1.

DISORT
The DISORT algorithm was run for the full SW spectrum of the Kato et al. (1999) correlated-k algorithm with absorption coefficients from the HITRAN 2000 database and an angular discretisation of 30 streams.For experiment 2 we used the pseudo-spherical version of DISORT (Dahlback and Stamnes, 1991) in order to account for the sphericity of Earth's atmosphere at higher solar zenith angles.
We used the Hu and Stamnes (1993) parametrization that is accurate to within 1% with respect to full spectrum SW fluxes to calculate the liquid cloud optical properties.We used the Fu (1996) parametrization to calculate the ice optical properties for DISORT.

Default IFS radiation scheme in HARMONIE
The default IFS radiation scheme in HARMONIE has six SW spectral bands: 3 bands in the ultraviolet and visible spectral ranges and 3 bands in the solar infrared spectral range.Specifically, these six spectral bands are defined by the wavelengths: 0.185-0.25-0.44-0.69-1.19-2.38-4.00µm (Mascart and Bougeault, 2011).At the top of the atmosphere each of these bands have 0.19%, 13.57%, 32.21%, 32.62%, 18.06% and 3.35% of the SW flux, respectively.The top of atmosphere SW flux is calculated using the formulae of Paltridge and Platt (1976).
The radiative transfer calculations were done using the delta-Eddington approximation (Joseph et al., 1976;Fouquart and Bonnel, 1980).In HARMONIE version 37h1, the default cloud liquid optical properties for the IFS scheme are calculated using the Fouquart (1987) parametrization.It is also possible to use the Slingo (1989) parametrization.The cloud ice optical properties are calculated using the Fu and Liou (1993) parametrization.Alternatively, the Fu (1996) parametrization can be used.Here we test each of these four cloud optical property options.
Additionally, we propose a new cloud liquid optical property scheme (hereafter, referred to as the Nielsen scheme) based on the Mie solution to Maxwell's equations (Mie, 1908;Wiscombe, 1980) for each of the 6 SW spectral bands used in the HARMONIE IFS scheme.The new parametrization is described in Eqs.3-5 and in Table 2. Details on the derivation of this new parametrization can be found in: "Supplement 1: Mie calculations" (reference to the supplements).In the IFS radiative transfer calculations the key input variable is the scaled optical depth τ * i , which is a product of the optical depth τ i and the delta-Eddington scaling factor (1 − ω i g 2 i ) (Joseph et al., 1976): where, according to the suggested new scheme, Here, i denotes the wavelength band index, ω i is the single scattering albedo, g i is the asymmetry factor, β i is the mass extinction coefficient in m 2 /g and C is the mass load in g/m 2 , r e, liq is the cloud droplet effective radius, and a i , b i , c i , d i , e i , f i , h i and j i are coefficients as detailed in Table 2.
As can be seen in Eq. 2, the asymmetry factor g i is important for the scaling factor.For this reason we have made a more comprehensive empirical formula for g i in Eq. 5 than those used in the Fouquart (1987) and Slingo (1989) parametrizations, where g i is given as constants in the former and as linear functions of r e in the latter.The parametrization for wavelength band 1 (0.185-0.25 µm) is irrelevant, as virtually no SW irradiance at these wavelengths reaches the levels of the lower atmosphere, where clouds occur.This is why it is omitted in Table 2.

HIRLAM radiation scheme (hlradia)
A slightly modified version of the HIRLAM radiation scheme hlradia as compared to the original version (Undén et al., 2002)) was used in our tests in the HARMONIE framework.In hlradia one SW spectral band is considered and both the clear sky and cloudy sky transmittances and absorptances are parametrized in order to make the scheme very fast (Savijärvi, 1990;Wyser et al., 1999).The impact of ozone, oxygen and carbon dioxide, as well as aerosols, on SW irradiance is assumed to be constant over time and space.Cloud transmittance and absorptance are calculated based on empirically derived functions of cloud condensate content and cloud particle effective radia.They are integrated vertically from the top of the atmosphere to each level, combining cloud liquid and ice particles.Thus within hlradia, the calculation of cloud optical properties and the radiative transfer calculation are not formally separable.
3 Experiments and initial atmospheric conditions MUSC, the single column version of HARMONIE, based on Malardel et al. (2006), was used as a framework to provide the radiation schemes with the atmospheric and surface input data.A vertical resolution of 41 full model levels (in pressure-based hybrid coordinates) between the surface and an elevation of 100 km was used (Undén, 2010).Hybrid coordinates follow the topography at the lowest levels and converge towards isobaric levels higher in the atmosphere (Holton, 1992).The same model levels were used in the DISORT calculations.Only output from the first time-step of MUSC was used in order to exclude the influence of an evolving atmospheric state on the results.
The reference (initial) definitions for the HARMONIE and DISORT experiments are: -Day of year = 79, i.e. the 20 t h of March or vernal equinox -Altitude = 0.0 km above sea level -Solar zenith angle = 56 • -Surface albedo = 0.18 -Clear sky -"Mid-latitude summer" atmospheric profile from Anderson et al. (1986), including geopotential heights, temperature, air densities, mixing ratios of water vapour, ozone, oxygen and carbon dioxide from the surface to the top of the atmosphere assumed to be at 120 km.
Water vapour input data were interpolated to the model levels, keeping the vertical distribution unchanged when modifying the vertically integrated load in the experiments.In experiment 1 (integrated water vapour) we also tested the effect of using the "mid-latitude winter" atmospheric profile (Anderson et al., 1986).
In order to make the sensitivity experiments as simple as possible, all clouds were assumed to be homogeneous, to be situated between 1 km and 2 km above the surface, and to have prescribed properties.The prescribed properties comprise vertically integrated specific cloud liquid and ice content, a constant liquid droplet effective radius and a constant ice crystal equivalent radius (Table 1).Such clouds may be unrealistic and in contradiction to the assumed water vapour distribution.However, our aim here is not to test realistic experiments, but rather to make focussed tests of the radiation schemes in HARMONIE.Using idealised clouds also enables us to test a wider variety of cloud types than if we, for instance, had used cloud data from measurement campaigns.
Aerosols are described differently in each of the three radiation schemes: in libRadtran/DISORT the Shettle (1989) rural aerosol profile is used, in the HARMONIE IFS scheme an aerosol climatology (Tanré et al., 1984;Mascart and Bougeault, 2011) is applied and in hlradia aerosols are represented by constant coefficients (Savijärvi, 1990).As these cannot be compared in detail, i.e. for various aerosol types and concentrations, we have excluded aerosols from this study.

Results and discussion
In the following sub-sections the total downward SW irradiance at the surface on a horizontal plane, i.e. the global radiation [W/m 2 ], is presented for each of the 11 radiation sensitivity experiments (Table 1).In addition, net fluxes on model levels are presented for a select number of the experiments.Net fluxes, defined as the downward SW fluxes minus the upward SW fluxes, are presented as % differences relative to the accurate DISORT model.Clear sky experiments are discussed in section 4.1 and cloud liquid and cloud ice experiments are discussed in sections 4.2 and 4.3, respectively.In these experiments when the cloud condensate is > 0 kg/m 2 the fractional cloud cover was set to 1.

Clear sky experiments
Three clear sky experiments were carried out to test the SW radiation sensitivity to water vapour (experiment 1), solar zenith angle (experiment 2) and surface albedo (experiment 3).

Experiment 1: Water vapour
In the water vapour experiment (experiment 1), the agreement between the IFS radiation scheme and DISORT is very good for the highest water vapour loads but has a negative bias of 2% or 15 W/m 2 at the lowest water vapour load.Hlradia has a positive bias of 3.5% or 20 W/m 2 relative to DISORT at the highest water vapour load, but lower bias at the lower water vapour loads (Fig. 1).When the "midlatitude summer" atmospheric profile is replaced with the "mid-latitude winter" atmospheric profile a decrease of approximately 1% in global radiation is seen in all three models in Fig. 1.The primary reason for this decrease is pressure broadening of the water vapour absorption lines (Thomas and Stamnes, 2002).This shows that the vertical distribution of the water vapour in the atmosphere is also important for the clear sky transmittance.Since the "mid-latitude winter" atmosphere cannot contain the highest water vapour loads there are fewer points for this in Fig. 1.The height level net flux deviation of IFS from DISORT (Fig. 2) varies from approximately 0.0% at the top of the atmosphere to the negative biases also seen in Fig. 1 for the lowest water vapour loads at the surface.The hlradia net flux differences are between 0.0% and +3.0% (Fig. 3).In all of the following experiments the "mid-latitude summer" profile was used.
Gases other than water vapour affect the SW radiation.The most important of these is ozone.We have tested the SW sensitivity to ozone by running DISORT.Here the ozone cross section parametrization of Bass and Paur (1985) was used.When ozone was reduced from 500 Dobson Units (DU) to 100 DU, the global radiation at the surface increased by 2.3%.Since 500 DU and 100 DU are extreme values of the ozone layer thickness, we find that there is little error in using climatological ozone data (as is done in the IFS radiation scheme), or in using constant ozone SW absorption (as is done in hlradia).In an hlradia experiment, a complete re- The variations in carbon dioxide and oxygen are practically irrelevant for SW radiation in NWP modelling, because they only affect the extreme fringes of the IR and UV parts of the solar spectrum, respectively.As an example we tested the sensitivity of SW radiation to carbon dioxide concentration with DISORT.The carbon dioxide concentration is set to 330 ppm in the "mid-latitude summer" profile from Anderson et al. (1986).In both IFS and hlradia, this is set to be 353 ppm.Increasing the carbon dioxide concentration in DISORT from 330 ppm to 353 ppm causes a relative decrease of 0.01% in the surface downward global SW radiation.The relative impact of changing the oxygen concentation on the surface downward global SW radiation is even less.In an hlradia experiment, the effect of doubling the CO 2 concentration on the SW radiation fluxes was negligible.

Experiment 2: Solar zenith angle
The solar zenith angle (SZA) experiment (experiment 2, Fig. 4) shows agreement between the IFS global radiation and DISORT of better than 1.0% for most of the SZAs.Hlradia shows a systematic overestimation of +2% to +5%.

Experiment 3: Surface albedo -clear sky case
In experiment 3 the sensitivity of varying the surface albedo in clear sky conditions was tested (Fig. 5).The relative difference for the IFS radiation scheme compared to DISORT is less than 1% for all surface albedos.For the hlradia radiation scheme the relative differences range from +3% to +4.5%, increasing with increasing albedo.The differences in net fluxes as a function of height in the atmosphere are shown in Figs. 6  and 7 for the IFS and the hlradia radiation schemes, respectively, relative to DISORT.These differences are mostly less than 2% for the IFS radiation scheme and less than +5% for the hlradia scheme but they decrease for the higher surface albedos.

Liquid cloud experiments
Four liquid cloud experiments were run to test the SW radiation sensitivity to cloud water load, cloud drop effective radius (r e, liq ) and surface albedo.Tests of the surface albedo sensitivity were run for both the case of a cloud of medium thickness (cloud water load 0.1 kg/m 2 ) and the case of a thick cloud (cloud water load 1.0 kg/m 2 ).

Experiment 4: Cloud water load
Experiment 4 was run to test the cloud water load sensitivity.A fixed r e, liq of 10 µm was used, which is a typical value for stratocumulus clouds (e.g.Martin et al., 1994).Before being used for SW radiation calculations the cloud water load is multiplied by a so-called cloud SW inhomogeneity factor of 0.7 in the default SW scheme of HARMONIE 37h1.In hlradia, a SW inhomogeneity factor of 0.8 is used by default.After several tests, it became clear that the SW inhomogenity  factor is very significant.We therefore performed an additional test showing the effect of changing the inhomogeneity factor (Fig. 8 and Table 3).
For cloud water loads of 0.05 kg/m 2 to 0.1 kg/m 2 the transmitted global radiation is almost 50 W/m 2 higher with the SW inhomogeneity factor set to 0.7, compared to when it is set to 1.0, for the case where the default Harmonie radiation settings are used i.e.IFS with the Fouquart cloud optical property parametrization (Fouquart, 1987).Similarly, when the hlradia scheme is used with a SW inhomogeneity factor of 1.0 rather than 0.8, the global radiation is more than 30 W/m 2 less.Note that in our experiments, where cloud water load is given as an input value, application of an inhomogeneity factor means reduction of this predefined input by 20 % (hlradia) or 30 % (IFS).
Overall, when the cloud SW inhomogeneity factor is set to 1.0, the results are closer to those from DISORT.This is not surprising, as the DISORT tests were done for horizontally homogeneous clouds.In both the Meso-NH model (Mascart and Bougeault, 2011) and in the most recent cycles of the global IFS model (ECMWF, 2012) the cloud SW inhomogeneity factor is set to 1.0.Some meteorological institutes also run hlradia (within the operational HIRLAM) with this setting.In the following cloud experiments (4-11) we set the cloud inhomogeneity factor to 1.0.In Sect.4.2.2 we discuss the cloud SW inhomogeneity at a more general level.
The discrepancies between using the IFS radiation scheme with the HARMONIE 37h1 default cloud liquid optical property scheme (Fouquart, 1987), hlradia and the accurate DIS-ORT model are shown in Fig. 9.In this figure the results obtained when using two alternative schemes are also shown.The first of these schemes is the Slingo (1989) scheme that is optionally available in HARMONIE 37h1, while the second is the Nielsen scheme as described by Eqs.(3-5) and Table 2.In Table 4 the results are shown in the form of absolute differences compared to the DISORT calculations.
It can be seen that the default IFS radiation scheme in HARMONIE, with the Fouquart (1987) liquid cloud optical Table 3. Absolute differences (unit: W/m 2 ) for experiment 4. IFS-F refers to the IFS radiation scheme using the Fouquart (1987)   Table 4. Absolute differences (unit: W/m 2 ) for experiment 4. IFS-F, IFS-S and IFS-N refer to the IFS radiation scheme using the Fouquart (1987), the Slingo (1989)  property scheme, shows a positive bias of up to +41 W/m 2 (+10%) for clouds with a 0.02 kg/m 2 cloud water load but shows a smaller bias for more water-rich clouds.For the alternative Slingo (1989) scheme the results show a negative bias of up to -29 W/m 2 (-14%) for clouds with a 0.1 kg/m 2 cloud water load.The Nielsen scheme has the lowest overall bias in this experiment, with a maximum bias of +12 W/m 2 (+3%) for a cloud water load of 0.02 kg/m 2 .There is a high positive bias of up to +45 W/m 2 (+10%) for a cloud water load of 0.01 kg/m 2 for the hlradia scheme.For cloud water loads of 0.1 W/m 2 and higher, the hlradia scheme has practically no bias.Here the positive bias of the hlradia scheme in clear sky conditions (Figs. 1,3,4,5 & 7) should be kept in mind.If this were to be corrected, the hlradia results would be around 4% lower overall in this experiment.
In Figs.10-13 the relative differences between the net fluxes as a function of height are shown for each of the four parametrizations tested.As for the global radiation, the net fluxes mostly have positive biases both below and above the clouds when the Fouquart parametrization is used.An increasingly negative bias is seen below increasingly thicker clouds (Fig. 10).Also for the Slingo and Nielsen parametrizations increasingly negative biases are seen below increasingly thicker clouds.The results for the hlradia scheme do not display this pattern (Fig. 13).Instead, a clear negative bias (< -20%) is seen in the net flux at the top of the thicker clouds, but not below and above these.Such a pattern can be explained by the cloud absorptance being too large and the cloud reflectance being too low in the hlradia scheme.

General discussion on cloud inhomogeneity
In an NWP model grid box clouds may be inhomogeneously distributed.It can easily be shown that representing clouds as a linear average covering a certain fraction of the grid box does not produce the same SW transmittance and reflectance as when the clouds are resolved.This is because the cloud transmittance and reflectance do not vary linearly with cloud optical depth.On the other hand, in cases where clouds are homogeneously distributed within a gridbox the linear average is correct.Therefore, it is incorrect to multiply clouds in all grid boxes by the same cloud inhomogeneity factor, as is currently done in HARMONIE.To deal with these issues properly, and also to deal with the effects of 3D cloud radiation transfer, a more sophisticated scheme like the one recently proposed by Hogan and Shonk (2013) is required.

Experiment 5: Cloud drop effective radius
Experiment 5 was run to study the sensitivity of the schemes to varying r e, liq , keeping the integrated cloud water load at 0.1 kg/m 2 .The global radiation results are shown in Fig. 14 and in Tables 5 and 6.The net irradiance results relative to DISORT are shown in Figs.15-18.
The new Nielsen parametrization clearly performs best compared to the DISORT calculations (Fig. 14 and Tables 5-6).The Slingo parametrization performs considerably worse when r e, liq is varied, which is in contrast to experiment 4 where the liquid cloud water load was varied (Fig. 9).This parametrization was designed for a limited cloud droplet size range of 4.2-16.6µm (Slingo, 1989).Therefore, it is not surprising that its performance deteriorates for large cloud droplet sizes.Previously, Dobbie et al. (1999) have shown that the Slingo parametrization does not agree with calculations based on Mie theory for large cloud droplets.The Fouquart parametrization has a consistently positive bias Table 5.Relative differences for experiment 5. IFS-F, IFS-S and IFS-N refer to the IFS radiation scheme using the Fouquart (1987), the Slingo (1989) and the Nielsen liquid cloud optical property parametrizations, respectively.as a function of r e, liq relative to DISORT (+20 W/m 2 to +34 W/m 2 i.e. +8% to +22%).The hlradia parametrization performs very well for the small cloud droplets but has a positive bias for cloud droplets which are larger than 15 µm.In Figs.15-18 the corresponding results for the net flux differences are shown.Again, the bias in the net fluxes is      consistently positive for all height levels when the Fouquart parametrization is used (Fig. 15).The biases of the Slingo (Fig. 16) and Nielsen (Fig. 17) schemes are mainly below the cloud rather than above.For the Nielsen scheme the negative bias increases with decreasing r e, liq .As in Fig. 14, the results for the hlradia scheme are also very good (Fig. 18) for the atmospheric net fluxes in this experiment.

Experiments 6 and 7: Surface albedo under liquid cloud conditions
Experiment 6 was run for a fixed integrated cloud water load of 0.1 kg/m 2 and a fixed r e, liq of 10 µm.The results of varying the surface albedo can be seen in Fig. 19, which shows that the three cloud optical property parametrizations in the IFS scheme have similar sensitivities to surface albedo.The increase in global radiation as a function of surface albedo is less in the case of hlradia.When DISORT is used the increase in global radiation as a function of the surface albedo is slightly larger than for each of the HARMONIE schemes.Experiment 7 was run for a fixed integrated cloud water load of 1.0 kg/m 2 and was otherwise similar to experiment 6.As in experiment 6, it can be seen that the global radiation sensitivity in DISORT as a function of the surface albedo is larger than for the HARMONIE cases (Fig. 20).

Ice cloud experiments
Four ice cloud experiments were run to test the SW radiation sensitivity to cloud ice load and cloud ice particle equivalent radius (r e, ice ) in the case of a typical land surface albedo of 0.18 and a typical snow surface albedo of 0.7.Experiments 8 and 9 were run for a fixed r e, ice of 50 µm.This is a typical value of r e, ice for pure ice clouds (Fu, 1996;Nielsen, 2011).

Experiments 8 and 9: Cloud ice load
For the IFS radiation scheme, we tested the Fu and Liou (1993) ice cloud parametrization, which is the default in Harmonie 37h1, and the more recent Fu (1996) parametrization.As can be seen in Fig. 21 both of these parametrizations compare very well to the DISORT calculations in the cloud ice load sensitivity test with low albedo (experiment 8).In the high albedo case (experiment 9) both parametrizations show a significant negative bias for cloud ice loads of 0.2 kg/m 2 and above as shown in Fig. 22.The hlradia scheme performs worse than the IFS schemes.Again, this can to a large extent be explained by the clear sky positive bias of hlradia.With the clear sky bias removed, the hlradia results in experiments 8 and 9 (Figs.21 and 22) would be around 4% lower and    It can be seen that the Fu (1996) parametrization shows better agreement with the DISORT calculations than the Fu and Liou (1993) parametrization.
Both of these show a negative bias of approximately -10% for the smallest r e, ice in the low and high surface albedo cases (Fig. 23 and 24).In Figs.25-27 the net flux differences relative to DISORT are shown for the Fu and Liou (1993) parametrization, the Fu (1996) parametrization and the hlradia scheme (Wyser et al., 1999), respectively.For the latter a strong positive bias is seen in the net fluxes at all atmospheric levels.In the two first cases the largest biases are seen under the clouds with the largest cloud ice loads.This is a similar pattern to that seen in the cloud liquid water load tests (Figs. 11 and 12).
In all parametrizations of cloud optical properties studied here, cloud ice was assumed to consist of hexagonal crystals.In reality, cloud ice particles come in multiple shapes (Baker and Lawson, 2006;Lawson et al., 2006).As shown by Kahnert et al. (2008), these shapes significantly affect the SW forcing of the cloud.Thus, the benchmarking quality of the DISORT run, in this case, is no better than the correctness of the basic assumption made on the cloud ice particle shape.It is only correct for ice clouds consisting of hexagonal crystals.In the libRadtran library it is possible to use alternative cloud ice particle shapes.Thus the DISORT runs could have been made for these, however, as all the cloud ice parametrizations in HARMONIE assume only hexagonal crystals, we have only performed tests for this cloud ice particle shape.

Discussion on the angular distribution of diffuse radiances
In the IFS scheme all scattered irradiance, including that transmitted diffusely through clouds, is modelled as having a zenith angle of 53.0 • , i.e. the irradiance transmitted diffusely through or reflected from clouds is assumed to be on average 60% (i.e cos(53.0• )) of the irradiance normal to the surface (Mascart and Bougeault, 2011).An assumption like this is necessary in the delta-Eddington scheme or any scheme that does not resolve the angular distribution of the radiance (Thomas and Stamnes, 2002).In the 30-stream DISORT simulations the angular radiance distribution below clouds is resolved.Even for cases with conservative scattering the use of a fixed average solar zenith angle has been shown to be imprecise (Hopf, 1936;King, 1960;King et al., 1965).By comparing results of DISORT with results from the delta-Eddington radiative transfer scheme where both are given the  In Fig. 28 a few of these test results are shown.These are for a single scattering albedo of 0.99, an asymmetry factor of 0.8 and a solar zenith angle of 53.0 • .A non-reflecting surface is assumed.It can be seen that the IFS delta-Eddington scheme performs well for optical thicknesses of up to 0.5.An increasing positive bias is seen for higher optical thicknesses, peaking for τ = 5 with +4% bias.For very high optical thicknesses, τ > 50, large relative negative biases are seen; however, these are insignificant as can be seen in the upper plot of Fig. 28 since the absolute transmittance, and thereby the absolute error, is very low.In the supplement, it can be seen that this pattern of errors is similar also when the asymmetry factor and the solar zenith angle are varied.
The results of the delta-Eddington tests can now be compared with the previous results.First let us consider the cloud ice test results (experiment 8) shown in Figs.21 and  26.In this experiment both DISORT and the IFS radiation scheme used the cloud optical property parametrization of Fu (1996).Thus the differences must be partly due to the delta-Eddington approximation and partly due to the approximation of using a limited number of spectral bands.In the IFS scheme only 5 bands are used in practice as mentioned in Sect.2.2.The maximum positive bias seen at a cloud ice load of 0.05 kg/m 2 is similar to the one seen in Fig. 28.The increasingly negative biases, seen at increasing cloud ice loads, however, appear to be much more significant than those seen when testing the delta-Eddington scheme directly.These are thus likely to arise from the error made by using the 5 spectral bands.A similar conclusion can be drawn from studying the results when the IFS-Nielsen cloud liquid optical property parametrization was used (Figs. 9 and 12).
The surface reflectance is also important for the angular distribution of the diffuse irradiances.In each of the experiments the surfaces were assumed to have angular distributions of scattering that follow Lambert's cosine law (Lambert, 1760).In reality, such surfaces do not exist (Lommel, 1889;Bhandari et al., 2011), but they are used in NWP and climate models due to the lack of computational resources needed to resolve the angular radiance distribution.In the DISORT version that is currently available in libRadtran the assumption of Lambertian surface reflectance is also made.Thus, we cannot study the effect of using more realistic bidirectional reflectance distribution functions.
We suggest that implementing a variable average solar zenith angle for scattered irradiance could improve the radiative transfer calculations with the delta-Eddington scheme.This suggestion has previously been made by Thomas and Stamnes (2002).

Conclusions and outlook
Overall we have demonstrated an effective method for testing the radiative transfer computations performed in an NWP model.By defining simplified atmospheric and surface conditions in a single-column model, we have full control of the input and can make clean comparisons of the different parametrizations.Such baseline testing is a necessary first step towards studying the sensitivity of NWP model results to the radiation parametrizations, for which integration in time in a realistic evolving 3D atmospheric environment is required.We have found strengths and weaknesses of the IFS and hlradia parametrizations in HARMONIE in the controlled experiments and suggest the improvements outlined below.
Regarding the IFS SW radiation scheme with HAR-MONIE default settings in comparison to highly detailed radiative transfer calculations we found that: -The clear sky computations show that IFS understimates the atmospheric transmittance for the lowest water vapour loads.
-The Fu (1996) cloud ice optical property parametrization compares better than the Fu and Liou (1993) parametrization.However, both compare fairly well, given the assumption that cloud ice particles are hexagonal columns.
-A new optical property parametrization for liquid clouds has been developed.We have shown that this is better than the parametrizations currently available in HARMONIE.
-Assuming climatological ozone profiles induces a SW error of a few percent at most.
Regarding the much simpler hlradia SW scheme which has only one SW spectral band we have found that: -In the clear sky test cases a bias of +3% to +5% is found in most of the experiments.
-The cloud liquid transmittance formula currently used in hlradia performs impressively well, especially considering the simplicity of the hlradia parametrization.
-The cloud ice transmittance calculated by hlradia is within 7% of that calculated by DISORT.
-Assuming a constant ozone SW absorptance induces an error of a few percent at most.Based on the above we propose the following future work: -The current choice of 6 spectral bands in HARMONIE/IFS should be re-assessed, as the first spectral band is irrelevant for NWP modelling.
-The effect of changing SW cloud inhomogeneity factor from 0.7 (0.8) to 1.0 in all schemes should be tested in the framework of 3D HARMONIE experiments.
-The effects on the general 3D NWP results of using the Nielsen cloud liquid optical property parametrization within the IFS scheme should be tested.
-In order to improve the delta-Eddington radiative transfer calculations, the possibility of using a variable average solar zenith angle for diffuse irradiances should be investigated.
-The hlradia gaseous transmission coefficients should be tuned to the DISORT clear sky results presented here and for the other AFGL atmospheric profiles.
-The impact of aerosols needs to be investigated further.

Fig. 1 .
Fig. 1.Experiment 1: Global radiation as a function of integrated water vapour.The results from DISORT with the "mid-latitude summer" (MS) and "mid-latitude winter" (MW) atmospheric profiles are shown with the black and cyan lines, respectively.The corresponding results for hlradia are shown with red and green lines, and the results for IFS are shown with blue and magenta lines.The vertical dashed line marks the reference integrated water vapour used in the other experiments. .

Fig. 2 .
Fig. 2. Experiment 1: Relative differences (%) in net fluxes between the IFS radiation scheme and DISORT shown as a function of integrated water vapour and height.

Fig. 4 .
Fig. 4. As for Fig. 1 but for the solar zenith angle (experiment 2).The vertical dashed line marks the reference solar zenith angle used in the other experiments.The subplot shows the corresponding relative differences defined as (X-DISORT)/DISORT•100%.

Fig. 8 .
Fig. 8. Global radiation as a function of cloud water load (experiment 4).The results from DISORT (green dashed curve), IFS Fouquart with SW inhomogeneity factors of 0.7 (magenta dotted curve) and 1.0 (red curve), and hlradia with SW inhomogeneity factors of 0.8 (blue dashed curve) and 1.0 (black curve) are shown.

Fig. 9 .
Fig. 9.As for Fig. 8 but here the SW inhomogeneity factor is 1.0 in all cases.The results for DISORT (red curve with +), IFS Fouquart (green curve with x), IFS Slingo (blue curve with *), IFS Nielsen (magenta curve with boxes) and hlradia (cyan curve with filled boxes) are shown.The vertical dashed line at 0.1 kg/m 2 marks the cloud water load used in experiments 5 and 6 while the other line at 1.0 marks that used in experiment 7.

Fig. 19 .
Fig. 19.As for Fig. 9 but for surface albedo for the case of a cloud of medium thickness (experiment 6).The vertical dashed line at 0.18 marks the reference albedo used in experiments 4 and 5.

Fig. 20 .
Fig. 20.As for Fig. 19 but for surface albedo for the case of a thick cloud (experiment 7).

Fig. 21 .
Fig. 21.As for Fig. 1 but for varying cloud ice load and fixed surface albedo = 0.18 (experiment 8).The vertical dashed line at 0.1 kg/m 2 marks the cloud ice load used in experiments 10 and 11.

Fig. 23 .
Fig. 23.As for Fig. 1 but for varying re, ice and fixed surface albedo = 0.18 (experiment 10).The vertical dashed line at 50 µm marks the re, ice value used in experiments 8 and 9.
4.3.2Experiments 10 and 11: Cloud ice particle equivalent radius In Figs.23 and 24 the results from the two tests of SW sensitivity to r e, ice are shown.

Fig. 25 .
Fig. 25.As for Fig. 2 but for the IFS radiation scheme using the Fu and Liou (1993) ice cloud optical property parametrization vs DISORT (experiment 8)

Fig. 28 .
Fig. 28.The results in this figure are calculated for an asymmetry factor g = 0.8 and a cosine solar zenith angle µ0 = 0.60.The upper plot shows comparison of DISORT and IFS transmittances for a single scattering albedo of 0.99 as a function of optical thickness τ .The lower plot shows the corresponding relative differences.

Table 1 .
The benchmark radiative transfer experiments.

Table 2 .
Coefficients used for the Nielsen liquid cloud optical parametrization.The columns show the wavelength bands (i) 2-6.

Table 6 .
As for Table5but for experiment 5 absolute differences in W/m 2 .