source: LMDZ6/trunk/libf/phylmd/ecrad.v1.5.1/radiation_scheme.F90 @ 5451

Last change on this file since 5451 was 4677, checked in by idelkadi, 2 years ago

Implementation in the LMDZ code of the double call of the ECRAD radiative transfer code to estimate the 3D radiative effect of clouds.

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