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

Last change on this file since 4660 was 4570, checked in by idelkadi, 18 months ago

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