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

Last change on this file since 4073 was 4031, checked in by idelkadi, 3 years ago

Corrections in the LMDZ-ECRAN interface:

  • gas concentrations in VMR instead of MMR
  • use of a constant value for each gas instead of a table of values
  • specify in the 14 SW bands of ECRAD those corresponding to the 6 SW bands for diffuse and direct albedo in LMDZ
  • passage of all processing of input fields for Ecrad (conversions, calculations, ...) from radiation.F90 to radlw.F90
File size: 30.0 KB
Line 
1! AI mars 2021
2! ====================== Interface between ECRAD and LMDZ ====================
3! radiation_scheme.F90 appelee dans radlwsw_m.F90 si iflag_rttm = 2
4! revoir toutes les parties avec "AI ATTENTION"
5! Mars 2021 :
6!             - Revoir toutes les parties commentees AI ATTENTION
7!             1. Traitement des aerosols
8!             2. Verifier les parametres times issus de LMDZ (calcul issed)
9!             3. Configuration a partir de namelist
10!             4. frac_std = 0.75     
11! ============================================================================
12
13SUBROUTINE RADIATION_SCHEME &
14! Inputs
15     & (KIDIA, KFDIA, KLON, KLEV, KAEROLMDZ, NSW, &
16     &  IDAY, TIME, &
17     &  PSOLAR_IRRADIANCE, &
18     &  PMU0, PTEMPERATURE_SKIN, PALBEDO_DIF, PALBEDO_DIR, &
19     &  PEMIS, PEMIS_WINDOW, &
20     &  PCCN_LAND, PCCN_SEA, &
21     &  PGELAM, PGEMU, PLAND_SEA_MASK, &
22     &  PPRESSURE, PTEMPERATURE, &
23     &  PPRESSURE_H, PTEMPERATURE_H, PQ, PQSAT, &
24     &  PCO2, PCH4, PN2O, PNO2, PCFC11, PCFC12, PHCFC22, &
25     &  PCCL4, PO3, PO2, &
26     &  PCLOUD_FRAC, PQ_LIQUID, PQ_ICE, PQ_RAIN, PQ_SNOW, &
27     &  ZRE_LIQUID_UM, ZRE_ICE_UM, &
28     &  PAEROSOL_OLD, PAEROSOL, &
29! Outputs
30     &  PFLUX_SW, PFLUX_LW, PFLUX_SW_CLEAR, PFLUX_LW_CLEAR, &
31     &  PFLUX_SW_DN, PFLUX_LW_DN, PFLUX_SW_DN_CLEAR, PFLUX_LW_DN_CLEAR, &
32     &  PFLUX_SW_UP, PFLUX_LW_UP, PFLUX_SW_UP_CLEAR, PFLUX_LW_UP_CLEAR, &
33     &  PFLUX_DIR, PFLUX_DIR_CLEAR, PFLUX_DIR_INTO_SUN, &
34     &  PFLUX_UV, PFLUX_PAR, PFLUX_PAR_CLEAR, &
35     &  PEMIS_OUT, PLWDERIVATIVE, &
36     &  PSWDIFFUSEBAND, PSWDIRECTBAND)
37
38
39! RADIATION_SCHEME - Interface to modular radiation scheme
40!
41! (C) Copyright 2015- ECMWF.
42!
43! This software is licensed under the terms of the Apache Licence Version 2.0
44! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
45!
46! In applying this licence, ECMWF does not waive the privileges and immunities
47! granted to it by virtue of its status as an intergovernmental organisation
48! nor does it submit to any jurisdiction.
49!
50! PURPOSE
51! -------
52!   The modular radiation scheme is contained in a separate
53!   library. This routine puts the the IFS arrays into appropriate
54!   objects, computing the additional data that is required, and sends
55!   it to the radiation scheme.  It returns net fluxes and surface
56!   flux components needed by the rest of the model.
57!
58!   Lower case is used for variables and types taken from the
59!   radiation library
60!
61! INTERFACE
62! ---------
63!    RADIATION_SCHEME is called from RADLSWR. The
64!    SETUP_RADIATION_SCHEME routine (in the RADIATION_SETUP module)
65!    should have been run first.
66!
67! AUTHOR
68! ------
69!   Robin Hogan, ECMWF
70!   Original: 2015-09-16
71!
72! MODIFICATIONS
73! -------------
74!
75! TO DO
76! -----
77!
78!-----------------------------------------------------------------------
79
80! Modules from ifs or ifsaux libraries
81USE PARKIND1 , ONLY : JPIM, JPRB
82USE YOMHOOK  , ONLY : LHOOK, DR_HOOK
83! AI ATTENTION module propre a ifs commentes
84!USE YOERAD   , ONLY : YRERAD
85USE RADIATION_SETUP, ONLY : SETUP_RADIATION_SCHEME, rad_config, &
86!USE RADIATION_SETUP, ONLY : &
87     &  NWEIGHT_UV,  IBAND_UV,  WEIGHT_UV, &
88     &  NWEIGHT_PAR, IBAND_PAR, WEIGHT_PAR, &
89     &  ITYPE_TROP_BG_AER,  TROP_BG_AER_MASS_EXT, &
90     &  ITYPE_STRAT_BG_AER, STRAT_BG_AER_MASS_EXT
91! Commentes : jour, date de la simulation
92!USE YOMRIP0  , ONLY : NINDAT
93!USE YOMCT3   , ONLY : NSTEP
94!USE YOMRIP   , ONLY : YRRIP
95USE YOMCST   , ONLY : RSIGMA ! Stefan-Boltzmann constant
96
97! Modules from radiation library
98! AI ATTENTION
99!use radiation_config,         only : config_type
100USE radiation_single_level,   ONLY : single_level_type
101USE radiation_thermodynamics, ONLY : thermodynamics_type
102USE radiation_gas
103USE radiation_cloud,          ONLY : cloud_type
104USE radiation_aerosol,        ONLY : aerosol_type
105USE radiation_flux,           ONLY : flux_type
106USE radiation_interface,      ONLY : radiation, set_gas_units
107USE radiation_save,           ONLY : save_inputs
108
109IMPLICIT NONE
110
111! INPUT ARGUMENTS
112! *** Array dimensions and ranges
113INTEGER(KIND=JPIM),INTENT(IN) :: KIDIA    ! Start column to process
114INTEGER(KIND=JPIM),INTENT(IN) :: KFDIA    ! End column to process
115!INTEGER, INTENT(IN) :: KIDIA, KFDIA
116INTEGER(KIND=JPIM),INTENT(IN) :: KLON     ! Number of columns
117INTEGER(KIND=JPIM),INTENT(IN) :: KLEV     ! Number of levels
118!INTEGER, INTENT(IN) :: KLON, KLEV
119INTEGER(KIND=JPIM),INTENT(IN) :: KAEROLMDZ ! Number of aerosol types
120INTEGER(KIND=JPIM),INTENT(IN) :: NSW ! Numbe of bands
121
122! AI ATTENTION
123INTEGER, PARAMETER :: KAEROSOL = 12
124
125! *** Single-level fields
126REAL(KIND=JPRB),   INTENT(IN) :: PSOLAR_IRRADIANCE ! (W m-2)
127REAL(KIND=JPRB),   INTENT(IN) :: PMU0(KLON) ! Cosine of solar zenith ang
128REAL(KIND=JPRB),   INTENT(IN) :: PTEMPERATURE_SKIN(KLON) ! (K)
129! Diffuse and direct components of surface shortwave albedo
130!REAL(KIND=JPRB),   INTENT(IN) :: PALBEDO_DIF(KLON,YRERAD%NSW)
131!REAL(KIND=JPRB),   INTENT(IN) :: PALBEDO_DIR(KLON,YRERAD%NSW)
132REAL(KIND=JPRB),   INTENT(IN) :: PALBEDO_DIF(KLON,NSW)
133REAL(KIND=JPRB),   INTENT(IN) :: PALBEDO_DIR(KLON,NSW)
134! Longwave emissivity outside and inside the window region
135REAL(KIND=JPRB),   INTENT(IN) :: PEMIS(KLON)
136REAL(KIND=JPRB),   INTENT(IN) :: PEMIS_WINDOW(KLON)
137! Longitude (radians), sine of latitude
138REAL(KIND=JPRB),   INTENT(IN) :: PGELAM(KLON)
139REAL(KIND=JPRB),   INTENT(IN) :: PGEMU(KLON)
140! Land-sea mask
141REAL(KIND=JPRB),   INTENT(IN) :: PLAND_SEA_MASK(KLON)
142
143! *** Variables on full levels
144REAL(KIND=JPRB),   INTENT(IN) :: PPRESSURE(KLON,KLEV)    ! (Pa)
145REAL(KIND=JPRB),   INTENT(IN) :: PTEMPERATURE(KLON,KLEV) ! (K)
146! *** Variables on half levels
147REAL(KIND=JPRB),   INTENT(IN) :: PPRESSURE_H(KLON,KLEV+1)    ! (Pa)
148REAL(KIND=JPRB),   INTENT(IN) :: PTEMPERATURE_H(KLON,KLEV+1) ! (K)
149
150! *** Gas mass mixing ratios on full levels
151REAL(KIND=JPRB),   INTENT(IN) :: PQ(KLON,KLEV)
152! AI
153REAL(KIND=JPRB),   INTENT(IN) :: PQSAT(KLON,KLEV)
154REAL(KIND=JPRB),   INTENT(IN) :: PCO2
155REAL(KIND=JPRB),   INTENT(IN) :: PCH4
156REAL(KIND=JPRB),   INTENT(IN) :: PN2O
157REAL(KIND=JPRB),   INTENT(IN) :: PNO2
158REAL(KIND=JPRB),   INTENT(IN) :: PCFC11
159REAL(KIND=JPRB),   INTENT(IN) :: PCFC12
160REAL(KIND=JPRB),   INTENT(IN) :: PHCFC22
161REAL(KIND=JPRB),   INTENT(IN) :: PCCL4
162REAL(KIND=JPRB),   INTENT(IN) :: PO3(KLON,KLEV) ! AI (kg/kg) ATTENTION (Pa*kg/kg)
163REAL(KIND=JPRB),   INTENT(IN) :: PO2
164
165! *** Cloud fraction and hydrometeor mass mixing ratios
166REAL(KIND=JPRB),   INTENT(IN) :: PCLOUD_FRAC(KLON,KLEV)
167REAL(KIND=JPRB),   INTENT(IN) :: PQ_LIQUID(KLON,KLEV)
168REAL(KIND=JPRB),   INTENT(IN) :: PQ_ICE(KLON,KLEV)
169REAL(KIND=JPRB),   INTENT(IN) :: PQ_RAIN(KLON,KLEV)
170REAL(KIND=JPRB),   INTENT(IN) :: PQ_SNOW(KLON,KLEV)
171
172! *** Aerosol mass mixing ratios
173REAL(KIND=JPRB),   INTENT(IN) :: PAEROSOL_OLD(KLON,6,KLEV)
174REAL(KIND=JPRB),   INTENT(IN) :: PAEROSOL(KLON,KLEV,KAEROSOL)
175
176REAL(KIND=JPRB),   INTENT(IN) :: PCCN_LAND(KLON)
177REAL(KIND=JPRB),   INTENT(IN) :: PCCN_SEA(KLON)
178
179!AI mars 2021
180INTEGER(KIND=JPIM), INTENT(IN) :: IDAY
181REAL(KIND=JPRB), INTENT(IN)    :: TIME
182
183
184! OUTPUT ARGUMENTS
185
186! *** Net fluxes on half-levels (W m-2)
187REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW(KLON,KLEV+1)
188REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW(KLON,KLEV+1)
189REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_CLEAR(KLON,KLEV+1)
190REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_CLEAR(KLON,KLEV+1)
191
192!*** DN and UP flux on half-levels (W m-2)
193REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_DN(KLON,KLEV+1)
194REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_DN(KLON,KLEV+1)
195REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_DN_CLEAR(KLON,KLEV+1)
196REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_DN_CLEAR(KLON,KLEV+1)
197REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_UP(KLON,KLEV+1)
198REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_UP(KLON,KLEV+1)
199REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_UP_CLEAR(KLON,KLEV+1)
200REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_UP_CLEAR(KLON,KLEV+1)
201
202! *** Surface flux components (W m-2)
203! AI ATTENTION
204!REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_DN_SURF(KLON)
205!REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_DN_SURF(KLON)
206!REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_DN_CLEAR_SURF(KLON)
207!REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_DN_CLEAR_SURF(KLON)
208!REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_UP_SURF(KLON)
209!REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_UP_SURF(KLON)
210!REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_UP_CLEAR_SURF(KLON)
211!REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_UP_CLEAR_SURF(KLON)
212
213! Direct component of surface flux into horizontal plane
214REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_DIR(KLON)
215REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_DIR_CLEAR(KLON)
216! As PFLUX_DIR but into a plane perpendicular to the sun
217REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_DIR_INTO_SUN(KLON)
218
219! *** Ultraviolet and photosynthetically active radiation (W m-2)
220REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_UV(KLON)
221REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_PAR(KLON)
222REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_PAR_CLEAR(KLON)
223
224! *** Other single-level diagnostics
225! Top-of-atmosphere incident solar flux (W m-2)
226! AI ATTENTION
227!REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_DN_TOA(KLON)
228!REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_UP_TOA(KLON)
229!REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_DN_TOA(KLON)
230!REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_UP_TOA(KLON)
231
232! Diagnosed longwave surface emissivity across the whole spectrum
233REAL(KIND=JPRB),  INTENT(OUT) :: PEMIS_OUT(KLON)   
234
235! Partial derivative of total-sky longwave upward flux at each level
236! with respect to upward flux at surface, used to correct heating
237! rates at gridpoints/timesteps between calls to the full radiation
238! scheme.  Note that this version uses the convention of level index
239! increasing downwards, unlike the local variable ZLwDerivative that
240! is returned from the LW radiation scheme.
241REAL(KIND=JPRB),  INTENT(OUT) :: PLWDERIVATIVE(KLON,KLEV+1)
242
243! Surface diffuse and direct downwelling shortwave flux in each
244! shortwave albedo band, used in RADINTG to update the surface fluxes
245! accounting for high-resolution albedo information
246REAL(KIND=JPRB),  INTENT(OUT) :: PSWDIFFUSEBAND(KLON,NSW)
247REAL(KIND=JPRB),  INTENT(OUT) :: PSWDIRECTBAND (KLON,NSW)
248
249! LOCAL VARIABLES
250! AI ATTENTION
251!type(config_type)         :: rad_config
252TYPE(single_level_type)   :: single_level
253TYPE(thermodynamics_type) :: thermodynamics
254TYPE(gas_type)            :: gas
255TYPE(cloud_type)          :: cloud
256TYPE(aerosol_type)        :: aerosol
257TYPE(flux_type)           :: flux
258
259! Mass mixing ratio of ozone (kg/kg)
260REAL(KIND=JPRB)           :: ZO3(KLON,KLEV)
261
262! Cloud effective radii in microns
263REAL(KIND=JPRB)           :: ZRE_LIQUID_UM(KLON,KLEV)
264REAL(KIND=JPRB)           :: ZRE_ICE_UM(KLON,KLEV)
265
266! Cloud overlap decorrelation length for cloud boundaries in km
267REAL(KIND=JPRB)           :: ZDECORR_LEN_KM(KLON)
268
269! Ratio of cloud overlap decorrelation length for cloud water
270! inhomogeneities to that for cloud boundaries (typically 0.5)
271REAL(KIND=JPRB)           :: ZDECORR_LEN_RATIO
272
273! The surface net longwave flux if the surface was a black body, used
274! to compute the effective broadband surface emissivity
275REAL(KIND=JPRB)           :: ZBLACK_BODY_NET_LW(KIDIA:KFDIA)
276
277! Layer mass in kg m-2
278REAL(KIND=JPRB)           :: ZLAYER_MASS(KIDIA:KFDIA,KLEV)
279
280! Time integers
281INTEGER :: ITIM
282
283! Loop indices
284INTEGER :: JLON, JLEV, JBAND, JB_ALBEDO, JAER
285
286REAL(KIND=JPRB) :: ZHOOK_HANDLE
287
288! AI ATTENTION traitement aerosols
289INTEGER, PARAMETER :: NAERMACC = 1
290
291! AI ATTENTION
292real(jprb), parameter    :: frac_std = 0.75
293
294! Name of file names specified on command line
295character(len=512) :: file_name
296
297logical :: loutput=.true.
298logical :: lprint_input=.false.
299logical :: lprint_config=.true.
300
301! Import time functions for iseed calculation
302! AI ATTENTION propre a ifs
303!#include "fcttim.func.h"
304!#include "liquid_effective_radius.intfb.h"
305!#include "ice_effective_radius.intfb.h"
306!#include "cloud_overlap_decorr_len.intfb.h"
307!#include "satur.intfb.h"
308
309! Verifier les inputs
310print*,'=============== dans radiation_scheme : ==================='
311if (lprint_input) then
312  print*,'********** Verification des entrees *************'
313  print*,'KIDIA, KFDIA, KLON, KLEV, KAEROLMDZ, NSW =', &
314        KIDIA, KFDIA, KLON, KLEV, KAEROLMDZ, NSW
315  print*,'IDAY, TIME =', IDAY, TIME
316  print*,'PSOLAR_IRRADIANCE =', PSOLAR_IRRADIANCE
317  print*,'PMU0 =', PMU0
318  print*,'PTEMPERATURE_SKIN =',PTEMPERATURE_SKIN
319  print*,'PEMIS, PEMIS_WINDOW =', PEMIS, PEMIS_WINDOW
320  print*,'PCCN_LAND, PCCN_SEA =', PCCN_LAND, PCCN_SEA
321  print*,'PGELAM, PGEMU =', PGELAM, PGEMU
322  print*,'PPRESSURE =', PPRESSURE
323  print*,'PTEMPERATURE =', PTEMPERATURE
324  print*,'PPRESSURE_H =', PPRESSURE_H
325  print*,'PTEMPERATURE_H =', PTEMPERATURE_H
326  print*,'PQ =', PQ
327  print*,'PQSAT=',PQSAT
328  print*,'PCO2, PCH4, PN2O, PNO2, PCFC11, PCFC12, PHCFC22, PCCL4 =', &
329        PCO2, PCH4, PN2O, PNO2, PCFC11, PCFC12, PHCFC22, PCCL4
330  print*,'PO3 =',PO3
331  print*,'PCLOUD_FRAC, PQ_LIQUID, PQ_ICE, PQ_RAIN, PQ_SNOW =', &
332        PCLOUD_FRAC, PQ_LIQUID, PQ_ICE, PQ_RAIN, PQ_SNOW
333  print*,'ZRE_LIQUID_UM, ZRE_ICE_UM =', &
334        ZRE_LIQUID_UM, ZRE_ICE_UM
335  print*,'PAEROSOL_OLD, PAEROSOL =', PAEROSOL_OLD, PAEROSOL
336endif
337! AI ATTENTION lecture de namelist
338! alternative a l appel de radiation_setup ifs
339!file_name="namelist_ecrad"
340!call rad_config%read(file_name=file_name)
341! Setup the radiation scheme: load the coefficients for gas and
342! cloud optics, currently from RRTMG
343!call setup_radiation(rad_config)
344
345IF (LHOOK) CALL DR_HOOK('RADIATION_SCHEME',0,ZHOOK_HANDLE)
346print*,'Entree dans radiation_scheme'
347! AI appel radiation_setup
348call SETUP_RADIATION_SCHEME(loutput)
349!! Les 6 bandes SW pour l'albedo :
350!! 0.185-0.25, 0.25-0.4, 0.4-0.69 , 0.69-1.19, 1.19-2.38, 2.38-4.00 micro-metre
351  call rad_config%define_sw_albedo_intervals(6, &
352             &  (/ 0.25e-6_jprb, 0.44e-6_jprb, 1.19e-6_jprb, &
353             &     2.38e-6_jprb, 4.00e-6_jprb /),  (/ 1,2,3,4,5,6 /))
354
355if (lprint_config) then
356  print*,'************* Parametres de configuration  ********************'
357  print*,'rad_config%iverbosesetup = ',rad_config%iverbosesetup
358  print*,'rad_config%iverbose = ',rad_config%iverbose
359  print*,'rad_config%directory_name =', rad_config%directory_name
360  print*,'rad_config%do_lw_derivatives =',rad_config%do_lw_derivatives
361  print*,'rad_config%do_surface_sw_spectral_flux =', &
362        rad_config%do_surface_sw_spectral_flux
363  print*,'rad_config%do_setup_ifsrrtm =', rad_config%do_setup_ifsrrtm
364  print*,'rad_config%i_liq_model =',rad_config%i_liq_model
365  print*,'rad_config%i_ice_model =',rad_config%i_ice_model
366  print*,'rad_config%i_overlap_scheme =', rad_config%i_overlap_scheme
367  print*,'rad_config%use_aerosols = ', rad_config%use_aerosols
368  print*,'rad_config%n_aerosol_types = ', rad_config%n_aerosol_types
369  print*,'rad_config%i_solver_lw =',rad_config%i_solver_lw
370  print*,'rad_config%i_solver_sw =',rad_config%i_solver_sw
371  print*,'rad_config%do_3d_effects =', rad_config%do_3d_effects
372  print*,'rad_config%do_sw_delta_scaling_with_gases =', &
373       rad_config%do_sw_delta_scaling_with_gases
374  print*,'rad_config%do_lw_aerosol_scattering =', &
375       rad_config%do_lw_aerosol_scattering
376  print*,'rad_config%i_albedo_from_band_sw = ', &
377       rad_config%i_albedo_from_band_sw
378  print*,'n_bands_lw =', rad_config%n_bands_lw
379  print*,'rad_config%i_emiss_from_band_lw =', rad_config%i_emiss_from_band_lw
380endif
381
382! AI ATTENTION
383! Allocate memory in radiation objects
384CALL single_level%allocate(KLON, NSW, 1, &
385     &                     use_sw_albedo_direct=.TRUE.)
386
387print*,'************* THERMO (allocate + input) ************************************'
388! Set thermodynamic profiles: simply copy over the half-level
389! pressure and temperature
390print*,'Appel allocate thermo'
391CALL thermodynamics%allocate(KLON, KLEV, use_h2o_sat=.true.)
392print*,'Definir les champs thermo'
393! AI
394! pressure_hl > paprs
395! temperature_hl calculee dans radlsw de la meme facon que pour RRTM
396thermodynamics%pressure_hl   (KIDIA:KFDIA,:) = PPRESSURE_H   (KIDIA:KFDIA,:)
397thermodynamics%temperature_hl(KIDIA:KFDIA,:) = PTEMPERATURE_H(KIDIA:KFDIA,:)
398
399! IFS currently sets the half-level temperature at the surface to be
400! equal to the skin temperature. The radiation scheme takes as input
401! only the half-level temperatures and assumes the Planck function to
402! vary linearly in optical depth between half levels. In the lowest
403! atmospheric layer, where the atmospheric temperature can be much
404! cooler than the skin temperature, this can lead to significant
405! differences between the effective temperature of this lowest layer
406! and the true value in the model.
407!
408! We may approximate the temperature profile in the lowest model level
409! as piecewise linear between the top of the layer T[k-1/2], the
410! centre of the layer T[k] and the base of the layer Tskin.  The mean
411! temperature of the layer is then 0.25*T[k-1/2] + 0.5*T[k] +
412! 0.25*Tskin, which can be achieved by setting the atmospheric
413! temperature at the half-level corresponding to the surface as
414! follows:
415!thermodynamics%temperature_hl(KIDIA:KFDIA,KLEV+1) &
416!     &  = PTEMPERATURE(KIDIA:KFDIA,KLEV) &
417!     &  + 0.5_JPRB * (PTEMPERATURE_H(KIDIA:KFDIA,KLEV+1) &
418!     &               -PTEMPERATURE_H(KIDIA:KFDIA,KLEV))
419
420! Alternatively we respect the model's atmospheric temperature in the
421! lowest model level by setting the temperature at the lowest
422! half-level such that the mean temperature of the layer is correct:
423!thermodynamics%temperature_hl(KIDIA:KFDIA,KLEV+1) &
424!     &  = 2.0_JPRB * PTEMPERATURE(KIDIA:KFDIA,KLEV) &
425!     &             - PTEMPERATURE_H(KIDIA:KFDIA,KLEV)
426
427! Compute saturation specific humidity, used to hydrate aerosols. The
428! "2" for the last argument indicates that the routine is not being
429! called from within the convection scheme.
430!CALL SATUR(KIDIA, KFDIA, KLON, 1, KLEV, &
431!     &  PPRESSURE, PTEMPERATURE, thermodynamics%h2o_sat_liq, 2) 
432! Alternative approximate version using temperature and pressure from
433! the thermodynamics structure
434print*,'Compute saturation specific humidity'
435CALL thermodynamics%calc_saturation_wrt_liquid(KIDIA, KFDIA)
436
437print*,'********** SINGLE LEVEL VARS **********************************'
438!AI ATTENTION
439!thermodynamics%h2o_sat_liq = PQSAT
440! Set single-level fileds
441single_level%solar_irradiance              = PSOLAR_IRRADIANCE
442!allocate(single_level%cos_sza(KIDIA:KFDIA))
443single_level%cos_sza(KIDIA:KFDIA)          = PMU0(KIDIA:KFDIA)
444!allocate(single_level%skin_temperature(KIDIA:KFDIA))
445single_level%skin_temperature(KIDIA:KFDIA) = PTEMPERATURE_SKIN(KIDIA:KFDIA)
446!allocate(single_level%sw_albedo(KIDIA:KFDIA,1))
447single_level%sw_albedo(KIDIA:KFDIA,:)      = PALBEDO_DIF(KIDIA:KFDIA,:)
448!single_level%sw_albedo(KIDIA:KFDIA,:)      = 0.2_JPRB
449!allocate(single_level%sw_albedo_direct(KIDIA:KFDIA,1))
450single_level%sw_albedo_direct(KIDIA:KFDIA,:)=PALBEDO_DIR(KIDIA:KFDIA,:)
451!single_level%sw_albedo_direct(KIDIA:KFDIA,:)=0.2_JPRB
452! Longwave emissivity is in two bands
453!allocate(single_level%lw_emissivity(KIDIA:KFDIA,1))
454!single_level%lw_emissivity(KIDIA:KFDIA,1)  = 1.0_JPRB
455single_level%lw_emissivity(KIDIA:KFDIA,1)  = PEMIS(KIDIA:KFDIA)
456!single_level%lw_emissivity(KIDIA:KFDIA,2)  = PEMIS_WINDOW(KIDIA:KFDIA)
457
458! Create the relevant seed from date and time get the starting day
459! and number of minutes since start
460!IDAY = NDD(NINDAT)
461!cur_day
462!ITIM = NINT(NSTEP * YRRIP%TSTEP / 60.0_JPRB)
463ITIM = NINT(TIME / 60.0_JPRB)
464!current_time
465!allocate(single_level%iseed(KIDIA:KFDIA))
466DO JLON = KIDIA, KFDIA
467  ! This method gives a unique value for roughly every 1-km square
468  ! on the globe and every minute.  ASIN(PGEMU)*60 gives rough
469  ! latitude in degrees, which we multiply by 100 to give a unique
470  ! value for roughly every km. PGELAM*60*100 gives a unique number
471  ! for roughly every km of longitude around the equator, which we
472  ! multiply by 180*100 so there is no overlap with the latitude
473  ! values.  The result can be contained in a 32-byte integer (but
474  ! since random numbers are generated with the help of integer
475  ! overflow, it should not matter if the number did overflow).
476  single_level%iseed(JLON) = ITIM + IDAY &
477       &  +  NINT(PGELAM(JLON)*108000000.0_JPRB &
478       &          + ASIN(PGEMU(JLON))*6000.0_JPRB)
479ENDDO
480
481print*,'********** CLOUDS (allocate + input) *******************************************'
482print*,'Appel Allocate clouds'
483CALL cloud%allocate(KLON, KLEV)
484! Set cloud fields
485cloud%q_liq(KIDIA:KFDIA,:)    = PQ_LIQUID(KIDIA:KFDIA,:)
486cloud%q_ice(KIDIA:KFDIA,:)    = PQ_ICE(KIDIA:KFDIA,:) + PQ_SNOW(KIDIA:KFDIA,:)
487cloud%fraction(KIDIA:KFDIA,:) = PCLOUD_FRAC(KIDIA:KFDIA,:)
488
489! Compute effective radii and convert to metres
490!CALL LIQUID_EFFECTIVE_RADIUS(KIDIA, KFDIA, KLON, KLEV, &
491!     &  PPRESSURE, PTEMPERATURE, PCLOUD_FRAC, PQ_LIQUID, PQ_RAIN, &
492!     &  PLAND_SEA_MASK, PCCN_LAND, PCCN_SEA, &
493!     &  ZRE_LIQUID_UM)
494cloud%re_liq(KIDIA:KFDIA,:) = ZRE_LIQUID_UM(KIDIA:KFDIA,:)
495
496!CALL ICE_EFFECTIVE_RADIUS(KIDIA, KFDIA, KLON, KLEV, &
497!     &  PPRESSURE, PTEMPERATURE, PCLOUD_FRAC, PQ_ICE, PQ_SNOW, PGEMU, &
498!     &  ZRE_ICE_UM)
499cloud%re_ice(KIDIA:KFDIA,:) = ZRE_ICE_UM(KIDIA:KFDIA,:)
500
501! Get the cloud overlap decorrelation length (for cloud boundaries),
502! in km, according to the parameterization specified by NDECOLAT,
503! and insert into the "cloud" object. Also get the ratio of
504! decorrelation lengths for cloud water content inhomogeneities and
505! cloud boundaries, and set it in the "rad_config" object.
506!CALL CLOUD_OVERLAP_DECORR_LEN(KIDIA, KFDIA, KLON, PGEMU, YRERAD%NDECOLAT, &
507!     &    ZDECORR_LEN_KM, PDECORR_LEN_RATIO=ZDECORR_LEN_RATIO)
508! AI ATTENTION a revoir
509ZDECORR_LEN_RATIO = 0.5_JPRB
510rad_config%cloud_inhom_decorr_scaling = ZDECORR_LEN_RATIO
511ZDECORR_LEN_KM = 2000.0_JPRB
512DO JLON = KIDIA,KFDIA
513  CALL cloud%set_overlap_param(thermodynamics, &
514       &                       ZDECORR_LEN_KM(JLON), &
515       &                       istartcol=JLON, iendcol=JLON)
516ENDDO
517
518! Cloud water content fractional standard deviation is configurable
519! from namelist NAERAD but must be globally constant. Before it was
520! hard coded at 1.0.
521!CALL cloud%create_fractional_std(KLON, KLEV, YRERAD%RCLOUD_FRAC_STD)
522CALL cloud%create_fractional_std(KLON, KLEV, frac_std)
523
524! By default mid and high cloud effective size is 10 km
525CALL cloud%create_inv_cloud_effective_size(KLON,KLEV,1.0_JPRB/10000.0_JPRB)
526! But for boundary clouds (eta > 0.8) we set it to 1 km
527DO JLEV = 1,KLEV
528  DO JLON = KIDIA,KFDIA
529    IF (PPRESSURE(JLON,JLEV) > 0.8_JPRB * PPRESSURE_H(JLON,KLEV+1)) THEN
530      cloud%inv_cloud_effective_size(JLON,JLEV) = 1.0e-3_JPRB
531    ENDIF
532  ENDDO
533ENDDO
534!AI ATTENTION meme traitement dans le version offline
535!call cloud%create_inv_cloud_effective_size_eta(ncol, nlev, &
536!               &  thermodynamics%pressure_hl, &
537!               &  low_inv_effective_size, &
538!               &  middle_inv_effective_size, &
539!               &  high_inv_effective_size, 0.8_jprb, 0.45_jprb)     
540
541print*,'******** AEROSOLS (allocate + input) **************************************'
542IF (NAERMACC > 0) THEN
543  CALL aerosol%allocate(KLON, 1, KLEV, KAEROSOL) ! MACC climatology
544ELSE
545  CALL aerosol%allocate(KLON, 1, KLEV, 6) ! Tegen climatology
546ENDIF
547! Compute the dry mass of each layer neglecting humidity effects, in
548! kg m-2, needed to scale some of the aerosol inputs
549! AI commente ATTENTION
550!CALL thermodynamics%get_layer_mass(ZLAYER_MASS)
551
552! Copy over aerosol mass mixing ratio
553!IF (NAERMACC > 0) THEN
554
555  ! MACC aerosol climatology - this is already in mass mixing ratio
556  ! units with the required array orientation so we can copy it over
557  ! directly
558  aerosol%mixing_ratio(KIDIA:KFDIA,:,:) = PAEROSOL(KIDIA:KFDIA,:,:)
559
560  ! Add the tropospheric and stratospheric backgrounds contained in the
561  ! old Tegen arrays - this is very ugly!
562! AI ATTENTION
563!  IF (TROP_BG_AER_MASS_EXT > 0.0_JPRB) THEN
564!    aerosol%mixing_ratio(KIDIA:KFDIA,:,ITYPE_TROP_BG_AER) &
565!         &  = aerosol%mixing_ratio(KIDIA:KFDIA,:,ITYPE_TROP_BG_AER) &
566!         &  + PAEROSOL_OLD(KIDIA:KFDIA,1,:) &
567!         &  / (ZLAYER_MASS * TROP_BG_AER_MASS_EXT)
568!  ENDIF
569!  IF (STRAT_BG_AER_MASS_EXT > 0.0_JPRB) THEN
570!    aerosol%mixing_ratio(KIDIA:KFDIA,:,ITYPE_STRAT_BG_AER) &
571!         &  = aerosol%mixing_ratio(KIDIA:KFDIA,:,ITYPE_STRAT_BG_AER) &
572!         &  + PAEROSOL_OLD(KIDIA:KFDIA,6,:) &
573!         &  / (ZLAYER_MASS * STRAT_BG_AER_MASS_EXT)
574!  ENDIF
575
576!ELSE
577
578  ! Tegen aerosol climatology - the array PAEROSOL_OLD contains the
579  ! 550-nm optical depth in each layer. The optics data file
580  ! aerosol_ifs_rrtm_tegen.nc does not contain mass extinction
581  ! coefficient, but a scaling factor that the 550-nm optical depth
582  ! should be multiplied by to obtain the optical depth in each
583  ! spectral band.  Therefore, in order for the units to work out, we
584  ! need to divide by the layer mass (in kg m-2) to obtain the 550-nm
585  ! cross-section per unit mass of dry air (so in m2 kg-1).  We also
586  ! need to permute the array.
587!  DO JLEV = 1,KLEV
588!    DO JAER = 1,6
589!      aerosol%mixing_ratio(KIDIA:KFDIA,JLEV,JAER) &
590!         &  = PAEROSOL_OLD(KIDIA:KFDIA,JAER,JLEV) &
591!         &  / ZLAYER_MASS(KIDIA:KFDIA,JLEV)
592!    ENDDO
593!  ENDDO
594
595!ENDIF
596
597print*,'********** GAS (allocate + input) ************************************************'
598print*,'Appel Allocate gas'
599CALL gas%allocate(KLON, KLEV)
600
601! Convert ozone Pa*kg/kg to kg/kg
602! AI ATTENTION
603!DO JLEV = 1,KLEV
604!  DO JLON = KIDIA,KFDIA
605!    ZO3(JLON,JLEV) = PO3_DP(JLON,JLEV) &
606!         &         / (PPRESSURE_H(JLON,JLEV+1)-PPRESSURE_H(JLON,JLEV))
607!  ENDDO
608!ENDDO
609
610!  Insert gas mixing ratios
611print*,'Insert gas mixing ratios'
612CALL gas%put(IH2O,    IMassMixingRatio, PQ)
613CALL gas%put(IO3,     IMassMixingRatio, PO3)
614CALL gas%put_well_mixed(ICO2,    IVolumeMixingRatio, PCO2)
615CALL gas%put_well_mixed(ICH4,    IVolumeMixingRatio, PCH4)
616CALL gas%put_well_mixed(IN2O,    IVolumeMixingRatio, PN2O)
617CALL gas%put_well_mixed(ICFC11,  IVolumeMixingRatio, PCFC11)
618CALL gas%put_well_mixed(ICFC12,  IVolumeMixingRatio, PCFC12)
619CALL gas%put_well_mixed(IHCFC22, IVolumeMixingRatio, PHCFC22)
620CALL gas%put_well_mixed(ICCL4,   IVolumeMixingRatio, PCCL4)
621CALL gas%put_well_mixed(IO2,     IVolumeMixingRatio, PO2)
622! Ensure the units of the gas mixing ratios are what is required by
623! the gas absorption model
624call set_gas_units(rad_config, gas)
625
626print*,'************** FLUX (allocate) ***********************'
627CALL flux%allocate(rad_config, 1, KLON, KLEV)
628
629! Call radiation scheme
630print*,'******** Appel radiation scheme **************************'
631CALL radiation(KLON, KLEV, KIDIA, KFDIA, rad_config, &
632     &  single_level, thermodynamics, gas, cloud, aerosol, flux)
633
634! Compute required output fluxes
635! DN and UP flux
636PFLUX_SW_DN(KIDIA:KFDIA,:) = flux%sw_dn(KIDIA:KFDIA,:)
637PFLUX_SW_UP(KIDIA:KFDIA,:) = flux%sw_up(KIDIA:KFDIA,:)
638PFLUX_LW_DN(KIDIA:KFDIA,:) = flux%lw_dn(KIDIA:KFDIA,:)
639PFLUX_LW_UP(KIDIA:KFDIA,:) = flux%lw_up(KIDIA:KFDIA,:)
640PFLUX_SW_DN_CLEAR(KIDIA:KFDIA,:) = flux%sw_dn_clear(KIDIA:KFDIA,:)
641PFLUX_SW_UP_CLEAR(KIDIA:KFDIA,:) = flux%sw_up_clear(KIDIA:KFDIA,:)
642PFLUX_LW_DN_CLEAR(KIDIA:KFDIA,:) = flux%lw_dn_clear(KIDIA:KFDIA,:)
643PFLUX_LW_UP_CLEAR(KIDIA:KFDIA,:) = flux%lw_up_clear(KIDIA:KFDIA,:)
644
645! First the net fluxes
646PFLUX_SW(KIDIA:KFDIA,:) = flux%sw_dn(KIDIA:KFDIA,:) - flux%sw_up(KIDIA:KFDIA,:)
647PFLUX_LW(KIDIA:KFDIA,:) = flux%lw_dn(KIDIA:KFDIA,:) - flux%lw_up(KIDIA:KFDIA,:)
648PFLUX_SW_CLEAR(KIDIA:KFDIA,:) &
649     &  = flux%sw_dn_clear(KIDIA:KFDIA,:) - flux%sw_up_clear(KIDIA:KFDIA,:)
650PFLUX_LW_CLEAR(KIDIA:KFDIA,:) &
651     &  = flux%lw_dn_clear(KIDIA:KFDIA,:) - flux%lw_up_clear(KIDIA:KFDIA,:)
652
653! Now the surface fluxes
654!PFLUX_SW_DN_SURF(KIDIA:KFDIA) = flux%sw_dn(KIDIA:KFDIA,KLEV+1)
655!PFLUX_LW_DN_SURF(KIDIA:KFDIA) = flux%lw_dn(KIDIA:KFDIA,KLEV+1)
656!PFLUX_SW_UP_SURF(KIDIA:KFDIA) = flux%sw_up(KIDIA:KFDIA,KLEV+1)
657!PFLUX_LW_UP_SURF(KIDIA:KFDIA) = flux%lw_up(KIDIA:KFDIA,KLEV+1)
658!PFLUX_SW_DN_CLEAR_SURF(KIDIA:KFDIA) = flux%sw_dn_clear(KIDIA:KFDIA,KLEV+1)
659!PFLUX_LW_DN_CLEAR_SURF(KIDIA:KFDIA) = flux%lw_dn_clear(KIDIA:KFDIA,KLEV+1)
660!PFLUX_SW_UP_CLEAR_SURF(KIDIA:KFDIA) = flux%sw_up_clear(KIDIA:KFDIA,KLEV+1)
661!PFLUX_LW_UP_CLEAR_SURF(KIDIA:KFDIA) = flux%lw_up_clear(KIDIA:KFDIA,KLEV+1)
662PFLUX_DIR(KIDIA:KFDIA) = flux%sw_dn_direct(KIDIA:KFDIA,KLEV+1)
663PFLUX_DIR_CLEAR(KIDIA:KFDIA) = flux%sw_dn_direct_clear(KIDIA:KFDIA,KLEV+1)
664PFLUX_DIR_INTO_SUN(KIDIA:KFDIA) = 0.0_JPRB
665WHERE (PMU0(KIDIA:KFDIA) > EPSILON(1.0_JPRB))
666  PFLUX_DIR_INTO_SUN(KIDIA:KFDIA) = PFLUX_DIR(KIDIA:KFDIA) / PMU0(KIDIA:KFDIA)
667END WHERE
668
669! Top-of-atmosphere downwelling flux
670!PFLUX_SW_DN_TOA(KIDIA:KFDIA) = flux%sw_dn(KIDIA:KFDIA,1)
671!PFLUX_SW_UP_TOA(KIDIA:KFDIA) = flux%sw_up(KIDIA:KFDIA,1)
672!PFLUX_LW_DN_TOA(KIDIA:KFDIA) = flux%lw_dn(KIDIA:KFDIA,1)
673!PFLUX_LW_UP_TOA(KIDIA:KFDIA) = flux%lw_up(KIDIA:KFDIA,1)
674
675! Compute UV fluxes as weighted sum of appropriate shortwave bands
676PFLUX_UV       (KIDIA:KFDIA) = 0.0_JPRB
677! AI ATTENTION
678!DO JBAND = 1,NWEIGHT_UV
679!  PFLUX_UV(KIDIA:KFDIA) = PFLUX_UV(KIDIA:KFDIA) + WEIGHT_UV(JBAND) &
680!       &  * flux%sw_dn_surf_band(IBAND_UV(JBAND),KIDIA:KFDIA)
681!ENDDO
682
683! Compute photosynthetically active radiation similarly
684PFLUX_PAR      (KIDIA:KFDIA) = 0.0_JPRB
685PFLUX_PAR_CLEAR(KIDIA:KFDIA) = 0.0_JPRB
686!AI ATTENTION
687!DO JBAND = 1,NWEIGHT_PAR
688!  PFLUX_PAR(KIDIA:KFDIA) = PFLUX_PAR(KIDIA:KFDIA) + WEIGHT_PAR(JBAND) &
689!       &  * flux%sw_dn_surf_band(IBAND_PAR(JBAND),KIDIA:KFDIA)
690!  PFLUX_PAR_CLEAR(KIDIA:KFDIA) = PFLUX_PAR_CLEAR(KIDIA:KFDIA) &
691!       &  + WEIGHT_PAR(JBAND) &
692!       &  * flux%sw_dn_surf_clear_band(IBAND_PAR(JBAND),KIDIA:KFDIA)
693!ENDDO
694
695! Compute effective broadband emissivity
696ZBLACK_BODY_NET_LW = flux%lw_dn(KIDIA:KFDIA,KLEV+1) &
697     &  - RSIGMA*PTEMPERATURE_SKIN(KIDIA:KFDIA)**4
698PEMIS_OUT(KIDIA:KFDIA) = PEMIS(KIDIA:KFDIA)
699WHERE (ABS(ZBLACK_BODY_NET_LW) > 1.0E-5)
700  PEMIS_OUT(KIDIA:KFDIA) = PFLUX_LW(KIDIA:KFDIA,KLEV+1) / ZBLACK_BODY_NET_LW
701END WHERE
702
703! Copy longwave derivatives
704!IF (YRERAD%LAPPROXLWUPDATE) THEN
705IF (rad_config%do_lw_derivatives) THEN
706  PLWDERIVATIVE(KIDIA:KFDIA,:) = flux%lw_derivatives(KIDIA:KFDIA,:)
707END IF
708
709! Store the shortwave downwelling fluxes in each albedo band
710!IF (YRERAD%LAPPROXSWUPDATE) THEN
711IF (rad_config%do_surface_sw_spectral_flux) THEN
712  PSWDIFFUSEBAND(KIDIA:KFDIA,:) = 0.0_JPRB
713  PSWDIRECTBAND (KIDIA:KFDIA,:) = 0.0_JPRB
714  DO JBAND = 1,rad_config%n_bands_sw
715    JB_ALBEDO = rad_config%i_albedo_from_band_sw(JBAND)
716    DO JLON = KIDIA,KFDIA
717      PSWDIFFUSEBAND(JLON,JB_ALBEDO) = PSWDIFFUSEBAND(JLON,JB_ALBEDO) &
718           &  + flux%sw_dn_surf_band(JBAND,JLON) &
719           &  - flux%sw_dn_direct_surf_band(JBAND,JLON)
720      PSWDIRECTBAND(JLON,JB_ALBEDO)  = PSWDIRECTBAND(JLON,JB_ALBEDO) &
721           &  + flux%sw_dn_direct_surf_band(JBAND,JLON)
722    ENDDO
723  ENDDO
724ENDIF
725
726CALL single_level%deallocate
727CALL thermodynamics%deallocate
728CALL gas%deallocate
729CALL cloud%deallocate
730CALL aerosol%deallocate
731CALL flux%deallocate
732
733IF (LHOOK) CALL DR_HOOK('RADIATION_SCHEME',1,ZHOOK_HANDLE)
734
735END SUBROUTINE RADIATION_SCHEME
Note: See TracBrowser for help on using the repository browser.