source: LMDZ6/trunk/libf/phylmd/ecrad/ifs/radiation_scheme.F90 @ 5423

Last change on this file since 5423 was 4773, checked in by idelkadi, 13 months ago
  • Update of Ecrad in LMDZ The same organization of the Ecrad offline version is retained in order to facilitate the updating of Ecrad in LMDZ and the comparison between online and offline results. version 1.6.1 of Ecrad (​https://github.com/lguez/ecrad.git)
  • Implementation of the double call of Ecrad in LMDZ


File size: 26.6 KB
Line 
1SUBROUTINE RADIATION_SCHEME &
2     & (YRADIATION,KIDIA, KFDIA, KLON, KLEV, KAEROSOL, &
3     &  PSOLAR_IRRADIANCE, &
4     &  PMU0, PTEMPERATURE_SKIN, PALBEDO_DIF, PALBEDO_DIR, &
5     &  PSPECTRALEMISS, &
6     &  PCCN_LAND, PCCN_SEA, &
7     &  PGELAM, PGEMU, PLAND_SEA_MASK, &
8     &  PPRESSURE, PTEMPERATURE, &
9     &  PPRESSURE_H, PTEMPERATURE_H, &
10     &  PQ, PCO2, PCH4, PN2O, PNO2, PCFC11, PCFC12, PHCFC22, PCCL4, PO3, &
11     &  PCLOUD_FRAC, PQ_LIQUID, PQ_ICE, PQ_RAIN, PQ_SNOW, &
12     &  PAEROSOL_OLD, PAEROSOL, &
13     &  PFLUX_SW, PFLUX_LW, PFLUX_SW_CLEAR, PFLUX_LW_CLEAR, &
14     &  PFLUX_SW_DN, PFLUX_LW_DN, PFLUX_SW_DN_CLEAR, PFLUX_LW_DN_CLEAR, &
15     &  PFLUX_DIR, PFLUX_DIR_CLEAR, PFLUX_DIR_INTO_SUN, &
16     &  PFLUX_UV, PFLUX_PAR, PFLUX_PAR_CLEAR, &
17     &  PFLUX_SW_DN_TOA, PEMIS_OUT, PLWDERIVATIVE, &
18     &  PSWDIFFUSEBAND, PSWDIRECTBAND)
19
20! RADIATION_SCHEME - Interface to modular radiation scheme
21!
22! (C) Copyright 2015- ECMWF.
23!
24! This software is licensed under the terms of the Apache Licence Version 2.0
25! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
26!
27! In applying this licence, ECMWF does not waive the privileges and immunities
28! granted to it by virtue of its status as an intergovernmental organisation
29! nor does it submit to any jurisdiction.
30!
31! PURPOSE
32! -------
33!   The modular radiation scheme is contained in a separate
34!   library. This routine puts the the IFS arrays into appropriate
35!   objects, computing the additional data that is required, and sends
36!   it to the radiation scheme.  It returns net fluxes and surface
37!   flux components needed by the rest of the model.
38!
39!   Lower case is used for variables and types taken from the
40!   radiation library
41!
42! INTERFACE
43! ---------
44!    RADIATION_SCHEME is called from RADLSWR. The
45!    SETUP_RADIATION_SCHEME routine (in the RADIATION_SETUP module)
46!    populates the YRADIATION object, and should have been run first.
47!
48! AUTHOR
49! ------
50!   Robin Hogan, ECMWF
51!   Original: 2015-09-16
52!
53! MODIFICATIONS
54! -------------
55!   2017-03-03  R. Hogan  Read configuration data from YRADIATION object
56!   2017-05-11  R. Hogan  Pass KIDIA,KFDIA to get_layer_mass
57!   2018-01-11  R. Hogan  Capability to scale solar spectrum in each band
58!   2017-11-11  M. Ahlgrimm add variable FSD for cloud heterogeneity
59!   2017-11-29  R. Hogan  Check fluxes in physical bounds
60!   2019-01-22  R. Hogan  Use fluxes in albedo bands from ecRad
61!   2019-01-23  R. Hogan  Spectral longwave emissivity in NLWEMISS bands
62!   2019-02-04  R. Hogan  Pass out surface longwave downwelling in each emissivity interval
63!   2019-02-07  R. Hogan  SPARTACUS cloud size from PARAM_CLOUD_EFFECTIVE_SEPARATION_ETA
64!
65!-----------------------------------------------------------------------
66
67! Modules from ifs or ifsaux libraries
68USE PARKIND1       , ONLY : JPIM, JPRB, JPRD
69USE YOMHOOK        , ONLY : LHOOK, DR_HOOK, JPHOOK
70USE YOMCST         , ONLY : RPI, RSIGMA ! Stefan-Boltzmann constant
71USE YOMLUN         , ONLY : NULERR
72USE RADIATION_SETUP, ONLY : ITYPE_TROP_BG_AER, ITYPE_STRAT_BG_AER, TRADIATION
73
74! Modules from ecRad radiation library
75USE RADIATION_CONFIG,         ONLY : ISOLVERSPARTACUS
76USE RADIATION_SINGLE_LEVEL,   ONLY : SINGLE_LEVEL_TYPE
77USE RADIATION_THERMODYNAMICS, ONLY : THERMODYNAMICS_TYPE
78USE RADIATION_GAS,            ONLY : GAS_TYPE,&
79     &                               IMASSMIXINGRATIO, IVOLUMEMIXINGRATIO,&
80     &                               IH2O, ICO2, ICH4, IN2O, ICFC11, ICFC12, IHCFC22, ICCL4, IO3, IO2
81USE RADIATION_CLOUD,          ONLY : CLOUD_TYPE
82USE RADIATION_AEROSOL,        ONLY : AEROSOL_TYPE
83USE RADIATION_FLUX,           ONLY : FLUX_TYPE
84USE RADIATION_INTERFACE,      ONLY : RADIATION, SET_GAS_UNITS
85USE RADIATION_SAVE,           ONLY : SAVE_INPUTS, SAVE_FLUXES
86
87IMPLICIT NONE
88
89! INPUT ARGUMENTS
90
91TYPE(TRADIATION), INTENT(IN)    :: YRADIATION
92
93! *** Array dimensions and ranges
94INTEGER(KIND=JPIM),INTENT(IN)   :: KIDIA    ! Start column to process
95INTEGER(KIND=JPIM),INTENT(IN)   :: KFDIA    ! End column to process
96INTEGER(KIND=JPIM),INTENT(IN)   :: KLON     ! Number of columns
97INTEGER(KIND=JPIM),INTENT(IN)   :: KLEV     ! Number of levels
98INTEGER(KIND=JPIM),INTENT(IN)   :: KAEROSOL ! Number of aerosol types
99
100! *** Single-level fields
101REAL(KIND=JPRB),   INTENT(IN) :: PSOLAR_IRRADIANCE ! (W m-2)
102REAL(KIND=JPRB),   INTENT(IN) :: PMU0(KLON) ! Cosine of solar zenith ang
103REAL(KIND=JPRB),   INTENT(IN) :: PTEMPERATURE_SKIN(KLON) ! (K)
104! Diffuse and direct components of surface shortwave albedo
105REAL(KIND=JPRB),   INTENT(IN) :: PALBEDO_DIF(KLON,YRADIATION%YRERAD%NSW)
106REAL(KIND=JPRB),   INTENT(IN) :: PALBEDO_DIR(KLON,YRADIATION%YRERAD%NSW)
107! Longwave spectral emissivity
108REAL(KIND=JPRB),   INTENT(IN) :: PSPECTRALEMISS(KLON,YRADIATION%YRERAD%NLWEMISS)
109! Longitude (radians), sine of latitude
110REAL(KIND=JPRB),   INTENT(IN) :: PGELAM(KLON)
111REAL(KIND=JPRB),   INTENT(IN) :: PGEMU(KLON)
112! Land-sea mask
113REAL(KIND=JPRB),   INTENT(IN) :: PLAND_SEA_MASK(KLON)
114
115! *** Variables on full levels
116REAL(KIND=JPRB),   INTENT(IN) :: PPRESSURE(KLON,KLEV)    ! (Pa)
117REAL(KIND=JPRB),   INTENT(IN) :: PTEMPERATURE(KLON,KLEV) ! (K)
118! *** Variables on half levels
119REAL(KIND=JPRB),   INTENT(IN) :: PPRESSURE_H(KLON,KLEV+1)    ! (Pa)
120REAL(KIND=JPRB),   INTENT(IN) :: PTEMPERATURE_H(KLON,KLEV+1) ! (K)
121
122! *** Gas mass mixing ratios on full levels
123REAL(KIND=JPRB),   INTENT(IN) :: PQ(KLON,KLEV)
124REAL(KIND=JPRB),   INTENT(IN) :: PCO2(KLON,KLEV)
125REAL(KIND=JPRB),   INTENT(IN) :: PCH4(KLON,KLEV)
126REAL(KIND=JPRB),   INTENT(IN) :: PN2O(KLON,KLEV)
127REAL(KIND=JPRB),   INTENT(IN) :: PNO2(KLON,KLEV)
128REAL(KIND=JPRB),   INTENT(IN) :: PCFC11(KLON,KLEV)
129REAL(KIND=JPRB),   INTENT(IN) :: PCFC12(KLON,KLEV)
130REAL(KIND=JPRB),   INTENT(IN) :: PHCFC22(KLON,KLEV)
131REAL(KIND=JPRB),   INTENT(IN) :: PCCL4(KLON,KLEV)
132REAL(KIND=JPRB),   INTENT(IN) :: PO3(KLON,KLEV)
133
134! *** Cloud fraction and hydrometeor mass mixing ratios
135REAL(KIND=JPRB),   INTENT(IN) :: PCLOUD_FRAC(KLON,KLEV)
136REAL(KIND=JPRB),   INTENT(IN) :: PQ_LIQUID(KLON,KLEV)
137REAL(KIND=JPRB),   INTENT(IN) :: PQ_ICE(KLON,KLEV)
138REAL(KIND=JPRB),   INTENT(IN) :: PQ_RAIN(KLON,KLEV)
139REAL(KIND=JPRB),   INTENT(IN) :: PQ_SNOW(KLON,KLEV)
140
141! *** Aerosol mass mixing ratios
142REAL(KIND=JPRB),   INTENT(IN) :: PAEROSOL_OLD(KLON,6,KLEV)
143REAL(KIND=JPRB),   INTENT(IN) :: PAEROSOL(KLON,KLEV,KAEROSOL)
144
145REAL(KIND=JPRB),   INTENT(IN) :: PCCN_LAND(KLON)
146REAL(KIND=JPRB),   INTENT(IN) :: PCCN_SEA(KLON)
147
148! OUTPUT ARGUMENTS
149
150! *** Net fluxes on half-levels (W m-2)
151REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW(KLON,KLEV+1)
152REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW(KLON,KLEV+1)
153REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_CLEAR(KLON,KLEV+1)
154REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_CLEAR(KLON,KLEV+1)
155
156! *** Surface flux components (W m-2)
157REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_DN(KLON)
158REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_DN(KLON)
159REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_DN_CLEAR(KLON)
160REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_DN_CLEAR(KLON)
161! Direct component of surface flux into horizontal plane
162REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_DIR(KLON)
163REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_DIR_CLEAR(KLON)
164! As PFLUX_DIR but into a plane perpendicular to the sun
165REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_DIR_INTO_SUN(KLON)
166
167! *** Ultraviolet and photosynthetically active radiation (W m-2)
168REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_UV(KLON)
169REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_PAR(KLON)
170REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_PAR_CLEAR(KLON)
171
172! *** Other single-level diagnostics
173! Top-of-atmosphere incident solar flux (W m-2)
174REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_DN_TOA(KLON)
175! Diagnosed longwave surface emissivity across the whole spectrum
176REAL(KIND=JPRB),  INTENT(OUT) :: PEMIS_OUT(KLON)
177
178! Partial derivative of total-sky longwave upward flux at each level
179! with respect to upward flux at surface, used to correct heating
180! rates at gridpoints/timesteps between calls to the full radiation
181! scheme.  Note that this version uses the convention of level index
182! increasing downwards, unlike the local variable ZLwDerivative that
183! is returned from the LW radiation scheme.
184REAL(KIND=JPRB),  INTENT(OUT) :: PLWDERIVATIVE(KLON,KLEV+1)
185
186! Surface diffuse and direct downwelling shortwave flux in each
187! shortwave albedo band, used in RADINTG to update the surface fluxes
188! accounting for high-resolution albedo information
189REAL(KIND=JPRB),  INTENT(OUT) :: PSWDIFFUSEBAND(KLON,YRADIATION%YRERAD%NSW)
190REAL(KIND=JPRB),  INTENT(OUT) :: PSWDIRECTBAND (KLON,YRADIATION%YRERAD%NSW)
191
192! LOCAL VARIABLES
193TYPE(SINGLE_LEVEL_TYPE)   :: SINGLE_LEVEL
194TYPE(THERMODYNAMICS_TYPE) :: THERMODYNAMICS
195TYPE(GAS_TYPE)            :: GAS
196TYPE(CLOUD_TYPE)          :: YLCLOUD
197TYPE(AEROSOL_TYPE)        :: AEROSOL
198TYPE(FLUX_TYPE)           :: FLUX
199
200! Cloud effective radii in microns
201REAL(KIND=JPRB)           :: ZRE_LIQUID_UM(KLON,KLEV)
202REAL(KIND=JPRB)           :: ZRE_ICE_UM(KLON,KLEV)
203
204! Cloud overlap decorrelation length for cloud boundaries in km
205REAL(KIND=JPRB)           :: ZDECORR_LEN_KM(KLON)
206
207! Ratio of cloud overlap decorrelation length for cloud water
208! inhomogeneities to that for cloud boundaries (typically 0.5)
209REAL(KIND=JPRB)           :: ZDECORR_LEN_RATIO
210
211! The surface net longwave flux if the surface was a black body, used
212! to compute the effective broadband surface emissivity
213REAL(KIND=JPRB)           :: ZBLACK_BODY_NET_LW(KIDIA:KFDIA)
214
215! Layer mass in kg m-2
216REAL(KIND=JPRB)           :: ZLAYER_MASS(KIDIA:KFDIA,KLEV)
217
218! Time integers
219! INTEGER(KIND=JPIM) :: ITIM, IDAY
220
221! Loop indices
222INTEGER(KIND=JPIM) :: JLON, JLEV, JBAND, JAER
223
224! Have any fluxes been returned that are out of a physically
225! reasonable range? This integer stores the number of blocks of fluxes
226! that have contained a bad value so far, for this task.  NetCDF files
227! will be written up to the value of NAERAD:NDUMPBADINPUTS.
228INTEGER(KIND=JPIM), SAVE :: N_BAD_FLUXES = 0
229
230! For debugging it can be useful to save input profiles and output
231! fluxes without the condition that the fluxes are out of a reasonable
232! range. NetCDF files will be written up to the value of
233! NAERAD:NDUMPINPUTS.
234INTEGER(KIND=JPIM), SAVE :: N_OUTPUT_FLUXES = 0
235
236! NetCDF file name in case of bad fluxes
237CHARACTER(LEN=512) :: CL_FILE_NAME
238
239REAL(KIND=JPHOOK) :: ZHOOK_HANDLE
240
241! Dummy from YOMCT3
242! INTEGER(KIND=JPIM) :: NSTEP = 0
243
244! Dummy from MPL_MYRANK_MOD
245INTEGER(KIND=JPIM) :: MPL_MYRANK
246MPL_MYRANK() = 1
247
248! Import time functions for iseed calculation
249#include "fcttim.func.h"
250
251#include "liquid_effective_radius.intfb.h"
252#include "ice_effective_radius.intfb.h"
253#include "cloud_overlap_decorr_len.intfb.h"
254!#include "satur.intfb.h"
255!#include "abor1.intfb.h"
256
257IF (LHOOK) CALL DR_HOOK('RADIATION_SCHEME',0,ZHOOK_HANDLE)
258
259ASSOCIATE(YRERAD    =>YRADIATION%YRERAD, &
260     &    RAD_CONFIG=>YRADIATION%RAD_CONFIG, &
261     &    NWEIGHT_UV=>YRADIATION%NWEIGHT_UV, &
262     &    IBAND_UV  =>YRADIATION%IBAND_UV(:), &
263     &    WEIGHT_UV =>YRADIATION%WEIGHT_UV(:), &
264     &    NWEIGHT_PAR=>YRADIATION%NWEIGHT_PAR, &
265     &    IBAND_PAR =>YRADIATION%IBAND_PAR(:), &
266     &    WEIGHT_PAR=>YRADIATION%WEIGHT_PAR(:), &
267     &    TROP_BG_AER_MASS_EXT=>YRADIATION%TROP_BG_AER_MASS_EXT, &
268     &    STRAT_BG_AER_MASS_EXT=>YRADIATION%STRAT_BG_AER_MASS_EXT)
269! Allocate memory in radiation objects
270CALL SINGLE_LEVEL%ALLOCATE(KLON, YRERAD%NSW, YRERAD%NLWEMISS, &
271     &                     USE_SW_ALBEDO_DIRECT=.TRUE.)
272CALL THERMODYNAMICS%ALLOCATE(KLON, KLEV, USE_H2O_SAT=.TRUE.)
273CALL GAS%ALLOCATE(KLON, KLEV)
274CALL YLCLOUD%ALLOCATE(KLON, KLEV)
275IF (YRERAD%NAERMACC == 1) THEN
276  CALL AEROSOL%ALLOCATE(KLON, 1, KLEV, KAEROSOL) ! MACC aerosols
277ELSE
278  CALL AEROSOL%ALLOCATE(KLON, 1, KLEV, 6) ! Tegen climatology
279ENDIF
280CALL FLUX%ALLOCATE(RAD_CONFIG, 1, KLON, KLEV)
281
282! Set thermodynamic profiles: simply copy over the half-level
283! pressure and temperature
284THERMODYNAMICS%PRESSURE_HL   (KIDIA:KFDIA,:) = PPRESSURE_H   (KIDIA:KFDIA,:)
285THERMODYNAMICS%TEMPERATURE_HL(KIDIA:KFDIA,:) = PTEMPERATURE_H(KIDIA:KFDIA,:)
286
287! IFS currently sets the half-level temperature at the surface to be
288! equal to the skin temperature. The radiation scheme takes as input
289! only the half-level temperatures and assumes the Planck function to
290! vary linearly in optical depth between half levels. In the lowest
291! atmospheric layer, where the atmospheric temperature can be much
292! cooler than the skin temperature, this can lead to significant
293! differences between the effective temperature of this lowest layer
294! and the true value in the model.
295!
296! We may approximate the temperature profile in the lowest model level
297! as piecewise linear between the top of the layer T[k-1/2], the
298! centre of the layer T[k] and the base of the layer Tskin.  The mean
299! temperature of the layer is then 0.25*T[k-1/2] + 0.5*T[k] +
300! 0.25*Tskin, which can be achieved by setting the atmospheric
301! temperature at the half-level corresponding to the surface as
302! follows:
303THERMODYNAMICS%TEMPERATURE_HL(KIDIA:KFDIA,KLEV+1)&
304     &  = PTEMPERATURE(KIDIA:KFDIA,KLEV)&
305     &  + 0.5_JPRB * (PTEMPERATURE_H(KIDIA:KFDIA,KLEV+1)&
306     &               -PTEMPERATURE_H(KIDIA:KFDIA,KLEV))
307
308! Alternatively we respect the model's atmospheric temperature in the
309! lowest model level by setting the temperature at the lowest
310! half-level such that the mean temperature of the layer is correct:
311!thermodynamics%temperature_hl(KIDIA:KFDIA,KLEV+1) &
312!     &  = 2.0_JPRB * PTEMPERATURE(KIDIA:KFDIA,KLEV) &
313!     &             - PTEMPERATURE_H(KIDIA:KFDIA,KLEV)
314
315! Compute saturation specific humidity, used to hydrate aerosols. The
316! "2" for the last argument indicates that the routine is not being
317! called from within the convection scheme.
318!CALL SATUR(KIDIA, KFDIA, KLON, 1, KLEV, .false., &
319!     &  PPRESSURE, PTEMPERATURE, THERMODYNAMICS%H2O_SAT_LIQ, 2)
320! Alternative approximate version using temperature and pressure from
321! the thermodynamics structure
322CALL thermodynamics%calc_saturation_wrt_liquid(KIDIA, KFDIA)
323
324! Set single-level fileds
325SINGLE_LEVEL%SOLAR_IRRADIANCE              = PSOLAR_IRRADIANCE
326SINGLE_LEVEL%COS_SZA(KIDIA:KFDIA)          = PMU0(KIDIA:KFDIA)
327SINGLE_LEVEL%SKIN_TEMPERATURE(KIDIA:KFDIA) = PTEMPERATURE_SKIN(KIDIA:KFDIA)
328SINGLE_LEVEL%SW_ALBEDO(KIDIA:KFDIA,:)      = PALBEDO_DIF(KIDIA:KFDIA,:)
329SINGLE_LEVEL%SW_ALBEDO_DIRECT(KIDIA:KFDIA,:)=PALBEDO_DIR(KIDIA:KFDIA,:)
330! Spectral longwave emissivity
331SINGLE_LEVEL%LW_EMISSIVITY(KIDIA:KFDIA,:)  = PSPECTRALEMISS(KIDIA:KFDIA,:)
332
333! Create the relevant seed from date and time get the starting day
334! and number of minutes since start
335! IDAY = NDD(NINDAT)
336! ITIM = NINT(NSTEP * YDMODEL%YRML_GCONF%YRRIP%TSTEP / 60.0_JPRB)
337! DO JLON = KIDIA, KFDIA
338!   ! This method gives a unique value for roughly every 1-km square
339!   ! on the globe and every minute.  ASIN(PGEMU)*60 gives rough
340!   ! latitude in degrees, which we multiply by 100 to give a unique
341!   ! value for roughly every km. PGELAM*60*100 gives a unique number
342!   ! for roughly every km of longitude around the equator, which we
343!   ! multiply by 180*100 so there is no overlap with the latitude
344!   ! values.  The result can be contained in a 32-byte integer (but
345!   ! since random numbers are generated with the help of integer
346!   ! overflow, it should not matter if the number did overflow).
347!   SINGLE_LEVEL%ISEED(JLON) = ITIM + IDAY &
348!        &  +  NINT(PGELAM(JLON)*108000000.0_JPRD &
349!        &          + ASIN(PGEMU(JLON))*6000.0_JPRD)
350! ENDDO
351
352! Simple initialization of the seeds for the Monte Carlo scheme
353call single_level%init_seed_simple(kidia, kfdia)
354
355! Set the solar spectrum scaling, if required
356IF (YRERAD%NSOLARSPECTRUM == 1) THEN
357  ALLOCATE(SINGLE_LEVEL%SPECTRAL_SOLAR_SCALING(RAD_CONFIG%N_BANDS_SW))
358  ! Ratio of SORCE (Coddington et al. 2016) and Kurucz solar spectra
359  SINGLE_LEVEL%SPECTRAL_SOLAR_SCALING &
360       &  = (/  1.0, 1.0, 1.0, 1.0478, 1.0404, 1.0317, 1.0231, &
361       &        1.0054, 0.98413, 0.99863, 0.99907, 0.90589, 0.92213, 1.0 /)
362ENDIF
363
364! Set cloud fields
365YLCLOUD%Q_LIQ(KIDIA:KFDIA,:)    = PQ_LIQUID(KIDIA:KFDIA,:)
366YLCLOUD%Q_ICE(KIDIA:KFDIA,:)    = PQ_ICE(KIDIA:KFDIA,:) + PQ_SNOW(KIDIA:KFDIA,:)
367YLCLOUD%FRACTION(KIDIA:KFDIA,:) = PCLOUD_FRAC(KIDIA:KFDIA,:)
368
369! Compute effective radii and convert to metres
370CALL LIQUID_EFFECTIVE_RADIUS(YRERAD, &
371     &  KIDIA, KFDIA, KLON, KLEV, &
372     &  PPRESSURE, PTEMPERATURE, PCLOUD_FRAC, PQ_LIQUID, PQ_RAIN, &
373     &  PLAND_SEA_MASK, PCCN_LAND, PCCN_SEA, &
374     &  ZRE_LIQUID_UM) !, PPERT=PPERT)
375YLCLOUD%RE_LIQ(KIDIA:KFDIA,:) = ZRE_LIQUID_UM(KIDIA:KFDIA,:) * 1.0E-6_JPRB
376
377CALL ICE_EFFECTIVE_RADIUS(YRERAD, KIDIA, KFDIA, KLON, KLEV, &
378     &  PPRESSURE, PTEMPERATURE, PCLOUD_FRAC, PQ_ICE, PQ_SNOW, PGEMU, &
379     &  ZRE_ICE_UM) !, PPERT=PPERT)
380YLCLOUD%RE_ICE(KIDIA:KFDIA,:) = ZRE_ICE_UM(KIDIA:KFDIA,:) * 1.0E-6_JPRB
381
382! Get the cloud overlap decorrelation length (for cloud boundaries),
383! in km, according to the parameterization specified by NDECOLAT,
384! and insert into the "cloud" object. Also get the ratio of
385! decorrelation lengths for cloud water content inhomogeneities and
386! cloud boundaries, and set it in the "rad_config" object.
387CALL CLOUD_OVERLAP_DECORR_LEN(KIDIA,KFDIA,KLON, &
388     &  PGEMU,YRERAD%NDECOLAT, &
389     &  PDECORR_LEN_EDGES_KM=ZDECORR_LEN_KM, PDECORR_LEN_RATIO=ZDECORR_LEN_RATIO)
390
391! Compute cloud overlap parameter from decorrelation length
392!RAD_CONFIG%CLOUD_INHOM_DECORR_SCALING = ZDECORR_LEN_RATIO
393DO JLON = KIDIA,KFDIA
394  CALL YLCLOUD%SET_OVERLAP_PARAM(THERMODYNAMICS,&
395      &                       ZDECORR_LEN_KM(JLON)*1000.0_JPRB,&
396      &                       ISTARTCOL=JLON, IENDCOL=JLON)
397ENDDO
398! Or we can call the routine on all columns at once
399!CALL YLCLOUD%SET_OVERLAP_PARAM(THERMODYNAMICS,&
400!     &                       ZDECORR_LEN_KM(KIDIA:KFDIA)*1000.0_JPRB,&
401!     &                       ISTARTCOL=KIDIA, IENDCOL=KFDIA)
402
403! Cloud water content fractional standard deviation is configurable
404! from namelist NAERAD but must be globally constant. Before it was
405! hard coded at 1.0.
406CALL YLCLOUD%CREATE_FRACTIONAL_STD(KLON, KLEV, YRERAD%RCLOUD_FRAC_STD)
407
408
409IF (         RAD_CONFIG%I_SOLVER_LW == ISOLVERSPARTACUS &
410     &  .OR. RAD_CONFIG%I_SOLVER_SW == ISOLVERSPARTACUS) THEN
411  ! We are using the SPARTACUS solver so need to specify cloud scale,
412  ! and use Mark Fielding's parameterization based on ARM data
413  CALL YLCLOUD%PARAM_CLOUD_EFFECTIVE_SEPARATION_ETA(KLON, KLEV, &
414       &  PPRESSURE_H, YRERAD%RCLOUD_SEPARATION_SCALE_SURF, &
415       &  YRERAD%RCLOUD_SEPARATION_SCALE_TOA, 3.5_JPRB, 0.75_JPRB, &
416       &  KIDIA, KFDIA)
417ENDIF
418
419! Compute the dry mass of each layer neglecting humidity effects, in
420! kg m-2, needed to scale some of the aerosol inputs
421CALL THERMODYNAMICS%GET_LAYER_MASS(KIDIA,KFDIA,ZLAYER_MASS)
422
423! Copy over aerosol mass mixing ratio
424IF (YRERAD%NAERMACC == 1) THEN
425
426
427  ! MACC aerosol from climatology or prognostic aerosol variables -
428  ! this is already in mass mixing ratio units with the required array
429  ! orientation so we can copy it over directly
430  ! AB need to cap the minimum mass mixing ratio/AOD to avoid instability
431  ! in case of negative values in input
432  DO JAER = 1,KAEROSOL
433    DO JLEV = 1,KLEV
434      DO JLON = KIDIA,KFDIA
435        AEROSOL%MIXING_RATIO(JLON,JLEV,JAER) = MAX(PAEROSOL(JLON,JLEV,JAER),0.0_JPRB)
436      ENDDO
437    ENDDO
438  ENDDO
439
440  IF (YRERAD%NAERMACC == 1) THEN
441    ! Add the tropospheric and stratospheric backgrounds contained in the
442    ! old Tegen arrays - this is very ugly!
443    IF (TROP_BG_AER_MASS_EXT > 0.0_JPRB) THEN
444      AEROSOL%MIXING_RATIO(KIDIA:KFDIA,:,ITYPE_TROP_BG_AER)&
445           &  = AEROSOL%MIXING_RATIO(KIDIA:KFDIA,:,ITYPE_TROP_BG_AER)&
446           &  + PAEROSOL_OLD(KIDIA:KFDIA,1,:)&
447           &  / (ZLAYER_MASS * TROP_BG_AER_MASS_EXT)
448    ENDIF
449    IF (STRAT_BG_AER_MASS_EXT > 0.0_JPRB) THEN
450      AEROSOL%MIXING_RATIO(KIDIA:KFDIA,:,ITYPE_STRAT_BG_AER)&
451           &  = AEROSOL%MIXING_RATIO(KIDIA:KFDIA,:,ITYPE_STRAT_BG_AER)&
452           &  + PAEROSOL_OLD(KIDIA:KFDIA,6,:)&
453           &  / (ZLAYER_MASS * STRAT_BG_AER_MASS_EXT)
454    ENDIF
455  ENDIF
456ELSE
457
458  ! Tegen aerosol climatology - the array PAEROSOL_OLD contains the
459  ! 550-nm optical depth in each layer. The optics data file
460  ! aerosol_ifs_rrtm_tegen.nc does not contain mass extinction
461  ! coefficient, but a scaling factor that the 550-nm optical depth
462  ! should be multiplied by to obtain the optical depth in each
463  ! spectral band.  Therefore, in order for the units to work out, we
464  ! need to divide by the layer mass (in kg m-2) to obtain the 550-nm
465  ! cross-section per unit mass of dry air (so in m2 kg-1).  We also
466  ! need to permute the array.
467  DO JLEV = 1,KLEV
468    DO JAER = 1,6
469      AEROSOL%MIXING_RATIO(KIDIA:KFDIA,JLEV,JAER)&
470         &  = PAEROSOL_OLD(KIDIA:KFDIA,JAER,JLEV)&
471         &  / ZLAYER_MASS(KIDIA:KFDIA,JLEV)
472    ENDDO
473  ENDDO
474
475ENDIF
476
477! Insert gas mixing ratios
478CALL GAS%PUT(IH2O,    IMASSMIXINGRATIO, PQ)
479CALL GAS%PUT(ICO2,    IMASSMIXINGRATIO, PCO2)
480CALL GAS%PUT(ICH4,    IMASSMIXINGRATIO, PCH4)
481CALL GAS%PUT(IN2O,    IMASSMIXINGRATIO, PN2O)
482CALL GAS%PUT(ICFC11,  IMASSMIXINGRATIO, PCFC11)
483CALL GAS%PUT(ICFC12,  IMASSMIXINGRATIO, PCFC12)
484CALL GAS%PUT(IHCFC22, IMASSMIXINGRATIO, PHCFC22)
485CALL GAS%PUT(ICCL4,   IMASSMIXINGRATIO, PCCL4)
486CALL GAS%PUT(IO3,     IMASSMIXINGRATIO, PO3)
487CALL GAS%PUT_WELL_MIXED(IO2, IVOLUMEMIXINGRATIO, 0.20944_JPRB)
488
489! Ensure the units of the gas mixing ratios are what is required by
490! the gas absorption model
491CALL SET_GAS_UNITS(RAD_CONFIG, GAS)
492
493!call save_inputs('inputs_ifs.nc', rad_config, single_level, thermodynamics, &
494!     &           gas, ylcloud, aerosol, &
495!     &           lat=spread(0.0_jprb,1,klon), &
496!     &           lon=spread(0.0_jprb,1,klon), &
497!     &           iverbose=2)
498
499! Call radiation scheme
500CALL RADIATION(KLON, KLEV, KIDIA, KFDIA, RAD_CONFIG,&
501     &  SINGLE_LEVEL, THERMODYNAMICS, GAS, YLCLOUD, AEROSOL, FLUX)
502
503! Check fluxes are within physical bounds
504IF (YRERAD%NDUMPBADINPUTS /= 0 &
505     &  .AND. (N_BAD_FLUXES == 0 .OR. N_BAD_FLUXES < YRERAD%NDUMPBADINPUTS)) THEN
506  IF (FLUX%OUT_OF_PHYSICAL_BOUNDS(KIDIA,KFDIA)) THEN
507!$OMP CRITICAL
508    N_BAD_FLUXES = N_BAD_FLUXES+1
509    WRITE(CL_FILE_NAME, '(A,I0,A,I0,A)') '/home/parr/ifs_dump/bad_inputs_', &
510         &  MPL_MYRANK(), '_', N_BAD_FLUXES, '.nc'
511    WRITE(NULERR,*) '  Writing ', TRIM(CL_FILE_NAME)
512    ! Implicit assumption that KFDIA==KLON
513    CALL SAVE_INPUTS(TRIM(CL_FILE_NAME), RAD_CONFIG, SINGLE_LEVEL, &
514         &  THERMODYNAMICS, GAS, YLCLOUD, AEROSOL, &
515         &  LAT=ASIN(PGEMU)*180.0/RPI, LON=PGELAM*180.0/RPI, IVERBOSE=3)
516    WRITE(CL_FILE_NAME, '(A,I0,A,I0,A)') '/home/parr/ifs_dump/bad_outputs_', &
517         &  MPL_MYRANK(), '_', N_BAD_FLUXES, '.nc'
518    WRITE(NULERR,*) '  Writing ', TRIM(CL_FILE_NAME)
519    CALL SAVE_FLUXES(TRIM(CL_FILE_NAME), RAD_CONFIG, THERMODYNAMICS, FLUX, IVERBOSE=3)
520    IF (YRERAD%NDUMPBADINPUTS < 0) THEN
521      ! Abort on the first set of bad fluxes
522      CALL ABOR1("RADIATION_SCHEME: ABORT DUE TO FLUXES OUT OF PHYSICAL BOUNDS")
523    ENDIF
524!$OMP END CRITICAL
525  ENDIF
526ENDIF
527
528! For debugging, do we store a certain number of inputs and outputs
529! regardless of whether bad fluxes have been detected?
530IF (N_OUTPUT_FLUXES < YRERAD%NDUMPINPUTS) THEN
531!$OMP CRITICAL
532  N_OUTPUT_FLUXES = N_OUTPUT_FLUXES+1
533  WRITE(CL_FILE_NAME, '(A,I0,A,I0,A)') '/home/parr/ifs_dump/inputs_', &
534       &  MPL_MYRANK(), '_', N_OUTPUT_FLUXES, '.nc'
535  WRITE(NULERR,*) '  Writing ', TRIM(CL_FILE_NAME)
536  ! Implicit assumption that KFDIA==KLON
537  CALL SAVE_INPUTS(TRIM(CL_FILE_NAME), RAD_CONFIG, SINGLE_LEVEL, &
538       &  THERMODYNAMICS, GAS, YLCLOUD, AEROSOL, &
539       &  LAT=ASIN(PGEMU)*180.0/RPI, LON=PGELAM*180.0/RPI, IVERBOSE=3)
540  WRITE(CL_FILE_NAME, '(A,I0,A,I0,A)') '/home/parr/ifs_dump/outputs_', &
541       &  MPL_MYRANK(), '_', N_OUTPUT_FLUXES, '.nc'
542  WRITE(NULERR,*) '  Writing ', TRIM(CL_FILE_NAME)
543  CALL SAVE_FLUXES(TRIM(CL_FILE_NAME), RAD_CONFIG, THERMODYNAMICS, FLUX, IVERBOSE=3)
544!$OMP END CRITICAL
545ENDIF
546
547! Compute required output fluxes
548! First the net fluxes
549PFLUX_SW(KIDIA:KFDIA,:) = FLUX%SW_DN(KIDIA:KFDIA,:) - FLUX%SW_UP(KIDIA:KFDIA,:)
550PFLUX_LW(KIDIA:KFDIA,:) = FLUX%LW_DN(KIDIA:KFDIA,:) - FLUX%LW_UP(KIDIA:KFDIA,:)
551PFLUX_SW_CLEAR(KIDIA:KFDIA,:)&
552     &  = FLUX%SW_DN_CLEAR(KIDIA:KFDIA,:) - FLUX%SW_UP_CLEAR(KIDIA:KFDIA,:)
553PFLUX_LW_CLEAR(KIDIA:KFDIA,:)&
554     &  = FLUX%LW_DN_CLEAR(KIDIA:KFDIA,:) - FLUX%LW_UP_CLEAR(KIDIA:KFDIA,:)
555! Now the surface fluxes
556PFLUX_SW_DN      (KIDIA:KFDIA) = FLUX%SW_DN             (KIDIA:KFDIA,KLEV+1)
557PFLUX_LW_DN      (KIDIA:KFDIA) = FLUX%LW_DN             (KIDIA:KFDIA,KLEV+1)
558PFLUX_SW_DN_CLEAR(KIDIA:KFDIA) = FLUX%SW_DN_CLEAR       (KIDIA:KFDIA,KLEV+1)
559PFLUX_LW_DN_CLEAR(KIDIA:KFDIA) = FLUX%LW_DN_CLEAR       (KIDIA:KFDIA,KLEV+1)
560PFLUX_DIR        (KIDIA:KFDIA) = FLUX%SW_DN_DIRECT      (KIDIA:KFDIA,KLEV+1)
561PFLUX_DIR_CLEAR  (KIDIA:KFDIA) = FLUX%SW_DN_DIRECT_CLEAR(KIDIA:KFDIA,KLEV+1)
562PFLUX_DIR_INTO_SUN(KIDIA:KFDIA) = 0.0_JPRB
563WHERE (PMU0(KIDIA:KFDIA) > EPSILON(1.0_JPRB))
564  PFLUX_DIR_INTO_SUN(KIDIA:KFDIA) = PFLUX_DIR(KIDIA:KFDIA) / PMU0(KIDIA:KFDIA)
565ENDWHERE
566! Top-of-atmosphere downwelling flux
567PFLUX_SW_DN_TOA(KIDIA:KFDIA) = FLUX%SW_DN(KIDIA:KFDIA,1)
568
569! Compute UV fluxes as weighted sum of appropriate shortwave bands
570PFLUX_UV       (KIDIA:KFDIA) = 0.0_JPRB
571DO JBAND = 1,NWEIGHT_UV
572!DEC$ IVDEP
573  PFLUX_UV(KIDIA:KFDIA) = PFLUX_UV(KIDIA:KFDIA) + WEIGHT_UV(JBAND)&
574       &  * FLUX%SW_DN_SURF_BAND(IBAND_UV(JBAND),KIDIA:KFDIA)
575ENDDO
576
577! Compute photosynthetically active radiation similarly
578PFLUX_PAR      (KIDIA:KFDIA) = 0.0_JPRB
579PFLUX_PAR_CLEAR(KIDIA:KFDIA) = 0.0_JPRB
580DO JBAND = 1,NWEIGHT_PAR
581!DEC$ IVDEP
582  PFLUX_PAR(KIDIA:KFDIA) = PFLUX_PAR(KIDIA:KFDIA) + WEIGHT_PAR(JBAND)&
583       &  * FLUX%SW_DN_SURF_BAND(IBAND_PAR(JBAND),KIDIA:KFDIA)
584!DEC$ IVDEP
585  PFLUX_PAR_CLEAR(KIDIA:KFDIA) = PFLUX_PAR_CLEAR(KIDIA:KFDIA)&
586       &  + WEIGHT_PAR(JBAND)&
587       &  * FLUX%SW_DN_SURF_CLEAR_BAND(IBAND_PAR(JBAND),KIDIA:KFDIA)
588ENDDO
589
590! Compute effective broadband emissivity. This is only approximate -
591! due to spectral variations in emissivity, it is not in general
592! possible to provide a broadband emissivity that can reproduce the
593! upwelling surface flux given the downwelling flux and the skin
594! temperature.
595ZBLACK_BODY_NET_LW = PFLUX_LW_DN(KIDIA:KFDIA) &
596     &  - RSIGMA*PTEMPERATURE_SKIN(KIDIA:KFDIA)**4
597PEMIS_OUT(KIDIA:KFDIA) = PSPECTRALEMISS(KIDIA:KFDIA,1) ! Default value
598WHERE (ABS(ZBLACK_BODY_NET_LW) > 1.0E-5)
599  ! This calculation can go outside the range of any individual
600  ! spectral emissivity value, so needs to be capped
601  PEMIS_OUT(KIDIA:KFDIA) = MAX(0.8_JPRB, MIN(0.99_JPRB, PFLUX_LW(KIDIA:KFDIA,KLEV+1) / ZBLACK_BODY_NET_LW))
602ENDWHERE
603
604! Copy longwave derivatives
605IF (YRERAD%LAPPROXLWUPDATE) THEN
606  PLWDERIVATIVE(KIDIA:KFDIA,:) = FLUX%LW_DERIVATIVES(KIDIA:KFDIA,:)
607ENDIF
608
609! Store the shortwave downwelling fluxes in each albedo band
610IF (YRERAD%LAPPROXSWUPDATE) THEN
611  PSWDIFFUSEBAND(KIDIA:KFDIA,:) = TRANSPOSE(FLUX%SW_DN_DIFFUSE_SURF_CANOPY(:,KIDIA:KFDIA))
612  PSWDIRECTBAND (KIDIA:KFDIA,:) = TRANSPOSE(FLUX%SW_DN_DIRECT_SURF_CANOPY (:,KIDIA:KFDIA))
613ENDIF
614
615CALL SINGLE_LEVEL%DEALLOCATE
616CALL THERMODYNAMICS%DEALLOCATE
617CALL GAS%DEALLOCATE
618CALL YLCLOUD%DEALLOCATE
619CALL AEROSOL%DEALLOCATE
620CALL FLUX%DEALLOCATE
621
622END ASSOCIATE
623
624IF (LHOOK) CALL DR_HOOK('RADIATION_SCHEME',1,ZHOOK_HANDLE)
625
626END SUBROUTINE RADIATION_SCHEME
Note: See TracBrowser for help on using the repository browser.