Context Navigation

radiation_scheme.F90 @ 5435

Last change on this file since 5435 was 4677, checked in by idelkadi, 16 months 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
15	SUBROUTINE 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
82	USE PARKIND1 , ONLY : JPIM, JPRB
83	USE YOMHOOK , ONLY : LHOOK, DR_HOOK
84	USE RADIATION_SETUP
85	USE YOMCST , ONLY : RSIGMA ! Stefan-Boltzmann constant
86	USE 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
95	USE radiation_single_level, ONLY : single_level_type
96	USE radiation_thermodynamics, ONLY : thermodynamics_type
97	USE radiation_gas
98	USE radiation_cloud, ONLY : cloud_type
99	USE radiation_aerosol, ONLY : aerosol_type
100	USE radiation_flux, ONLY : flux_type
101	USE radiation_interface, ONLY : radiation, set_gas_units
102	USE radiation_save, ONLY : save_inputs
103
104	USE mod_phys_lmdz_para
105
106	IMPLICIT NONE
107
108	! INPUT ARGUMENTS
109	! *** Array dimensions and ranges
110	INTEGER(KIND=JPIM),INTENT(IN) :: KIDIA ! Start column to process
111	INTEGER(KIND=JPIM),INTENT(IN) :: KFDIA ! End column to process
112	!INTEGER, INTENT(IN) :: KIDIA, KFDIA
113	INTEGER(KIND=JPIM),INTENT(IN) :: KLON ! Number of columns
114	INTEGER(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
117	INTEGER(KIND=JPIM),INTENT(IN) :: KAEROSOL
118	INTEGER(KIND=JPIM),INTENT(IN) :: NSW ! Numbe of bands
119
120	! AI ATTENTION
121	!INTEGER, PARAMETER :: KAEROSOL = 12
122
123	! *** Single-level fields
124	REAL(KIND=JPRB), INTENT(IN) :: PSOLAR_IRRADIANCE ! (W m-2)
125	REAL(KIND=JPRB), INTENT(IN) :: PMU0(KLON) ! Cosine of solar zenith ang
126	REAL(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)
130	REAL(KIND=JPRB), INTENT(IN) :: PALBEDO_DIF(KLON,NSW)
131	REAL(KIND=JPRB), INTENT(IN) :: PALBEDO_DIR(KLON,NSW)
132	! Longwave emissivity outside and inside the window region
133	REAL(KIND=JPRB), INTENT(IN) :: PEMIS(KLON)
134	REAL(KIND=JPRB), INTENT(IN) :: PEMIS_WINDOW(KLON)
135	! Longitude (radians), sine of latitude
136	REAL(KIND=JPRB), INTENT(IN) :: PGELAM(KLON)
137	REAL(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
142	REAL(KIND=JPRB), INTENT(IN) :: PPRESSURE_H(KLON,KLEV+1) ! (Pa)
143	REAL(KIND=JPRB), INTENT(IN) :: PTEMPERATURE_H(KLON,KLEV+1) ! (K)
144
145	! *** Gas mass mixing ratios on full levels
146	REAL(KIND=JPRB), INTENT(IN) :: PQ(KLON,KLEV)
147	! AI
148	REAL(KIND=JPRB), INTENT(IN) :: PQSAT(KLON,KLEV)
149	REAL(KIND=JPRB), INTENT(IN) :: PCO2
150	REAL(KIND=JPRB), INTENT(IN) :: PCH4
151	REAL(KIND=JPRB), INTENT(IN) :: PN2O
152	REAL(KIND=JPRB), INTENT(IN) :: PNO2
153	REAL(KIND=JPRB), INTENT(IN) :: PCFC11
154	REAL(KIND=JPRB), INTENT(IN) :: PCFC12
155	REAL(KIND=JPRB), INTENT(IN) :: PHCFC22
156	REAL(KIND=JPRB), INTENT(IN) :: PCCL4
157	REAL(KIND=JPRB), INTENT(IN) :: PO3(KLON,KLEV) ! AI (kg/kg) ATTENTION (Pa*kg/kg)
158	REAL(KIND=JPRB), INTENT(IN) :: PO2
159
160	! *** Cloud fraction and hydrometeor mass mixing ratios
161	REAL(KIND=JPRB), INTENT(IN) :: PCLOUD_FRAC(KLON,KLEV)
162	REAL(KIND=JPRB), INTENT(IN) :: PQ_LIQUID(KLON,KLEV)
163	REAL(KIND=JPRB), INTENT(IN) :: PQ_ICE(KLON,KLEV)
164	!REAL(KIND=JPRB), INTENT(IN) :: PQ_RAIN(KLON,KLEV)
165	REAL(KIND=JPRB), INTENT(IN) :: PQ_SNOW(KLON,KLEV)
166
167	! *** Aerosol mass mixing ratios
168	REAL(KIND=JPRB), INTENT(IN) :: PAEROSOL_OLD(KLON,6,KLEV)
169	REAL(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
175	INTEGER(KIND=JPIM), INTENT(IN) :: IDAY
176	REAL(KIND=JPRB), INTENT(IN) :: TIME
177
178
179	! OUTPUT ARGUMENTS
180
181	! *** Net fluxes on half-levels (W m-2)
182	REAL(KIND=JPRB), INTENT(OUT) :: PFLUX_SW(KLON,KLEV+1)
183	REAL(KIND=JPRB), INTENT(OUT) :: PFLUX_LW(KLON,KLEV+1)
184	REAL(KIND=JPRB), INTENT(OUT) :: PFLUX_SW_CLEAR(KLON,KLEV+1)
185	REAL(KIND=JPRB), INTENT(OUT) :: PFLUX_LW_CLEAR(KLON,KLEV+1)
186
187	!*** DN and UP flux on half-levels (W m-2)
188	REAL(KIND=JPRB), INTENT(OUT) :: PFLUX_SW_DN(KLON,KLEV+1)
189	REAL(KIND=JPRB), INTENT(OUT) :: PFLUX_LW_DN(KLON,KLEV+1)
190	REAL(KIND=JPRB), INTENT(OUT) :: PFLUX_SW_DN_CLEAR(KLON,KLEV+1)
191	REAL(KIND=JPRB), INTENT(OUT) :: PFLUX_LW_DN_CLEAR(KLON,KLEV+1)
192	REAL(KIND=JPRB), INTENT(OUT) :: PFLUX_SW_UP(KLON,KLEV+1)
193	REAL(KIND=JPRB), INTENT(OUT) :: PFLUX_LW_UP(KLON,KLEV+1)
194	REAL(KIND=JPRB), INTENT(OUT) :: PFLUX_SW_UP_CLEAR(KLON,KLEV+1)
195	REAL(KIND=JPRB), INTENT(OUT) :: PFLUX_LW_UP_CLEAR(KLON,KLEV+1)
196
197	! Direct component of surface flux into horizontal plane
198	REAL(KIND=JPRB), INTENT(OUT) :: PFLUX_DIR(KLON)
199	REAL(KIND=JPRB), INTENT(OUT) :: PFLUX_DIR_CLEAR(KLON)
200	! As PFLUX_DIR but into a plane perpendicular to the sun
201	REAL(KIND=JPRB), INTENT(OUT) :: PFLUX_DIR_INTO_SUN(KLON)
202
203	! *** Ultraviolet and photosynthetically active radiation (W m-2)
204	REAL(KIND=JPRB), INTENT(OUT) :: PFLUX_UV(KLON)
205	REAL(KIND=JPRB), INTENT(OUT) :: PFLUX_PAR(KLON)
206	REAL(KIND=JPRB), INTENT(OUT) :: PFLUX_PAR_CLEAR(KLON)
207
208	! Diagnosed longwave surface emissivity across the whole spectrum
209	REAL(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.
217	REAL(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
222	REAL(KIND=JPRB), INTENT(OUT) :: PSWDIFFUSEBAND(KLON,NSW)
223	REAL(KIND=JPRB), INTENT(OUT) :: PSWDIRECTBAND (KLON,NSW)
224
225	! LOCAL VARIABLES
226	! AI ATTENTION
227	type(config_type),save :: rad_config
228	!!$OMP THREADPRIVATE(rad_config)
229	type(driver_config_type),save :: driver_config
230	!!$OMP THREADPRIVATE(driver_config)
231	!type(config_type) :: rad_config
232	!type(driver_config_type) :: driver_config
233	TYPE(single_level_type) :: single_level
234	TYPE(thermodynamics_type) :: thermodynamics
235	TYPE(gas_type) :: gas
236	TYPE(cloud_type) :: cloud
237	TYPE(aerosol_type) :: aerosol
238	TYPE(flux_type) :: flux
239
240	! Mass mixing ratio of ozone (kg/kg)
241	REAL(KIND=JPRB) :: ZO3(KLON,KLEV)
242
243	! Cloud effective radii in microns
244	REAL(KIND=JPRB) :: ZRE_LIQUID_UM(KLON,KLEV)
245	REAL(KIND=JPRB) :: ZRE_ICE_UM(KLON,KLEV)
246
247	! Cloud overlap decorrelation length for cloud boundaries in km
248	REAL(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
256	REAL(KIND=JPRB) :: ZBLACK_BODY_NET_LW(KIDIA:KFDIA)
257
258	! Layer mass in kg m-2
259	REAL(KIND=JPRB) :: ZLAYER_MASS(KIDIA:KFDIA,KLEV)
260
261	! Time integers
262	INTEGER :: ITIM
263
264	! Loop indices
265	INTEGER :: JLON, JLEV, JBAND, JB_ALBEDO, JAER
266
267	REAL(KIND=JPRB) :: ZHOOK_HANDLE
268
269	! AI ATTENTION traitement aerosols
270	INTEGER, PARAMETER :: NAERMACC = 1
271
272	! Name of file names specified on command line
273	character(len=512) :: namelist_file
274
275	logical :: loutput=.true.
276	logical :: lprint_input=.false.
277	logical :: lprint_config=.true.
278	logical, save :: debut_ecrad=.true.
279	!$OMP THREADPRIVATE(debut_ecrad)
280	integer, save :: itap_ecrad=1
281	logical :: ok_3Deffect
282
283	IF (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
288	if (.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
293	else
294	call SETUP_RADIATION_SCHEME(loutput,namelist_file,rad_config,driver_config)
295	endif
296	!$OMP END MASTER
297	!$OMP BARRIER
298	! Fin partie initialisation et configuration
299
300	!AI juillet 2023 : verif des param de config :
301	if (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
383	endif
384
385	! AI ATTENTION
386	! Allocate memory in radiation objects
387	! emissivite avec une seule bande
388	CALL single_level%allocate(KLON, NSW, 1, &
389	& use_sw_albedo_direct=.TRUE.)
390
391	print,'********** THERMO (allocate + input) **********************************'
392	! Set thermodynamic profiles: simply copy over the half-level
393	! pressure and temperature
394	!print*,'Appel allocate thermo'
395	CALL 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
400	thermodynamics%pressure_hl (KIDIA:KFDIA,:) = PPRESSURE_H (KIDIA:KFDIA,:)
401	thermodynamics%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
413	thermodynamics%h2o_sat_liq = PQSAT
414
415	print,'******* SINGLE LEVEL VARS ********************************'
416	!AI ATTENTION
417	! Set single-level fileds
418	single_level%solar_irradiance = PSOLAR_IRRADIANCE
419	single_level%cos_sza(KIDIA:KFDIA) = PMU0(KIDIA:KFDIA)
420	single_level%skin_temperature(KIDIA:KFDIA) = PTEMPERATURE_SKIN(KIDIA:KFDIA)
421	single_level%sw_albedo(KIDIA:KFDIA,:) = PALBEDO_DIF(KIDIA:KFDIA,:)
422	single_level%sw_albedo_direct(KIDIA:KFDIA,:)=PALBEDO_DIR(KIDIA:KFDIA,:)
423	single_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)
431	ITIM = NINT(TIME / 60.0_JPRB)
432	!current_time
433	!allocate(single_level%iseed(KIDIA:KFDIA))
434	DO 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. PGELAM60100 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)
447	ENDDO
448
449	print,'******* CLOUDS (allocate + input) *****************************************'
450	!print*,'Appel Allocate clouds'
451	CALL cloud%allocate(KLON, KLEV)
452	! Set cloud fields
453	cloud%q_liq(KIDIA:KFDIA,:) = PQ_LIQUID(KIDIA:KFDIA,:)
454	cloud%q_ice(KIDIA:KFDIA,:) = PQ_ICE(KIDIA:KFDIA,:) + PQ_SNOW(KIDIA:KFDIA,:)
455	cloud%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
460	cloud%re_liq(KIDIA:KFDIA,:) = ZRE_LIQUID_UM(KIDIA:KFDIA,:)
461	cloud%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
476	ZDECORR_LEN_KM = driver_config%overlap_decorr_length
477	DO JLON = KIDIA,KFDIA
478	CALL cloud%set_overlap_param(thermodynamics, &
479	& ZDECORR_LEN_KM(JLON), &
480	& istartcol=JLON, iendcol=JLON)
481	ENDDO
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
489	CALL cloud%create_fractional_std(KLON, KLEV, driver_config%frac_std)
490
491	if (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
518	endif
519
520	print,'***** 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
576	print,'******* GAS (allocate + input) **********************************************'
577	!print*,'Appel Allocate gas'
578	CALL 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'
591	CALL gas%put(IH2O, IMassMixingRatio, PQ)
592	CALL gas%put(IO3, IMassMixingRatio, PO3)
593	CALL gas%put_well_mixed(ICO2, IMAssMixingRatio, PCO2)
594	CALL gas%put_well_mixed(ICH4, IMassMixingRatio, PCH4)
595	CALL gas%put_well_mixed(IN2O, IMassMixingRatio, PN2O)
596	CALL gas%put_well_mixed(ICFC11, IMassMixingRatio, PCFC11)
597	CALL gas%put_well_mixed(ICFC12, IMassMixingRatio, PCFC12)
598	CALL gas%put_well_mixed(IHCFC22, IMassMixingRatio, PHCFC22)
599	CALL gas%put_well_mixed(ICCL4, IMassMixingRatio, PCCL4)
600	CALL 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
603	call set_gas_units(rad_config, gas)
604
605	print,'*********** FLUX (allocate) *********************'
606	CALL flux%allocate(rad_config, 1, KLON, KLEV)
607
608	! Call radiation scheme
609	print,'***** Appel radiation scheme ************************'
610	CALL 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
615	PFLUX_SW_DN(KIDIA:KFDIA,:) = flux%sw_dn(KIDIA:KFDIA,:)
616	PFLUX_SW_UP(KIDIA:KFDIA,:) = flux%sw_up(KIDIA:KFDIA,:)
617	PFLUX_LW_DN(KIDIA:KFDIA,:) = flux%lw_dn(KIDIA:KFDIA,:)
618	PFLUX_LW_UP(KIDIA:KFDIA,:) = flux%lw_up(KIDIA:KFDIA,:)
619	PFLUX_SW_DN_CLEAR(KIDIA:KFDIA,:) = flux%sw_dn_clear(KIDIA:KFDIA,:)
620	PFLUX_SW_UP_CLEAR(KIDIA:KFDIA,:) = flux%sw_up_clear(KIDIA:KFDIA,:)
621	PFLUX_LW_DN_CLEAR(KIDIA:KFDIA,:) = flux%lw_dn_clear(KIDIA:KFDIA,:)
622	PFLUX_LW_UP_CLEAR(KIDIA:KFDIA,:) = flux%lw_up_clear(KIDIA:KFDIA,:)
623
624	! First the net fluxes
625	PFLUX_SW(KIDIA:KFDIA,:) = flux%sw_dn(KIDIA:KFDIA,:) - flux%sw_up(KIDIA:KFDIA,:)
626	PFLUX_LW(KIDIA:KFDIA,:) = flux%lw_dn(KIDIA:KFDIA,:) - flux%lw_up(KIDIA:KFDIA,:)
627	PFLUX_SW_CLEAR(KIDIA:KFDIA,:) &
628	& = flux%sw_dn_clear(KIDIA:KFDIA,:) - flux%sw_up_clear(KIDIA:KFDIA,:)
629	PFLUX_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)
641	PFLUX_DIR(KIDIA:KFDIA) = flux%sw_dn_direct(KIDIA:KFDIA,KLEV+1)
642	PFLUX_DIR_CLEAR(KIDIA:KFDIA) = flux%sw_dn_direct_clear(KIDIA:KFDIA,KLEV+1)
643	PFLUX_DIR_INTO_SUN(KIDIA:KFDIA) = 0.0_JPRB
644	WHERE (PMU0(KIDIA:KFDIA) > EPSILON(1.0_JPRB))
645	PFLUX_DIR_INTO_SUN(KIDIA:KFDIA) = PFLUX_DIR(KIDIA:KFDIA) / PMU0(KIDIA:KFDIA)
646	END 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
656	if (0.eq.1) then
657	PFLUX_UV (KIDIA:KFDIA) = 0.0_JPRB
658	DO 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)
661	ENDDO
662
663	! Compute photosynthetically active radiation similarly
664	PFLUX_PAR (KIDIA:KFDIA) = 0.0_JPRB
665	PFLUX_PAR_CLEAR(KIDIA:KFDIA) = 0.0_JPRB
666	DO 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)
672	ENDDO
673	endif
674	! Compute effective broadband emissivity
675	ZBLACK_BODY_NET_LW = flux%lw_dn(KIDIA:KFDIA,KLEV+1) &
676	& - RSIGMAPTEMPERATURE_SKIN(KIDIA:KFDIA)*4
677	PEMIS_OUT(KIDIA:KFDIA) = PEMIS(KIDIA:KFDIA)
678	WHERE (ABS(ZBLACK_BODY_NET_LW) > 1.0E-5)
679	PEMIS_OUT(KIDIA:KFDIA) = PFLUX_LW(KIDIA:KFDIA,KLEV+1) / ZBLACK_BODY_NET_LW
680	END WHERE
681
682	! Copy longwave derivatives
683	! AI ATTENTION
684	!IF (YRERAD%LAPPROXLWUPDATE) THEN
685	IF (rad_config%do_lw_derivatives) THEN
686	PLWDERIVATIVE(KIDIA:KFDIA,:) = flux%lw_derivatives(KIDIA:KFDIA,:)
687	END IF
688
689	! Store the shortwave downwelling fluxes in each albedo band
690	!AI ATTENTION
691	!IF (YRERAD%LAPPROXSWUPDATE) THEN
692	if (0.eq.1) then
693	IF (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
706	ENDIF
707	endif
708	CALL single_level%deallocate
709	CALL thermodynamics%deallocate
710	CALL gas%deallocate
711	CALL cloud%deallocate
712	CALL aerosol%deallocate
713	CALL flux%deallocate
714
715	IF (LHOOK) CALL DR_HOOK('RADIATION_SCHEME',1,ZHOOK_HANDLE)
716
717	END SUBROUTINE RADIATION_SCHEME

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format