source: LMDZ6/trunk/libf/phylmd/ecrad/lmdz/radiation_scheme_mod.F90 @ 4867

Last change on this file since 4867 was 4853, checked in by idelkadi, 20 months ago

Ecrad update in LMDZ / Implementation of Ecrad double call in LMDZ

  • version 1.6.1 (commit 210d7911380f53a788c3cad73b3cf9b4e022ef87)
  • interface routines between lmdz and ecrad grouped in a new "lmdz" directory
  • double call method redesigned so as to go through the Ecrad initialization and configuration part only once for the entire simulation
  • clean-up in the read.F routine: delete unitules arguments
  • modification of default gas model in namelist (default: ECCKD)
File size: 54.8 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! ============================================================================
14module interface_lmdz_ecrad
15
16IMPLICIT NONE       
17
18contains
19
20SUBROUTINE RADIATION_SCHEME &
21! Inputs
22     & (KIDIA, KFDIA, KLON, KLEV, KAEROSOL, NSW, &
23     &  namelist_file, ok_3Deffect, &
24     &  debut, ok_volcan, flag_aerosol_strat, &
25     &  IDAY, TIME, &
26     &  PSOLAR_IRRADIANCE, &
27     &  PMU0, PTEMPERATURE_SKIN, &
28     &  PALBEDO_DIF, PALBEDO_DIR, &
29     &  PEMIS, PEMIS_WINDOW, &
30     &  PGELAM, PGEMU, &
31     &  PPRESSURE_H, PTEMPERATURE_H, PQ, PQSAT, &
32     &  PCO2, PCH4, PN2O, PNO2, PCFC11, PCFC12, PHCFC22, &
33     &  PCCL4, PO3, PO2, &
34     &  PCLOUD_FRAC, PQ_LIQUID, PQ_ICE, PQ_SNOW, &
35     &  ZRE_LIQUID_UM, ZRE_ICE_UM, &
36     &  PAEROSOL_OLD, PAEROSOL, &
37! Outputs
38     &  PFLUX_SW, PFLUX_LW, PFLUX_SW_CLEAR, PFLUX_LW_CLEAR, &
39     &  PFLUX_SW_DN, PFLUX_LW_DN, PFLUX_SW_DN_CLEAR, PFLUX_LW_DN_CLEAR, &
40     &  PFLUX_SW_UP, PFLUX_LW_UP, PFLUX_SW_UP_CLEAR, PFLUX_LW_UP_CLEAR, &
41     &  PFLUX_DIR, PFLUX_DIR_CLEAR, PFLUX_DIR_INTO_SUN, &
42     &  PFLUX_UV, PFLUX_PAR, PFLUX_PAR_CLEAR, &
43     &  PEMIS_OUT, PLWDERIVATIVE, &
44     &  PSWDIFFUSEBAND, PSWDIRECTBAND, &
45     &  ecrad_cloud_cover_sw)
46
47! RADIATION_SCHEME - Interface to modular radiation scheme
48!
49! (C) Copyright 2015- ECMWF.
50!
51! This software is licensed under the terms of the Apache Licence Version 2.0
52! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
53!
54! In applying this licence, ECMWF does not waive the privileges and immunities
55! granted to it by virtue of its status as an intergovernmental organisation
56! nor does it submit to any jurisdiction.
57!
58! PURPOSE
59! -------
60!   The modular radiation scheme is contained in a separate
61!   library. This routine puts the the IFS arrays into appropriate
62!   objects, computing the additional data that is required, and sends
63!   it to the radiation scheme.  It returns net fluxes and surface
64!   flux components needed by the rest of the model.
65!
66!   Lower case is used for variables and types taken from the
67!   radiation library
68!
69! INTERFACE
70! ---------
71!    RADIATION_SCHEME is called from RADLSWR. The
72!    SETUP_RADIATION_SCHEME routine (in the RADIATION_SETUP module)
73!    should have been run first.
74!
75! AUTHOR
76! ------
77!   Robin Hogan, ECMWF
78!   Original: 2015-09-16
79!
80! MODIFICATIONS
81! -------------
82!
83! TO DO
84! -----
85!
86!-----------------------------------------------------------------------
87
88! Modules from ifs or ifsaux libraries
89USE PARKIND1 , ONLY : JPIM, JPRB
90USE YOMHOOK  , ONLY : LHOOK, DR_HOOK
91USE RADIATION_SETUP
92USE YOMCST   , ONLY : RSIGMA ! Stefan-Boltzmann constant
93!USE RADIATION_SETUP, ONLY : SETUP_RADIATION_SCHEME, &
94!                         &  config_type, driver_config_type, &
95!                         &  NWEIGHT_UV,  IBAND_UV,  WEIGHT_UV, &
96!                         &  NWEIGHT_PAR, IBAND_PAR, WEIGHT_PAR, &
97!                         &  ITYPE_TROP_BG_AER,  TROP_BG_AER_MASS_EXT, &
98!                         &  ITYPE_STRAT_BG_AER, STRAT_BG_AER_MASS_EXT, &
99!                         &  ISolverSpartacus
100
101! Modules from radiation library
102USE radiation_single_level,   ONLY : single_level_type
103USE radiation_thermodynamics, ONLY : thermodynamics_type
104USE radiation_gas
105USE radiation_cloud,          ONLY : cloud_type
106USE radiation_aerosol,        ONLY : aerosol_type
107USE radiation_flux,           ONLY : flux_type
108USE radiation_interface,      ONLY : radiation, set_gas_units
109USE radiation_save,           ONLY : save_inputs
110
111USE mod_phys_lmdz_para
112
113IMPLICIT NONE
114
115! INPUT ARGUMENTS
116! *** Array dimensions and ranges
117INTEGER(KIND=JPIM),INTENT(IN) :: KIDIA    ! Start column to process
118INTEGER(KIND=JPIM),INTENT(IN) :: KFDIA    ! End column to process
119!INTEGER, INTENT(IN) :: KIDIA, KFDIA
120INTEGER(KIND=JPIM),INTENT(IN) :: KLON     ! Number of columns
121INTEGER(KIND=JPIM),INTENT(IN) :: KLEV     ! Number of levels
122!INTEGER, INTENT(IN) :: KLON, KLEV
123!INTEGER(KIND=JPIM),INTENT(IN) :: KAEROLMDZ ! Number of aerosol types
124INTEGER(KIND=JPIM),INTENT(IN) :: KAEROSOL
125INTEGER(KIND=JPIM),INTENT(IN) :: NSW ! Numbe of bands
126
127! AI ATTENTION
128!INTEGER, PARAMETER :: KAEROSOL = 12
129
130! *** Single-level fields
131REAL(KIND=JPRB),   INTENT(IN) :: PSOLAR_IRRADIANCE ! (W m-2)
132REAL(KIND=JPRB),   INTENT(IN) :: PMU0(KLON) ! Cosine of solar zenith ang
133REAL(KIND=JPRB),   INTENT(IN) :: PTEMPERATURE_SKIN(KLON) ! (K)
134! Diffuse and direct components of surface shortwave albedo
135!REAL(KIND=JPRB),   INTENT(IN) :: PALBEDO_DIF(KLON,YRERAD%NSW)
136!REAL(KIND=JPRB),   INTENT(IN) :: PALBEDO_DIR(KLON,YRERAD%NSW)
137REAL(KIND=JPRB),   INTENT(IN) :: PALBEDO_DIF(KLON,NSW)
138REAL(KIND=JPRB),   INTENT(IN) :: PALBEDO_DIR(KLON,NSW)
139! Longwave emissivity outside and inside the window region
140REAL(KIND=JPRB),   INTENT(IN) :: PEMIS(KLON)
141REAL(KIND=JPRB),   INTENT(IN) :: PEMIS_WINDOW(KLON)
142! Longitude (radians), sine of latitude
143REAL(KIND=JPRB),   INTENT(IN) :: PGELAM(KLON)
144REAL(KIND=JPRB),   INTENT(IN) :: PGEMU(KLON)
145! Land-sea mask
146!REAL(KIND=JPRB),   INTENT(IN) :: PLAND_SEA_MASK(KLON)
147
148! *** Variables on half levels
149REAL(KIND=JPRB),   INTENT(IN) :: PPRESSURE_H(KLON,KLEV+1)    ! (Pa)
150REAL(KIND=JPRB),   INTENT(IN) :: PTEMPERATURE_H(KLON,KLEV+1) ! (K)
151
152! *** Gas mass mixing ratios on full levels
153REAL(KIND=JPRB),   INTENT(IN) :: PQ(KLON,KLEV)
154! AI
155REAL(KIND=JPRB),   INTENT(IN) :: PQSAT(KLON,KLEV)
156REAL(KIND=JPRB),   INTENT(IN) :: PCO2
157REAL(KIND=JPRB),   INTENT(IN) :: PCH4
158REAL(KIND=JPRB),   INTENT(IN) :: PN2O
159REAL(KIND=JPRB),   INTENT(IN) :: PNO2
160REAL(KIND=JPRB),   INTENT(IN) :: PCFC11
161REAL(KIND=JPRB),   INTENT(IN) :: PCFC12
162REAL(KIND=JPRB),   INTENT(IN) :: PHCFC22
163REAL(KIND=JPRB),   INTENT(IN) :: PCCL4
164REAL(KIND=JPRB),   INTENT(IN) :: PO3(KLON,KLEV) ! AI (kg/kg) ATTENTION (Pa*kg/kg)
165REAL(KIND=JPRB),   INTENT(IN) :: PO2
166
167! *** Cloud fraction and hydrometeor mass mixing ratios
168REAL(KIND=JPRB),   INTENT(IN) :: PCLOUD_FRAC(KLON,KLEV)
169REAL(KIND=JPRB),   INTENT(IN) :: PQ_LIQUID(KLON,KLEV)
170REAL(KIND=JPRB),   INTENT(IN) :: PQ_ICE(KLON,KLEV)
171!REAL(KIND=JPRB),   INTENT(IN) :: PQ_RAIN(KLON,KLEV)
172REAL(KIND=JPRB),   INTENT(IN) :: PQ_SNOW(KLON,KLEV)
173
174! *** Aerosol mass mixing ratios
175REAL(KIND=JPRB),   INTENT(IN) :: PAEROSOL_OLD(KLON,6,KLEV)
176REAL(KIND=JPRB),   INTENT(IN) :: PAEROSOL(KLON,KLEV,KAEROSOL)
177
178!REAL(KIND=JPRB),   INTENT(IN) :: PCCN_LAND(KLON)
179!REAL(KIND=JPRB),   INTENT(IN) :: PCCN_SEA(KLON)
180
181!AI mars 2021
182INTEGER(KIND=JPIM), INTENT(IN) :: IDAY
183REAL(KIND=JPRB), INTENT(IN)    :: TIME
184
185! Name of file names specified on command line
186character(len=512), INTENT(IN) :: namelist_file
187logical, INTENT(IN)            :: ok_3Deffect, debut, ok_volcan
188INTEGER(KIND=JPIM), INTENT(IN) :: flag_aerosol_strat
189
190
191! OUTPUT ARGUMENTS
192
193! *** Net fluxes on half-levels (W m-2)
194REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW(KLON,KLEV+1)
195REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW(KLON,KLEV+1)
196REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_CLEAR(KLON,KLEV+1)
197REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_CLEAR(KLON,KLEV+1)
198
199!*** DN and UP flux on half-levels (W m-2)
200REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_DN(KLON,KLEV+1)
201REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_DN(KLON,KLEV+1)
202REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_DN_CLEAR(KLON,KLEV+1)
203REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_DN_CLEAR(KLON,KLEV+1)
204REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_UP(KLON,KLEV+1)
205REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_UP(KLON,KLEV+1)
206REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_UP_CLEAR(KLON,KLEV+1)
207REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_UP_CLEAR(KLON,KLEV+1)
208
209! Direct component of surface flux into horizontal plane
210REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_DIR(KLON)
211REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_DIR_CLEAR(KLON)
212! As PFLUX_DIR but into a plane perpendicular to the sun
213REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_DIR_INTO_SUN(KLON)
214
215! *** Ultraviolet and photosynthetically active radiation (W m-2)
216REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_UV(KLON)
217REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_PAR(KLON)
218REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_PAR_CLEAR(KLON)
219
220! Diagnosed longwave surface emissivity across the whole spectrum
221REAL(KIND=JPRB),  INTENT(OUT) :: PEMIS_OUT(KLON)   
222
223! Partial derivative of total-sky longwave upward flux at each level
224! with respect to upward flux at surface, used to correct heating
225! rates at gridpoints/timesteps between calls to the full radiation
226! scheme.  Note that this version uses the convention of level index
227! increasing downwards, unlike the local variable ZLwDerivative that
228! is returned from the LW radiation scheme.
229REAL(KIND=JPRB),  INTENT(OUT) :: PLWDERIVATIVE(KLON,KLEV+1)
230
231! Surface diffuse and direct downwelling shortwave flux in each
232! shortwave albedo band, used in RADINTG to update the surface fluxes
233! accounting for high-resolution albedo information
234REAL(KIND=JPRB),  INTENT(OUT) :: PSWDIFFUSEBAND(KLON,NSW)
235REAL(KIND=JPRB),  INTENT(OUT) :: PSWDIRECTBAND (KLON,NSW)
236
237!AI Nov 2023
238REAL(KIND=JPRB),  INTENT(OUT) :: ecrad_cloud_cover_sw(KLON)
239
240! LOCAL VARIABLES
241! AI ATTENTION
242type(config_type),save         :: rad_config
243!!$OMP THREADPRIVATE(rad_config)
244type(driver_config_type),save  :: driver_config
245!!$OMP THREADPRIVATE(driver_config)
246!type(config_type)        :: rad_config
247!type(driver_config_type)  :: driver_config
248TYPE(single_level_type)   :: single_level
249TYPE(thermodynamics_type) :: thermodynamics
250TYPE(gas_type)            :: gas
251TYPE(cloud_type)          :: cloud
252TYPE(aerosol_type)        :: aerosol
253TYPE(flux_type)           :: flux
254
255! Mass mixing ratio of ozone (kg/kg)
256REAL(KIND=JPRB)           :: ZO3(KLON,KLEV)
257
258! Cloud effective radii in microns
259REAL(KIND=JPRB)           :: ZRE_LIQUID_UM(KLON,KLEV)
260REAL(KIND=JPRB)           :: ZRE_ICE_UM(KLON,KLEV)
261
262! Cloud overlap decorrelation length for cloud boundaries in km
263REAL(KIND=JPRB)           :: ZDECORR_LEN_KM(KLON)
264
265! Ratio of cloud overlap decorrelation length for cloud water
266! inhomogeneities to that for cloud boundaries (typically 0.5)
267!REAL(KIND=JPRB)           :: ZDECORR_LEN_RATIO = 0.5_jprb
268
269! The surface net longwave flux if the surface was a black body, used
270! to compute the effective broadband surface emissivity
271REAL(KIND=JPRB)           :: ZBLACK_BODY_NET_LW(KIDIA:KFDIA)
272
273! Layer mass in kg m-2
274REAL(KIND=JPRB)           :: ZLAYER_MASS(KIDIA:KFDIA,KLEV)
275
276! Time integers
277INTEGER :: ITIM
278
279! Loop indices
280INTEGER :: JLON, JLEV, JBAND, JB_ALBEDO, JAER
281
282REAL(KIND=JPRB) :: ZHOOK_HANDLE
283
284! AI ATTENTION traitement aerosols
285INTEGER, PARAMETER :: NAERMACC = 1
286
287logical :: loutput=.true.
288logical :: lprint_input=.false.
289logical :: lprint_config=.false.
290logical, save :: debut_ecrad=.true.
291!$OMP THREADPRIVATE(debut_ecrad)
292integer, save :: itap_ecrad=0
293!$OMP THREADPRIVATE(itap_ecrad)
294
295REAL(KIND=JPRB) ::  inv_cloud_effective_size(KLON,KLEV)
296REAL(KIND=JPRB) ::  inv_inhom_effective_size(KLON,KLEV)
297
298integer :: irang
299
300
301IF (LHOOK) CALL DR_HOOK('RADIATION_SCHEME',0,ZHOOK_HANDLE)
302
303  print*,'Entree radiation_scheme, ok_3Deffect, namelist_file = ', &
304          ok_3Deffect, namelist_file
305! A.I juillet 2023 :
306! Initialisation dans radiation_setup au 1er passage dans Ecrad
307!$OMP MASTER
308!if (.not.ok_3Deffect) then
309  if (debut_ecrad) then
310   call SETUP_RADIATION_SCHEME(loutput,namelist_file,rad_config,driver_config)
311   debut_ecrad=.false.
312  endif
313!else
314!   call SETUP_RADIATION_SCHEME(loutput,namelist_file,rad_config,driver_config)
315!endif
316!$OMP END MASTER
317!$OMP BARRIER
318! Fin partie initialisation et configuration
319
320!AI print fichiers namelist utilise
321!if (is_omp_root) then
322! itap_ecrad=itap_ecrad+1
323! print*,'Dans radiation_scheme itap_ecrad, mpi_rank, omp_rank, namelist_file : ', &
324!         itap_ecrad, mpi_rank, omp_rank, namelist_file
325!else
326! print*,'mpi_rank omp_rank, namelist_file :', mpi_rank, omp_rank, namelist_file
327!endif
328
329! AI 11 23 Allocates depplaces au debut
330print*,'*********** ALLOCATES *******************************'
331! AI ATTENTION
332! Allocate memory in radiation objects
333! emissivite avec une seule bande
334CALL single_level%allocate(KLON, NSW, 1, &
335     &                     use_sw_albedo_direct=.TRUE.)
336CALL thermodynamics%allocate(KLON, KLEV, use_h2o_sat=.true.)
337CALL cloud%allocate(KLON, KLEV)
338CALL aerosol%allocate(KLON, 1, KLEV, KAEROSOL)
339CALL gas%allocate(KLON, KLEV)
340CALL flux%allocate(rad_config, 1, KLON, KLEV)
341
342print*,'************* THERMO (input) ************************************'
343! Set thermodynamic profiles: simply copy over the half-level
344! pressure and temperature
345! AI
346! pressure_hl > paprs
347! temperature_hl calculee dans radlsw de la meme facon que pour RRTM
348thermodynamics%pressure_hl   (KIDIA:KFDIA,:) = PPRESSURE_H   (KIDIA:KFDIA,:)
349thermodynamics%temperature_hl(KIDIA:KFDIA,:) = PTEMPERATURE_H(KIDIA:KFDIA,:)
350!print*,'Compute saturation specific humidity'
351! Compute saturation specific humidity, used to hydrate aerosols. The
352! "2" for the last argument indicates that the routine is not being
353! called from within the convection scheme.
354!CALL SATUR(KIDIA, KFDIA, KLON, 1, KLEV, &
355!     &  PPRESSURE, PTEMPERATURE, thermodynamics%h2o_sat_liq, 2)
356! Alternative approximate version using temperature and pressure from
357! the thermodynamics structure
358!CALL thermodynamics%calc_saturation_wrt_liquid(KIDIA, KFDIA)
359!AI ATTENTION
360thermodynamics%h2o_sat_liq = PQSAT
361
362print*,'********** SINGLE LEVEL VARS **********************************'
363!AI ATTENTION
364! Set single-level fileds
365single_level%solar_irradiance              = PSOLAR_IRRADIANCE
366single_level%cos_sza(KIDIA:KFDIA)          = PMU0(KIDIA:KFDIA)
367single_level%skin_temperature(KIDIA:KFDIA) = PTEMPERATURE_SKIN(KIDIA:KFDIA)
368single_level%sw_albedo(KIDIA:KFDIA,:)      = PALBEDO_DIF(KIDIA:KFDIA,:)
369single_level%sw_albedo_direct(KIDIA:KFDIA,:)=PALBEDO_DIR(KIDIA:KFDIA,:)
370single_level%lw_emissivity(KIDIA:KFDIA,1)  = PEMIS(KIDIA:KFDIA)
371!single_level%lw_emissivity(KIDIA:KFDIA,2)  = PEMIS_WINDOW(KIDIA:KFDIA)
372
373! Create the relevant seed from date and time get the starting day
374! and number of minutes since start
375!IDAY = NDD(NINDAT)
376!cur_day
377!ITIM = NINT(NSTEP * YRRIP%TSTEP / 60.0_JPRB)
378!ITIM = NINT(TIME / 60.0_JPRB)
379!current_time
380!allocate(single_level%iseed(KIDIA:KFDIA))
381!DO JLON = KIDIA, KFDIA
382  ! This method gives a unique value for roughly every 1-km square
383  ! on the globe and every minute.  ASIN(PGEMU)*60 gives rough
384  ! latitude in degrees, which we multiply by 100 to give a unique
385  ! value for roughly every km. PGELAM*60*100 gives a unique number
386  ! for roughly every km of longitude around the equator, which we
387  ! multiply by 180*100 so there is no overlap with the latitude
388  ! values.  The result can be contained in a 32-byte integer (but
389  ! since random numbers are generated with the help of integer
390  ! overflow, it should not matter if the number did overflow).
391!  single_level%iseed(JLON) = ITIM + IDAY &
392!       &  +  NINT(PGELAM(JLON)*108000000.0_JPRB &
393!       &          + ASIN(PGEMU(JLON))*6000.0_JPRB)
394!ENDDO
395!AI Nov 23
396! Simple initialization of the seeds for the Monte Carlo scheme
397call single_level%init_seed_simple(kidia, kfdia)
398
399print*,'********** CLOUDS (allocate + input) *******************************************'
400!print*,'Appel Allocate clouds'
401! Set cloud fields
402cloud%q_liq(KIDIA:KFDIA,:)    = PQ_LIQUID(KIDIA:KFDIA,:)
403cloud%q_ice(KIDIA:KFDIA,:)    = PQ_ICE(KIDIA:KFDIA,:) + PQ_SNOW(KIDIA:KFDIA,:)
404cloud%fraction(KIDIA:KFDIA,:) = PCLOUD_FRAC(KIDIA:KFDIA,:)
405!!! ok AI ATTENTION a voir avec JL
406! Compute effective radi and convert to metres
407! AI. : on passe directement les champs de LMDZ
408cloud%re_liq(KIDIA:KFDIA,:) = ZRE_LIQUID_UM(KIDIA:KFDIA,:)
409cloud%re_ice(KIDIA:KFDIA,:) = ZRE_ICE_UM(KIDIA:KFDIA,:)
410! Get the cloud overlap decorrelation length (for cloud boundaries),
411! in km, according to the parameterization specified by NDECOLAT,
412! and insert into the "cloud" object. Also get the ratio of
413! decorrelation lengths for cloud water content inhomogeneities and
414! cloud boundaries, and set it in the "rad_config" object.
415! IFS :
416!CALL CLOUD_OVERLAP_DECORR_LEN(KIDIA, KFDIA, KLON, PGEMU, YRERAD%NDECOLAT, &
417!     &    ZDECORR_LEN_KM, PDECORR_LEN_RATIO=ZDECORR_LEN_RATIO)
418! AI valeur dans namelist
419! rad_config%cloud_inhom_decorr_scaling = ZDECORR_LEN_RATIO
420!AI ATTENTION meme valeur que dans offline
421! A mettre dans namelist
422ZDECORR_LEN_KM = driver_config%overlap_decorr_length
423DO JLON = KIDIA,KFDIA
424  CALL cloud%set_overlap_param(thermodynamics, &
425       &                 ZDECORR_LEN_KM(JLON), &
426       &                 istartcol=JLON, iendcol=JLON)
427ENDDO
428! IFS :
429! Cloud water content fractional standard deviation is configurable
430! from namelist NAERAD but must be globally constant. Before it was
431! hard coded at 1.0.
432!CALL cloud%create_fractional_std(KLON, KLEV, YRERAD%RCLOUD_FRAC_STD)
433! AI ATTENTION frac_std=0.75 meme valeur que dans la version offline
434CALL cloud%create_fractional_std(KLON, KLEV, driver_config%frac_std)
435
436!if (ok_3Deffect) then
437! if (driver_config%ok_effective_size) then
438!     call cloud%create_inv_cloud_effective_size_eta(klon, klev, &
439!               &  thermodynamics%pressure_hl, &
440!               &  driver_config%low_inv_effective_size, &
441!               &  driver_config%middle_inv_effective_size, &
442!               &  driver_config%high_inv_effective_size, 0.8_jprb, 0.45_jprb, &
443!               &  KIDIA, KFDIA)
444!   else if (driver_config%ok_separation) then
445!     call cloud%param_cloud_effective_separation_eta(klon, klev, &
446!               &  thermodynamics%pressure_hl, &
447!               &  driver_config%cloud_separation_scale_surface, &
448!               &  driver_config%cloud_separation_scale_toa, &
449!               &  driver_config%cloud_separation_scale_power, &
450!               &  driver_config%cloud_inhom_separation_factor, &
451!               &  KIDIA, KFDIA)
452!   endif     
453! else 
454 if (rad_config%i_solver_sw == ISolverSPARTACUS &
455      & .or.   rad_config%i_solver_lw == ISolverSPARTACUS) then
456   ! AI ! Read cloud properties needed by SPARTACUS
457   if (driver_config%ok_effective_size) then
458     call cloud%create_inv_cloud_effective_size_eta(klon, klev, &
459               &  thermodynamics%pressure_hl, &
460               &  driver_config%low_inv_effective_size, &
461               &  driver_config%middle_inv_effective_size, &
462               &  driver_config%high_inv_effective_size, 0.8_jprb, 0.45_jprb, &
463               &  KIDIA, KFDIA)
464   else if (driver_config%ok_separation) then
465     call cloud%param_cloud_effective_separation_eta(klon, klev, &
466               &  thermodynamics%pressure_hl, &
467               &  driver_config%cloud_separation_scale_surface, &
468               &  driver_config%cloud_separation_scale_toa, &
469               &  driver_config%cloud_separation_scale_power, &
470               &  driver_config%cloud_inhom_separation_factor, &
471               &  KIDIA, KFDIA)
472   endif
473 endif 
474!endif
475
476print*,'******** AEROSOLS (input) **************************************'
477!IF (NAERMACC > 0) THEN
478!ELSE
479!  CALL aerosol%allocate(KLON, 1, KLEV, 6) ! Tegen climatology
480!ENDIF
481! Compute the dry mass of each layer neglecting humidity effects, in
482! kg m-2, needed to scale some of the aerosol inputs
483! AI commente ATTENTION
484!CALL thermodynamics%get_layer_mass(ZLAYER_MASS)
485
486! Copy over aerosol mass mixing ratio
487!IF (NAERMACC > 0) THEN
488
489  ! MACC aerosol climatology - this is already in mass mixing ratio
490  ! units with the required array orientation so we can copy it over
491  ! directly
492  aerosol%mixing_ratio(KIDIA:KFDIA,:,:) = PAEROSOL(KIDIA:KFDIA,:,:)
493
494  ! Add the tropospheric and stratospheric backgrounds contained in the
495  ! old Tegen arrays - this is very ugly!
496! AI ATTENTION
497!  IF (TROP_BG_AER_MASS_EXT > 0.0_JPRB) THEN
498!    aerosol%mixing_ratio(KIDIA:KFDIA,:,ITYPE_TROP_BG_AER) &
499!         &  = aerosol%mixing_ratio(KIDIA:KFDIA,:,ITYPE_TROP_BG_AER) &
500!         &  + PAEROSOL_OLD(KIDIA:KFDIA,1,:) &
501!         &  / (ZLAYER_MASS * TROP_BG_AER_MASS_EXT)
502!  ENDIF
503!  IF (STRAT_BG_AER_MASS_EXT > 0.0_JPRB) THEN
504!    aerosol%mixing_ratio(KIDIA:KFDIA,:,ITYPE_STRAT_BG_AER) &
505!         &  = aerosol%mixing_ratio(KIDIA:KFDIA,:,ITYPE_STRAT_BG_AER) &
506!         &  + PAEROSOL_OLD(KIDIA:KFDIA,6,:) &
507!         &  / (ZLAYER_MASS * STRAT_BG_AER_MASS_EXT)
508!  ENDIF
509
510!ELSE
511
512  ! Tegen aerosol climatology - the array PAEROSOL_OLD contains the
513  ! 550-nm optical depth in each layer. The optics data file
514  ! aerosol_ifs_rrtm_tegen.nc does not contain mass extinction
515  ! coefficient, but a scaling factor that the 550-nm optical depth
516  ! should be multiplied by to obtain the optical depth in each
517  ! spectral band.  Therefore, in order for the units to work out, we
518  ! need to divide by the layer mass (in kg m-2) to obtain the 550-nm
519  ! cross-section per unit mass of dry air (so in m2 kg-1).  We also
520  ! need to permute the array.
521!  DO JLEV = 1,KLEV
522!    DO JAER = 1,6
523!      aerosol%mixing_ratio(KIDIA:KFDIA,JLEV,JAER) &
524!         &  = PAEROSOL_OLD(KIDIA:KFDIA,JAER,JLEV) &
525!         &  / ZLAYER_MASS(KIDIA:KFDIA,JLEV)
526!    ENDDO
527!  ENDDO
528!ENDIF
529
530print*,'********** GAS (input) ************************************************'
531!print*,'Appel Allocate gas'
532! Convert ozone Pa*kg/kg to kg/kg
533! AI ATTENTION
534!DO JLEV = 1,KLEV
535!  DO JLON = KIDIA,KFDIA
536!    ZO3(JLON,JLEV) = PO3_DP(JLON,JLEV) &
537!         &         / (PPRESSURE_H(JLON,JLEV+1)-PPRESSURE_H(JLON,JLEV))
538!  ENDDO
539!ENDDO
540!  Insert gas mixing ratios
541!print*,'Insert gas mixing ratios'
542CALL gas%put(IH2O,    IMassMixingRatio, PQ)
543CALL gas%put(IO3,     IMassMixingRatio, PO3)
544CALL gas%put_well_mixed(ICO2,    IMAssMixingRatio, PCO2)
545CALL gas%put_well_mixed(ICH4,    IMassMixingRatio, PCH4)
546CALL gas%put_well_mixed(IN2O,    IMassMixingRatio, PN2O)
547CALL gas%put_well_mixed(ICFC11,  IMassMixingRatio, PCFC11)
548CALL gas%put_well_mixed(ICFC12,  IMassMixingRatio, PCFC12)
549CALL gas%put_well_mixed(IHCFC22, IMassMixingRatio, PHCFC22)
550CALL gas%put_well_mixed(ICCL4,   IMassMixingRatio, PCCL4)
551CALL gas%put_well_mixed(IO2,     IMassMixingRatio, PO2)
552! Ensure the units of the gas mixing ratios are what is required by
553! the gas absorption model
554call set_gas_units(rad_config, gas)
555
556! Call radiation scheme
557!print*,'*** Appel radiation *** namelist **** omp_rank ****', &
558!         omp_rank, namelist_file
559!   if (rad_config%i_solver_sw == ISolverSPARTACUS) then
560!    if (driver_config%ok_separation) then
561!      print*,'Avant radiation, mpi_rank, omp_rank, size, chape inv_cloud = ',&
562!              mpi_rank, omp_rank, &
563!              shape(cloud%inv_cloud_effective_size), &
564!              size(cloud%inv_cloud_effective_size)     
565!      do jlon=KIDIA, KFDIA     
566!        do jlev=1,klev
567!         print*,' Avant radiation mpi_rank, omp_rank, jlon, jlev, &
568!            &     cloud%inv_cloud_effective_size =', mpi_rank, &
569!            &     omp_rank, jlon, jlev, &
570!            &     cloud%inv_cloud_effective_size(jlon,jlev)
571!        enddo
572!      enddo
573!       cloud%inv_cloud_effective_size=inv_cloud_effective_size
574!       cloud%inv_inhom_effective_size=inv_inhom_effective_size
575!    endif
576!   endif
577CALL radiation(KLON, KLEV, KIDIA, KFDIA, rad_config, &
578     &  single_level, thermodynamics, gas, cloud, aerosol, flux)
579
580 if (rad_config%use_aerosols) then
581    if (rad_config%i_gas_model_sw == IGasModelIFSRRTMG .or. &
582     &  rad_config%i_gas_model_lw == IGasModelIFSRRTMG) then     
583       CALL aeropt_5wv_ecrad(kidia, kfdia, 1, klev,     &
584     &                       rad_config,thermodynamics, aerosol)
585    endif                 
586
587    if (flag_aerosol_strat.eq.2) then
588       CALL readaerosolstrato_ecrad(rad_config, debut, ok_volcan)
589    endif
590 endif               
591
592print*,'*********** Sortie flux ****************'
593! Cloud cover
594ecrad_cloud_cover_sw = flux%cloud_cover_sw
595! Compute required output fluxes
596! DN and UP flux
597PFLUX_SW_DN(KIDIA:KFDIA,:) = flux%sw_dn(KIDIA:KFDIA,:)
598PFLUX_SW_UP(KIDIA:KFDIA,:) = flux%sw_up(KIDIA:KFDIA,:)
599PFLUX_LW_DN(KIDIA:KFDIA,:) = flux%lw_dn(KIDIA:KFDIA,:)
600PFLUX_LW_UP(KIDIA:KFDIA,:) = flux%lw_up(KIDIA:KFDIA,:)
601PFLUX_SW_DN_CLEAR(KIDIA:KFDIA,:) = flux%sw_dn_clear(KIDIA:KFDIA,:)
602PFLUX_SW_UP_CLEAR(KIDIA:KFDIA,:) = flux%sw_up_clear(KIDIA:KFDIA,:)
603PFLUX_LW_DN_CLEAR(KIDIA:KFDIA,:) = flux%lw_dn_clear(KIDIA:KFDIA,:)
604PFLUX_LW_UP_CLEAR(KIDIA:KFDIA,:) = flux%lw_up_clear(KIDIA:KFDIA,:)
605! First the net fluxes
606PFLUX_SW(KIDIA:KFDIA,:) = flux%sw_dn(KIDIA:KFDIA,:) - flux%sw_up(KIDIA:KFDIA,:)
607PFLUX_LW(KIDIA:KFDIA,:) = flux%lw_dn(KIDIA:KFDIA,:) - flux%lw_up(KIDIA:KFDIA,:)
608PFLUX_SW_CLEAR(KIDIA:KFDIA,:) &
609     &  = flux%sw_dn_clear(KIDIA:KFDIA,:) - flux%sw_up_clear(KIDIA:KFDIA,:)
610PFLUX_LW_CLEAR(KIDIA:KFDIA,:) &
611     &  = flux%lw_dn_clear(KIDIA:KFDIA,:) - flux%lw_up_clear(KIDIA:KFDIA,:)
612! Now the surface fluxes
613!PFLUX_SW_DN_SURF(KIDIA:KFDIA) = flux%sw_dn(KIDIA:KFDIA,KLEV+1)
614!PFLUX_LW_DN_SURF(KIDIA:KFDIA) = flux%lw_dn(KIDIA:KFDIA,KLEV+1)
615!PFLUX_SW_UP_SURF(KIDIA:KFDIA) = flux%sw_up(KIDIA:KFDIA,KLEV+1)
616!PFLUX_LW_UP_SURF(KIDIA:KFDIA) = flux%lw_up(KIDIA:KFDIA,KLEV+1)
617!PFLUX_SW_DN_CLEAR_SURF(KIDIA:KFDIA) = flux%sw_dn_clear(KIDIA:KFDIA,KLEV+1)
618!PFLUX_LW_DN_CLEAR_SURF(KIDIA:KFDIA) = flux%lw_dn_clear(KIDIA:KFDIA,KLEV+1)
619!PFLUX_SW_UP_CLEAR_SURF(KIDIA:KFDIA) = flux%sw_up_clear(KIDIA:KFDIA,KLEV+1)
620!PFLUX_LW_UP_CLEAR_SURF(KIDIA:KFDIA) = flux%lw_up_clear(KIDIA:KFDIA,KLEV+1)
621PFLUX_DIR(KIDIA:KFDIA) = flux%sw_dn_direct(KIDIA:KFDIA,KLEV+1)
622PFLUX_DIR_CLEAR(KIDIA:KFDIA) = flux%sw_dn_direct_clear(KIDIA:KFDIA,KLEV+1)
623PFLUX_DIR_INTO_SUN(KIDIA:KFDIA) = 0.0_JPRB
624WHERE (PMU0(KIDIA:KFDIA) > EPSILON(1.0_JPRB))
625  PFLUX_DIR_INTO_SUN(KIDIA:KFDIA) = PFLUX_DIR(KIDIA:KFDIA) / PMU0(KIDIA:KFDIA)
626END WHERE
627! Top-of-atmosphere downwelling flux
628!PFLUX_SW_DN_TOA(KIDIA:KFDIA) = flux%sw_dn(KIDIA:KFDIA,1)
629!PFLUX_SW_UP_TOA(KIDIA:KFDIA) = flux%sw_up(KIDIA:KFDIA,1)
630!PFLUX_LW_DN_TOA(KIDIA:KFDIA) = flux%lw_dn(KIDIA:KFDIA,1)
631!PFLUX_LW_UP_TOA(KIDIA:KFDIA) = flux%lw_up(KIDIA:KFDIA,1)
632!AI ATTENTION
633if (0.eq.1) then
634PFLUX_UV       (KIDIA:KFDIA) = 0.0_JPRB
635DO JBAND = 1,NWEIGHT_UV
636  PFLUX_UV(KIDIA:KFDIA) = PFLUX_UV(KIDIA:KFDIA) + WEIGHT_UV(JBAND) &
637       &  * flux%sw_dn_surf_band(IBAND_UV(JBAND),KIDIA:KFDIA)
638ENDDO
639! Compute photosynthetically active radiation similarly
640PFLUX_PAR      (KIDIA:KFDIA) = 0.0_JPRB
641PFLUX_PAR_CLEAR(KIDIA:KFDIA) = 0.0_JPRB
642DO JBAND = 1,NWEIGHT_PAR
643  PFLUX_PAR(KIDIA:KFDIA) = PFLUX_PAR(KIDIA:KFDIA) + WEIGHT_PAR(JBAND) &
644       &  * flux%sw_dn_surf_band(IBAND_PAR(JBAND),KIDIA:KFDIA)
645  PFLUX_PAR_CLEAR(KIDIA:KFDIA) = PFLUX_PAR_CLEAR(KIDIA:KFDIA) &
646       &  + WEIGHT_PAR(JBAND) &
647       &  * flux%sw_dn_surf_clear_band(IBAND_PAR(JBAND),KIDIA:KFDIA)
648ENDDO
649endif
650! Compute effective broadband emissivity
651ZBLACK_BODY_NET_LW = flux%lw_dn(KIDIA:KFDIA,KLEV+1) &
652     &  - RSIGMA*PTEMPERATURE_SKIN(KIDIA:KFDIA)**4
653PEMIS_OUT(KIDIA:KFDIA) = PEMIS(KIDIA:KFDIA)
654WHERE (ABS(ZBLACK_BODY_NET_LW) > 1.0E-5)
655  PEMIS_OUT(KIDIA:KFDIA) = PFLUX_LW(KIDIA:KFDIA,KLEV+1) / ZBLACK_BODY_NET_LW
656END WHERE
657! Copy longwave derivatives
658! AI ATTENTION
659!IF (YRERAD%LAPPROXLWUPDATE) THEN
660IF (rad_config%do_lw_derivatives) THEN
661  PLWDERIVATIVE(KIDIA:KFDIA,:) = flux%lw_derivatives(KIDIA:KFDIA,:)
662END IF
663! Store the shortwave downwelling fluxes in each albedo band
664!AI ATTENTION
665!IF (YRERAD%LAPPROXSWUPDATE) THEN
666if (0.eq.1) then
667IF (rad_config%do_surface_sw_spectral_flux) THEN
668  PSWDIFFUSEBAND(KIDIA:KFDIA,:) = 0.0_JPRB
669  PSWDIRECTBAND (KIDIA:KFDIA,:) = 0.0_JPRB
670  DO JBAND = 1,rad_config%n_bands_sw
671    JB_ALBEDO = rad_config%i_albedo_from_band_sw(JBAND)
672    DO JLON = KIDIA,KFDIA
673      PSWDIFFUSEBAND(JLON,JB_ALBEDO) = PSWDIFFUSEBAND(JLON,JB_ALBEDO) &
674           &  + flux%sw_dn_surf_band(JBAND,JLON) &
675           &  - flux%sw_dn_direct_surf_band(JBAND,JLON)
676      PSWDIRECTBAND(JLON,JB_ALBEDO)  = PSWDIRECTBAND(JLON,JB_ALBEDO) &
677           &  + flux%sw_dn_direct_surf_band(JBAND,JLON)
678    ENDDO
679  ENDDO
680ENDIF
681endif
682
683print*,'********** DEALLOCATIONS ************************'
684CALL single_level%deallocate
685CALL thermodynamics%deallocate
686CALL gas%deallocate
687CALL cloud%deallocate
688CALL aerosol%deallocate
689CALL flux%deallocate
690
691IF (LHOOK) CALL DR_HOOK('RADIATION_SCHEME',1,ZHOOK_HANDLE)
692
693END SUBROUTINE RADIATION_SCHEME
694
695SUBROUTINE RADIATION_SCHEME_S2 &
696! Inputs
697     & (KIDIA, KFDIA, KLON, KLEV, KAEROSOL, NSW, &
698     &  namelist_file, ok_3Deffect, &
699     &  debut, ok_volcan, flag_aerosol_strat, &
700     &  IDAY, TIME, &
701     &  PSOLAR_IRRADIANCE, &
702     &  PMU0, PTEMPERATURE_SKIN, &
703     &  PALBEDO_DIF, PALBEDO_DIR, &
704     &  PEMIS, PEMIS_WINDOW, &
705     &  PGELAM, PGEMU, &
706     &  PPRESSURE_H, PTEMPERATURE_H, PQ, PQSAT, &
707     &  PCO2, PCH4, PN2O, PNO2, PCFC11, PCFC12, PHCFC22, &
708     &  PCCL4, PO3, PO2, &
709     &  PCLOUD_FRAC, PQ_LIQUID, PQ_ICE, PQ_SNOW, &
710     &  ZRE_LIQUID_UM, ZRE_ICE_UM, &
711     &  PAEROSOL_OLD, PAEROSOL, &
712! Outputs
713     &  PFLUX_SW, PFLUX_LW, PFLUX_SW_CLEAR, PFLUX_LW_CLEAR, &
714     &  PFLUX_SW_DN, PFLUX_LW_DN, PFLUX_SW_DN_CLEAR, PFLUX_LW_DN_CLEAR, &
715     &  PFLUX_SW_UP, PFLUX_LW_UP, PFLUX_SW_UP_CLEAR, PFLUX_LW_UP_CLEAR, &
716     &  PFLUX_DIR, PFLUX_DIR_CLEAR, PFLUX_DIR_INTO_SUN, &
717     &  PFLUX_UV, PFLUX_PAR, PFLUX_PAR_CLEAR, &
718     &  PEMIS_OUT, PLWDERIVATIVE, &
719     &  PSWDIFFUSEBAND, PSWDIRECTBAND, &
720     &  ecrad_cloud_cover_sw)
721
722! RADIATION_SCHEME - Interface to modular radiation scheme
723!
724! (C) Copyright 2015- ECMWF.
725!
726! This software is licensed under the terms of the Apache Licence Version 2.0
727! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
728!
729! In applying this licence, ECMWF does not waive the privileges and immunities
730! granted to it by virtue of its status as an intergovernmental organisation
731! nor does it submit to any jurisdiction.
732!
733! PURPOSE
734! -------
735!   The modular radiation scheme is contained in a separate
736!   library. This routine puts the the IFS arrays into appropriate
737!   objects, computing the additional data that is required, and sends
738!   it to the radiation scheme.  It returns net fluxes and surface
739!   flux components needed by the rest of the model.
740!
741!   Lower case is used for variables and types taken from the
742!   radiation library
743!
744! INTERFACE
745! ---------
746!    RADIATION_SCHEME is called from RADLSWR. The
747!    SETUP_RADIATION_SCHEME routine (in the RADIATION_SETUP module)
748!    should have been run first.
749!
750! AUTHOR
751! ------
752!   Robin Hogan, ECMWF
753!   Original: 2015-09-16
754!
755! MODIFICATIONS
756! -------------
757!
758! TO DO
759! -----
760!
761!-----------------------------------------------------------------------
762
763! Modules from ifs or ifsaux libraries
764USE PARKIND1 , ONLY : JPIM, JPRB
765USE YOMHOOK  , ONLY : LHOOK, DR_HOOK
766USE RADIATION_SETUP
767USE YOMCST   , ONLY : RSIGMA ! Stefan-Boltzmann constant
768!USE RADIATION_SETUP, ONLY : SETUP_RADIATION_SCHEME, &
769!                         &  config_type, driver_config_type, &
770!                         &  NWEIGHT_UV,  IBAND_UV,  WEIGHT_UV, &
771!                         &  NWEIGHT_PAR, IBAND_PAR, WEIGHT_PAR, &
772!                         &  ITYPE_TROP_BG_AER,  TROP_BG_AER_MASS_EXT, &
773!                         &  ITYPE_STRAT_BG_AER, STRAT_BG_AER_MASS_EXT, &
774!                         &  ISolverSpartacus
775
776! Modules from radiation library
777USE radiation_single_level,   ONLY : single_level_type
778USE radiation_thermodynamics, ONLY : thermodynamics_type
779USE radiation_gas
780USE radiation_cloud,          ONLY : cloud_type
781USE radiation_aerosol,        ONLY : aerosol_type
782USE radiation_flux,           ONLY : flux_type
783USE radiation_interface,      ONLY : radiation, set_gas_units
784USE radiation_save,           ONLY : save_inputs
785
786USE mod_phys_lmdz_para
787
788IMPLICIT NONE
789
790! INPUT ARGUMENTS
791! *** Array dimensions and ranges
792INTEGER(KIND=JPIM),INTENT(IN) :: KIDIA    ! Start column to process
793INTEGER(KIND=JPIM),INTENT(IN) :: KFDIA    ! End column to process
794!INTEGER, INTENT(IN) :: KIDIA, KFDIA
795INTEGER(KIND=JPIM),INTENT(IN) :: KLON     ! Number of columns
796INTEGER(KIND=JPIM),INTENT(IN) :: KLEV     ! Number of levels
797!INTEGER, INTENT(IN) :: KLON, KLEV
798!INTEGER(KIND=JPIM),INTENT(IN) :: KAEROLMDZ ! Number of aerosol types
799INTEGER(KIND=JPIM),INTENT(IN) :: KAEROSOL
800INTEGER(KIND=JPIM),INTENT(IN) :: NSW ! Numbe of bands
801
802! AI ATTENTION
803!INTEGER, PARAMETER :: KAEROSOL = 12
804
805! *** Single-level fields
806REAL(KIND=JPRB),   INTENT(IN) :: PSOLAR_IRRADIANCE ! (W m-2)
807REAL(KIND=JPRB),   INTENT(IN) :: PMU0(KLON) ! Cosine of solar zenith ang
808REAL(KIND=JPRB),   INTENT(IN) :: PTEMPERATURE_SKIN(KLON) ! (K)
809! Diffuse and direct components of surface shortwave albedo
810!REAL(KIND=JPRB),   INTENT(IN) :: PALBEDO_DIF(KLON,YRERAD%NSW)
811!REAL(KIND=JPRB),   INTENT(IN) :: PALBEDO_DIR(KLON,YRERAD%NSW)
812REAL(KIND=JPRB),   INTENT(IN) :: PALBEDO_DIF(KLON,NSW)
813REAL(KIND=JPRB),   INTENT(IN) :: PALBEDO_DIR(KLON,NSW)
814! Longwave emissivity outside and inside the window region
815REAL(KIND=JPRB),   INTENT(IN) :: PEMIS(KLON)
816REAL(KIND=JPRB),   INTENT(IN) :: PEMIS_WINDOW(KLON)
817! Longitude (radians), sine of latitude
818REAL(KIND=JPRB),   INTENT(IN) :: PGELAM(KLON)
819REAL(KIND=JPRB),   INTENT(IN) :: PGEMU(KLON)
820! Land-sea mask
821!REAL(KIND=JPRB),   INTENT(IN) :: PLAND_SEA_MASK(KLON)
822
823! *** Variables on half levels
824REAL(KIND=JPRB),   INTENT(IN) :: PPRESSURE_H(KLON,KLEV+1)    ! (Pa)
825REAL(KIND=JPRB),   INTENT(IN) :: PTEMPERATURE_H(KLON,KLEV+1) ! (K)
826
827! *** Gas mass mixing ratios on full levels
828REAL(KIND=JPRB),   INTENT(IN) :: PQ(KLON,KLEV)
829! AI
830REAL(KIND=JPRB),   INTENT(IN) :: PQSAT(KLON,KLEV)
831REAL(KIND=JPRB),   INTENT(IN) :: PCO2
832REAL(KIND=JPRB),   INTENT(IN) :: PCH4
833REAL(KIND=JPRB),   INTENT(IN) :: PN2O
834REAL(KIND=JPRB),   INTENT(IN) :: PNO2
835REAL(KIND=JPRB),   INTENT(IN) :: PCFC11
836REAL(KIND=JPRB),   INTENT(IN) :: PCFC12
837REAL(KIND=JPRB),   INTENT(IN) :: PHCFC22
838REAL(KIND=JPRB),   INTENT(IN) :: PCCL4
839REAL(KIND=JPRB),   INTENT(IN) :: PO3(KLON,KLEV) ! AI (kg/kg) ATTENTION (Pa*kg/kg)
840REAL(KIND=JPRB),   INTENT(IN) :: PO2
841
842! *** Cloud fraction and hydrometeor mass mixing ratios
843REAL(KIND=JPRB),   INTENT(IN) :: PCLOUD_FRAC(KLON,KLEV)
844REAL(KIND=JPRB),   INTENT(IN) :: PQ_LIQUID(KLON,KLEV)
845REAL(KIND=JPRB),   INTENT(IN) :: PQ_ICE(KLON,KLEV)
846!REAL(KIND=JPRB),   INTENT(IN) :: PQ_RAIN(KLON,KLEV)
847REAL(KIND=JPRB),   INTENT(IN) :: PQ_SNOW(KLON,KLEV)
848
849! *** Aerosol mass mixing ratios
850REAL(KIND=JPRB),   INTENT(IN) :: PAEROSOL_OLD(KLON,6,KLEV)
851REAL(KIND=JPRB),   INTENT(IN) :: PAEROSOL(KLON,KLEV,KAEROSOL)
852
853!REAL(KIND=JPRB),   INTENT(IN) :: PCCN_LAND(KLON)
854!REAL(KIND=JPRB),   INTENT(IN) :: PCCN_SEA(KLON)
855
856!AI mars 2021
857INTEGER(KIND=JPIM), INTENT(IN) :: IDAY
858REAL(KIND=JPRB), INTENT(IN)    :: TIME
859
860! Name of file names specified on command line
861character(len=512), INTENT(IN) :: namelist_file
862logical, INTENT(IN)            :: ok_3Deffect, debut, ok_volcan
863INTEGER(KIND=JPIM), INTENT(IN) :: flag_aerosol_strat
864
865
866! OUTPUT ARGUMENTS
867
868! *** Net fluxes on half-levels (W m-2)
869REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW(KLON,KLEV+1)
870REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW(KLON,KLEV+1)
871REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_CLEAR(KLON,KLEV+1)
872REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_CLEAR(KLON,KLEV+1)
873
874!*** DN and UP flux on half-levels (W m-2)
875REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_DN(KLON,KLEV+1)
876REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_DN(KLON,KLEV+1)
877REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_DN_CLEAR(KLON,KLEV+1)
878REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_DN_CLEAR(KLON,KLEV+1)
879REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_UP(KLON,KLEV+1)
880REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_UP(KLON,KLEV+1)
881REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_UP_CLEAR(KLON,KLEV+1)
882REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_UP_CLEAR(KLON,KLEV+1)
883
884! Direct component of surface flux into horizontal plane
885REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_DIR(KLON)
886REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_DIR_CLEAR(KLON)
887! As PFLUX_DIR but into a plane perpendicular to the sun
888REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_DIR_INTO_SUN(KLON)
889
890! *** Ultraviolet and photosynthetically active radiation (W m-2)
891REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_UV(KLON)
892REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_PAR(KLON)
893REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_PAR_CLEAR(KLON)
894
895! Diagnosed longwave surface emissivity across the whole spectrum
896REAL(KIND=JPRB),  INTENT(OUT) :: PEMIS_OUT(KLON)   
897
898! Partial derivative of total-sky longwave upward flux at each level
899! with respect to upward flux at surface, used to correct heating
900! rates at gridpoints/timesteps between calls to the full radiation
901! scheme.  Note that this version uses the convention of level index
902! increasing downwards, unlike the local variable ZLwDerivative that
903! is returned from the LW radiation scheme.
904REAL(KIND=JPRB),  INTENT(OUT) :: PLWDERIVATIVE(KLON,KLEV+1)
905
906! Surface diffuse and direct downwelling shortwave flux in each
907! shortwave albedo band, used in RADINTG to update the surface fluxes
908! accounting for high-resolution albedo information
909REAL(KIND=JPRB),  INTENT(OUT) :: PSWDIFFUSEBAND(KLON,NSW)
910REAL(KIND=JPRB),  INTENT(OUT) :: PSWDIRECTBAND (KLON,NSW)
911
912!AI Nov 2023
913REAL(KIND=JPRB),  INTENT(OUT) :: ecrad_cloud_cover_sw(KLON)
914
915! LOCAL VARIABLES
916! AI ATTENTION
917type(config_type),save         :: rad_config
918!!$OMP THREADPRIVATE(rad_config)
919type(driver_config_type),save  :: driver_config
920!!$OMP THREADPRIVATE(driver_config)
921!type(config_type)        :: rad_config
922!type(driver_config_type)  :: driver_config
923TYPE(single_level_type)   :: single_level
924TYPE(thermodynamics_type) :: thermodynamics
925TYPE(gas_type)            :: gas
926TYPE(cloud_type)          :: cloud
927TYPE(aerosol_type)        :: aerosol
928TYPE(flux_type)           :: flux
929
930! Mass mixing ratio of ozone (kg/kg)
931REAL(KIND=JPRB)           :: ZO3(KLON,KLEV)
932
933! Cloud effective radii in microns
934REAL(KIND=JPRB)           :: ZRE_LIQUID_UM(KLON,KLEV)
935REAL(KIND=JPRB)           :: ZRE_ICE_UM(KLON,KLEV)
936
937! Cloud overlap decorrelation length for cloud boundaries in km
938REAL(KIND=JPRB)           :: ZDECORR_LEN_KM(KLON)
939
940! Ratio of cloud overlap decorrelation length for cloud water
941! inhomogeneities to that for cloud boundaries (typically 0.5)
942!REAL(KIND=JPRB)           :: ZDECORR_LEN_RATIO = 0.5_jprb
943
944! The surface net longwave flux if the surface was a black body, used
945! to compute the effective broadband surface emissivity
946REAL(KIND=JPRB)           :: ZBLACK_BODY_NET_LW(KIDIA:KFDIA)
947
948! Layer mass in kg m-2
949REAL(KIND=JPRB)           :: ZLAYER_MASS(KIDIA:KFDIA,KLEV)
950
951! Time integers
952INTEGER :: ITIM
953
954! Loop indices
955INTEGER :: JLON, JLEV, JBAND, JB_ALBEDO, JAER
956
957REAL(KIND=JPRB) :: ZHOOK_HANDLE
958
959! AI ATTENTION traitement aerosols
960INTEGER, PARAMETER :: NAERMACC = 1
961
962logical :: loutput=.true.
963logical :: lprint_input=.false.
964logical :: lprint_config=.false.
965logical, save :: debut_ecrad=.true.
966!$OMP THREADPRIVATE(debut_ecrad)
967integer, save :: itap_ecrad=0
968!$OMP THREADPRIVATE(itap_ecrad)
969
970REAL(KIND=JPRB) ::  inv_cloud_effective_size(KLON,KLEV)
971REAL(KIND=JPRB) ::  inv_inhom_effective_size(KLON,KLEV)
972
973integer :: irang
974
975
976IF (LHOOK) CALL DR_HOOK('RADIATION_SCHEME',0,ZHOOK_HANDLE)
977
978   print*,'Entree radiation_scheme_s2, ok_3Deffect, namelist_file = ', &
979          ok_3Deffect, namelist_file
980! A.I juillet 2023 :
981! Initialisation dans radiation_setup au 1er passage dans Ecrad
982!$OMP MASTER
983!if (.not.ok_3Deffect) then
984  if (debut_ecrad) then
985   call SETUP_RADIATION_SCHEME(loutput,namelist_file,rad_config,driver_config)
986   debut_ecrad=.false.
987  endif
988!else
989!   call SETUP_RADIATION_SCHEME(loutput,namelist_file,rad_config,driver_config)
990!endif
991!$OMP END MASTER
992!$OMP BARRIER
993! Fin partie initialisation et configuration
994
995!AI print fichiers namelist utilise
996!if (is_omp_root) then
997! itap_ecrad=itap_ecrad+1
998! print*,'Dans radiation_scheme itap_ecrad, mpi_rank, omp_rank, namelist_file : ', &
999!         itap_ecrad, mpi_rank, omp_rank, namelist_file
1000!else
1001! print*,'mpi_rank omp_rank, namelist_file :', mpi_rank, omp_rank, namelist_file
1002!endif
1003
1004! AI 11 23 Allocates depplaces au debut
1005print*,'*********** ALLOCATES *******************************'
1006! AI ATTENTION
1007! Allocate memory in radiation objects
1008! emissivite avec une seule bande
1009CALL single_level%allocate(KLON, NSW, 1, &
1010     &                     use_sw_albedo_direct=.TRUE.)
1011CALL thermodynamics%allocate(KLON, KLEV, use_h2o_sat=.true.)
1012CALL cloud%allocate(KLON, KLEV)
1013CALL aerosol%allocate(KLON, 1, KLEV, KAEROSOL)
1014CALL gas%allocate(KLON, KLEV)
1015CALL flux%allocate(rad_config, 1, KLON, KLEV)
1016
1017print*,'************* THERMO (input) ************************************'
1018! Set thermodynamic profiles: simply copy over the half-level
1019! pressure and temperature
1020! AI
1021! pressure_hl > paprs
1022! temperature_hl calculee dans radlsw de la meme facon que pour RRTM
1023thermodynamics%pressure_hl   (KIDIA:KFDIA,:) = PPRESSURE_H   (KIDIA:KFDIA,:)
1024thermodynamics%temperature_hl(KIDIA:KFDIA,:) = PTEMPERATURE_H(KIDIA:KFDIA,:)
1025!print*,'Compute saturation specific humidity'
1026! Compute saturation specific humidity, used to hydrate aerosols. The
1027! "2" for the last argument indicates that the routine is not being
1028! called from within the convection scheme.
1029!CALL SATUR(KIDIA, KFDIA, KLON, 1, KLEV, &
1030!     &  PPRESSURE, PTEMPERATURE, thermodynamics%h2o_sat_liq, 2)
1031! Alternative approximate version using temperature and pressure from
1032! the thermodynamics structure
1033!CALL thermodynamics%calc_saturation_wrt_liquid(KIDIA, KFDIA)
1034!AI ATTENTION
1035thermodynamics%h2o_sat_liq = PQSAT
1036
1037print*,'********** SINGLE LEVEL VARS **********************************'
1038!AI ATTENTION
1039! Set single-level fileds
1040single_level%solar_irradiance              = PSOLAR_IRRADIANCE
1041single_level%cos_sza(KIDIA:KFDIA)          = PMU0(KIDIA:KFDIA)
1042single_level%skin_temperature(KIDIA:KFDIA) = PTEMPERATURE_SKIN(KIDIA:KFDIA)
1043single_level%sw_albedo(KIDIA:KFDIA,:)      = PALBEDO_DIF(KIDIA:KFDIA,:)
1044single_level%sw_albedo_direct(KIDIA:KFDIA,:)=PALBEDO_DIR(KIDIA:KFDIA,:)
1045single_level%lw_emissivity(KIDIA:KFDIA,1)  = PEMIS(KIDIA:KFDIA)
1046!single_level%lw_emissivity(KIDIA:KFDIA,2)  = PEMIS_WINDOW(KIDIA:KFDIA)
1047
1048! Create the relevant seed from date and time get the starting day
1049! and number of minutes since start
1050!IDAY = NDD(NINDAT)
1051!cur_day
1052!ITIM = NINT(NSTEP * YRRIP%TSTEP / 60.0_JPRB)
1053!ITIM = NINT(TIME / 60.0_JPRB)
1054!current_time
1055!allocate(single_level%iseed(KIDIA:KFDIA))
1056!DO JLON = KIDIA, KFDIA
1057  ! This method gives a unique value for roughly every 1-km square
1058  ! on the globe and every minute.  ASIN(PGEMU)*60 gives rough
1059  ! latitude in degrees, which we multiply by 100 to give a unique
1060  ! value for roughly every km. PGELAM*60*100 gives a unique number
1061  ! for roughly every km of longitude around the equator, which we
1062  ! multiply by 180*100 so there is no overlap with the latitude
1063  ! values.  The result can be contained in a 32-byte integer (but
1064  ! since random numbers are generated with the help of integer
1065  ! overflow, it should not matter if the number did overflow).
1066!  single_level%iseed(JLON) = ITIM + IDAY &
1067!       &  +  NINT(PGELAM(JLON)*108000000.0_JPRB &
1068!       &          + ASIN(PGEMU(JLON))*6000.0_JPRB)
1069!ENDDO
1070!AI Nov 23
1071! Simple initialization of the seeds for the Monte Carlo scheme
1072call single_level%init_seed_simple(kidia, kfdia)
1073
1074print*,'********** CLOUDS (allocate + input) *******************************************'
1075!print*,'Appel Allocate clouds'
1076! Set cloud fields
1077cloud%q_liq(KIDIA:KFDIA,:)    = PQ_LIQUID(KIDIA:KFDIA,:)
1078cloud%q_ice(KIDIA:KFDIA,:)    = PQ_ICE(KIDIA:KFDIA,:) + PQ_SNOW(KIDIA:KFDIA,:)
1079cloud%fraction(KIDIA:KFDIA,:) = PCLOUD_FRAC(KIDIA:KFDIA,:)
1080!!! ok AI ATTENTION a voir avec JL
1081! Compute effective radi and convert to metres
1082! AI. : on passe directement les champs de LMDZ
1083cloud%re_liq(KIDIA:KFDIA,:) = ZRE_LIQUID_UM(KIDIA:KFDIA,:)
1084cloud%re_ice(KIDIA:KFDIA,:) = ZRE_ICE_UM(KIDIA:KFDIA,:)
1085! Get the cloud overlap decorrelation length (for cloud boundaries),
1086! in km, according to the parameterization specified by NDECOLAT,
1087! and insert into the "cloud" object. Also get the ratio of
1088! decorrelation lengths for cloud water content inhomogeneities and
1089! cloud boundaries, and set it in the "rad_config" object.
1090! IFS :
1091!CALL CLOUD_OVERLAP_DECORR_LEN(KIDIA, KFDIA, KLON, PGEMU, YRERAD%NDECOLAT, &
1092!     &    ZDECORR_LEN_KM, PDECORR_LEN_RATIO=ZDECORR_LEN_RATIO)
1093! AI valeur dans namelist
1094! rad_config%cloud_inhom_decorr_scaling = ZDECORR_LEN_RATIO
1095!AI ATTENTION meme valeur que dans offline
1096! A mettre dans namelist
1097ZDECORR_LEN_KM = driver_config%overlap_decorr_length
1098DO JLON = KIDIA,KFDIA
1099  CALL cloud%set_overlap_param(thermodynamics, &
1100       &                 ZDECORR_LEN_KM(JLON), &
1101       &                 istartcol=JLON, iendcol=JLON)
1102ENDDO
1103! IFS :
1104! Cloud water content fractional standard deviation is configurable
1105! from namelist NAERAD but must be globally constant. Before it was
1106! hard coded at 1.0.
1107!CALL cloud%create_fractional_std(KLON, KLEV, YRERAD%RCLOUD_FRAC_STD)
1108! AI ATTENTION frac_std=0.75 meme valeur que dans la version offline
1109CALL cloud%create_fractional_std(KLON, KLEV, driver_config%frac_std)
1110
1111!if (ok_3Deffect) then
1112! if (driver_config%ok_effective_size) then
1113!     call cloud%create_inv_cloud_effective_size_eta(klon, klev, &
1114!               &  thermodynamics%pressure_hl, &
1115!               &  driver_config%low_inv_effective_size, &
1116!               &  driver_config%middle_inv_effective_size, &
1117!               &  driver_config%high_inv_effective_size, 0.8_jprb, 0.45_jprb, &
1118!               &  KIDIA, KFDIA)
1119!   else if (driver_config%ok_separation) then
1120!     call cloud%param_cloud_effective_separation_eta(klon, klev, &
1121!               &  thermodynamics%pressure_hl, &
1122!               &  driver_config%cloud_separation_scale_surface, &
1123!               &  driver_config%cloud_separation_scale_toa, &
1124!               &  driver_config%cloud_separation_scale_power, &
1125!               &  driver_config%cloud_inhom_separation_factor, &
1126!               &  KIDIA, KFDIA)
1127!   endif     
1128! else 
1129 if (rad_config%i_solver_sw == ISolverSPARTACUS &
1130      & .or.   rad_config%i_solver_lw == ISolverSPARTACUS) then
1131   ! AI ! Read cloud properties needed by SPARTACUS
1132   if (driver_config%ok_effective_size) then
1133     call cloud%create_inv_cloud_effective_size_eta(klon, klev, &
1134               &  thermodynamics%pressure_hl, &
1135               &  driver_config%low_inv_effective_size, &
1136               &  driver_config%middle_inv_effective_size, &
1137               &  driver_config%high_inv_effective_size, 0.8_jprb, 0.45_jprb, &
1138               &  KIDIA, KFDIA)
1139   else if (driver_config%ok_separation) then
1140     call cloud%param_cloud_effective_separation_eta(klon, klev, &
1141               &  thermodynamics%pressure_hl, &
1142               &  driver_config%cloud_separation_scale_surface, &
1143               &  driver_config%cloud_separation_scale_toa, &
1144               &  driver_config%cloud_separation_scale_power, &
1145               &  driver_config%cloud_inhom_separation_factor, &
1146               &  KIDIA, KFDIA)
1147   endif
1148 endif 
1149!endif
1150
1151print*,'******** AEROSOLS (input) **************************************'
1152!IF (NAERMACC > 0) THEN
1153!ELSE
1154!  CALL aerosol%allocate(KLON, 1, KLEV, 6) ! Tegen climatology
1155!ENDIF
1156! Compute the dry mass of each layer neglecting humidity effects, in
1157! kg m-2, needed to scale some of the aerosol inputs
1158! AI commente ATTENTION
1159!CALL thermodynamics%get_layer_mass(ZLAYER_MASS)
1160
1161! Copy over aerosol mass mixing ratio
1162!IF (NAERMACC > 0) THEN
1163
1164  ! MACC aerosol climatology - this is already in mass mixing ratio
1165  ! units with the required array orientation so we can copy it over
1166  ! directly
1167  aerosol%mixing_ratio(KIDIA:KFDIA,:,:) = PAEROSOL(KIDIA:KFDIA,:,:)
1168
1169  ! Add the tropospheric and stratospheric backgrounds contained in the
1170  ! old Tegen arrays - this is very ugly!
1171! AI ATTENTION
1172!  IF (TROP_BG_AER_MASS_EXT > 0.0_JPRB) THEN
1173!    aerosol%mixing_ratio(KIDIA:KFDIA,:,ITYPE_TROP_BG_AER) &
1174!         &  = aerosol%mixing_ratio(KIDIA:KFDIA,:,ITYPE_TROP_BG_AER) &
1175!         &  + PAEROSOL_OLD(KIDIA:KFDIA,1,:) &
1176!         &  / (ZLAYER_MASS * TROP_BG_AER_MASS_EXT)
1177!  ENDIF
1178!  IF (STRAT_BG_AER_MASS_EXT > 0.0_JPRB) THEN
1179!    aerosol%mixing_ratio(KIDIA:KFDIA,:,ITYPE_STRAT_BG_AER) &
1180!         &  = aerosol%mixing_ratio(KIDIA:KFDIA,:,ITYPE_STRAT_BG_AER) &
1181!         &  + PAEROSOL_OLD(KIDIA:KFDIA,6,:) &
1182!         &  / (ZLAYER_MASS * STRAT_BG_AER_MASS_EXT)
1183!  ENDIF
1184
1185!ELSE
1186
1187  ! Tegen aerosol climatology - the array PAEROSOL_OLD contains the
1188  ! 550-nm optical depth in each layer. The optics data file
1189  ! aerosol_ifs_rrtm_tegen.nc does not contain mass extinction
1190  ! coefficient, but a scaling factor that the 550-nm optical depth
1191  ! should be multiplied by to obtain the optical depth in each
1192  ! spectral band.  Therefore, in order for the units to work out, we
1193  ! need to divide by the layer mass (in kg m-2) to obtain the 550-nm
1194  ! cross-section per unit mass of dry air (so in m2 kg-1).  We also
1195  ! need to permute the array.
1196!  DO JLEV = 1,KLEV
1197!    DO JAER = 1,6
1198!      aerosol%mixing_ratio(KIDIA:KFDIA,JLEV,JAER) &
1199!         &  = PAEROSOL_OLD(KIDIA:KFDIA,JAER,JLEV) &
1200!         &  / ZLAYER_MASS(KIDIA:KFDIA,JLEV)
1201!    ENDDO
1202!  ENDDO
1203!ENDIF
1204
1205print*,'********** GAS (input) ************************************************'
1206!print*,'Appel Allocate gas'
1207! Convert ozone Pa*kg/kg to kg/kg
1208! AI ATTENTION
1209!DO JLEV = 1,KLEV
1210!  DO JLON = KIDIA,KFDIA
1211!    ZO3(JLON,JLEV) = PO3_DP(JLON,JLEV) &
1212!         &         / (PPRESSURE_H(JLON,JLEV+1)-PPRESSURE_H(JLON,JLEV))
1213!  ENDDO
1214!ENDDO
1215!  Insert gas mixing ratios
1216!print*,'Insert gas mixing ratios'
1217CALL gas%put(IH2O,    IMassMixingRatio, PQ)
1218CALL gas%put(IO3,     IMassMixingRatio, PO3)
1219CALL gas%put_well_mixed(ICO2,    IMAssMixingRatio, PCO2)
1220CALL gas%put_well_mixed(ICH4,    IMassMixingRatio, PCH4)
1221CALL gas%put_well_mixed(IN2O,    IMassMixingRatio, PN2O)
1222CALL gas%put_well_mixed(ICFC11,  IMassMixingRatio, PCFC11)
1223CALL gas%put_well_mixed(ICFC12,  IMassMixingRatio, PCFC12)
1224CALL gas%put_well_mixed(IHCFC22, IMassMixingRatio, PHCFC22)
1225CALL gas%put_well_mixed(ICCL4,   IMassMixingRatio, PCCL4)
1226CALL gas%put_well_mixed(IO2,     IMassMixingRatio, PO2)
1227! Ensure the units of the gas mixing ratios are what is required by
1228! the gas absorption model
1229call set_gas_units(rad_config, gas)
1230
1231! Call radiation scheme
1232!print*,'*** Appel radiation *** namelist **** omp_rank ****', &
1233!         omp_rank, namelist_file
1234!   if (rad_config%i_solver_sw == ISolverSPARTACUS) then
1235!    if (driver_config%ok_separation) then
1236!      print*,'Avant radiation, mpi_rank, omp_rank, size, chape inv_cloud = ',&
1237!              mpi_rank, omp_rank, &
1238!              shape(cloud%inv_cloud_effective_size), &
1239!              size(cloud%inv_cloud_effective_size)     
1240!      do jlon=KIDIA, KFDIA     
1241!        do jlev=1,klev
1242!         print*,' Avant radiation mpi_rank, omp_rank, jlon, jlev, &
1243!            &     cloud%inv_cloud_effective_size =', mpi_rank, &
1244!            &     omp_rank, jlon, jlev, &
1245!            &     cloud%inv_cloud_effective_size(jlon,jlev)
1246!        enddo
1247!      enddo
1248!       cloud%inv_cloud_effective_size=inv_cloud_effective_size
1249!       cloud%inv_inhom_effective_size=inv_inhom_effective_size
1250!    endif
1251!   endif
1252CALL radiation(KLON, KLEV, KIDIA, KFDIA, rad_config, &
1253     &  single_level, thermodynamics, gas, cloud, aerosol, flux)
1254
1255 if (rad_config%use_aerosols) then
1256    if (rad_config%i_gas_model_sw == IGasModelIFSRRTMG .or. &
1257     &    rad_config%i_gas_model_lw == IGasModelIFSRRTMG) then     
1258       CALL aeropt_5wv_ecrad(kidia, kfdia, 1, klev, &
1259                             rad_config,thermodynamics,aerosol)
1260    endif                 
1261    if (flag_aerosol_strat.eq.2) then
1262       CALL readaerosolstrato_ecrad(rad_config, debut, ok_volcan)
1263    endif
1264 endif               
1265
1266print*,'*********** Sortie flux ****************'
1267! Cloud cover
1268ecrad_cloud_cover_sw = flux%cloud_cover_sw
1269! Compute required output fluxes
1270! DN and UP flux
1271PFLUX_SW_DN(KIDIA:KFDIA,:) = flux%sw_dn(KIDIA:KFDIA,:)
1272PFLUX_SW_UP(KIDIA:KFDIA,:) = flux%sw_up(KIDIA:KFDIA,:)
1273PFLUX_LW_DN(KIDIA:KFDIA,:) = flux%lw_dn(KIDIA:KFDIA,:)
1274PFLUX_LW_UP(KIDIA:KFDIA,:) = flux%lw_up(KIDIA:KFDIA,:)
1275PFLUX_SW_DN_CLEAR(KIDIA:KFDIA,:) = flux%sw_dn_clear(KIDIA:KFDIA,:)
1276PFLUX_SW_UP_CLEAR(KIDIA:KFDIA,:) = flux%sw_up_clear(KIDIA:KFDIA,:)
1277PFLUX_LW_DN_CLEAR(KIDIA:KFDIA,:) = flux%lw_dn_clear(KIDIA:KFDIA,:)
1278PFLUX_LW_UP_CLEAR(KIDIA:KFDIA,:) = flux%lw_up_clear(KIDIA:KFDIA,:)
1279! First the net fluxes
1280PFLUX_SW(KIDIA:KFDIA,:) = flux%sw_dn(KIDIA:KFDIA,:) - flux%sw_up(KIDIA:KFDIA,:)
1281PFLUX_LW(KIDIA:KFDIA,:) = flux%lw_dn(KIDIA:KFDIA,:) - flux%lw_up(KIDIA:KFDIA,:)
1282PFLUX_SW_CLEAR(KIDIA:KFDIA,:) &
1283     &  = flux%sw_dn_clear(KIDIA:KFDIA,:) - flux%sw_up_clear(KIDIA:KFDIA,:)
1284PFLUX_LW_CLEAR(KIDIA:KFDIA,:) &
1285     &  = flux%lw_dn_clear(KIDIA:KFDIA,:) - flux%lw_up_clear(KIDIA:KFDIA,:)
1286! Now the surface fluxes
1287!PFLUX_SW_DN_SURF(KIDIA:KFDIA) = flux%sw_dn(KIDIA:KFDIA,KLEV+1)
1288!PFLUX_LW_DN_SURF(KIDIA:KFDIA) = flux%lw_dn(KIDIA:KFDIA,KLEV+1)
1289!PFLUX_SW_UP_SURF(KIDIA:KFDIA) = flux%sw_up(KIDIA:KFDIA,KLEV+1)
1290!PFLUX_LW_UP_SURF(KIDIA:KFDIA) = flux%lw_up(KIDIA:KFDIA,KLEV+1)
1291!PFLUX_SW_DN_CLEAR_SURF(KIDIA:KFDIA) = flux%sw_dn_clear(KIDIA:KFDIA,KLEV+1)
1292!PFLUX_LW_DN_CLEAR_SURF(KIDIA:KFDIA) = flux%lw_dn_clear(KIDIA:KFDIA,KLEV+1)
1293!PFLUX_SW_UP_CLEAR_SURF(KIDIA:KFDIA) = flux%sw_up_clear(KIDIA:KFDIA,KLEV+1)
1294!PFLUX_LW_UP_CLEAR_SURF(KIDIA:KFDIA) = flux%lw_up_clear(KIDIA:KFDIA,KLEV+1)
1295PFLUX_DIR(KIDIA:KFDIA) = flux%sw_dn_direct(KIDIA:KFDIA,KLEV+1)
1296PFLUX_DIR_CLEAR(KIDIA:KFDIA) = flux%sw_dn_direct_clear(KIDIA:KFDIA,KLEV+1)
1297PFLUX_DIR_INTO_SUN(KIDIA:KFDIA) = 0.0_JPRB
1298WHERE (PMU0(KIDIA:KFDIA) > EPSILON(1.0_JPRB))
1299  PFLUX_DIR_INTO_SUN(KIDIA:KFDIA) = PFLUX_DIR(KIDIA:KFDIA) / PMU0(KIDIA:KFDIA)
1300END WHERE
1301! Top-of-atmosphere downwelling flux
1302!PFLUX_SW_DN_TOA(KIDIA:KFDIA) = flux%sw_dn(KIDIA:KFDIA,1)
1303!PFLUX_SW_UP_TOA(KIDIA:KFDIA) = flux%sw_up(KIDIA:KFDIA,1)
1304!PFLUX_LW_DN_TOA(KIDIA:KFDIA) = flux%lw_dn(KIDIA:KFDIA,1)
1305!PFLUX_LW_UP_TOA(KIDIA:KFDIA) = flux%lw_up(KIDIA:KFDIA,1)
1306!AI ATTENTION
1307if (0.eq.1) then
1308PFLUX_UV       (KIDIA:KFDIA) = 0.0_JPRB
1309DO JBAND = 1,NWEIGHT_UV
1310  PFLUX_UV(KIDIA:KFDIA) = PFLUX_UV(KIDIA:KFDIA) + WEIGHT_UV(JBAND) &
1311       &  * flux%sw_dn_surf_band(IBAND_UV(JBAND),KIDIA:KFDIA)
1312ENDDO
1313! Compute photosynthetically active radiation similarly
1314PFLUX_PAR      (KIDIA:KFDIA) = 0.0_JPRB
1315PFLUX_PAR_CLEAR(KIDIA:KFDIA) = 0.0_JPRB
1316DO JBAND = 1,NWEIGHT_PAR
1317  PFLUX_PAR(KIDIA:KFDIA) = PFLUX_PAR(KIDIA:KFDIA) + WEIGHT_PAR(JBAND) &
1318       &  * flux%sw_dn_surf_band(IBAND_PAR(JBAND),KIDIA:KFDIA)
1319  PFLUX_PAR_CLEAR(KIDIA:KFDIA) = PFLUX_PAR_CLEAR(KIDIA:KFDIA) &
1320       &  + WEIGHT_PAR(JBAND) &
1321       &  * flux%sw_dn_surf_clear_band(IBAND_PAR(JBAND),KIDIA:KFDIA)
1322ENDDO
1323endif
1324! Compute effective broadband emissivity
1325ZBLACK_BODY_NET_LW = flux%lw_dn(KIDIA:KFDIA,KLEV+1) &
1326     &  - RSIGMA*PTEMPERATURE_SKIN(KIDIA:KFDIA)**4
1327PEMIS_OUT(KIDIA:KFDIA) = PEMIS(KIDIA:KFDIA)
1328WHERE (ABS(ZBLACK_BODY_NET_LW) > 1.0E-5)
1329  PEMIS_OUT(KIDIA:KFDIA) = PFLUX_LW(KIDIA:KFDIA,KLEV+1) / ZBLACK_BODY_NET_LW
1330END WHERE
1331! Copy longwave derivatives
1332! AI ATTENTION
1333!IF (YRERAD%LAPPROXLWUPDATE) THEN
1334IF (rad_config%do_lw_derivatives) THEN
1335  PLWDERIVATIVE(KIDIA:KFDIA,:) = flux%lw_derivatives(KIDIA:KFDIA,:)
1336END IF
1337! Store the shortwave downwelling fluxes in each albedo band
1338!AI ATTENTION
1339!IF (YRERAD%LAPPROXSWUPDATE) THEN
1340if (0.eq.1) then
1341IF (rad_config%do_surface_sw_spectral_flux) THEN
1342  PSWDIFFUSEBAND(KIDIA:KFDIA,:) = 0.0_JPRB
1343  PSWDIRECTBAND (KIDIA:KFDIA,:) = 0.0_JPRB
1344  DO JBAND = 1,rad_config%n_bands_sw
1345    JB_ALBEDO = rad_config%i_albedo_from_band_sw(JBAND)
1346    DO JLON = KIDIA,KFDIA
1347      PSWDIFFUSEBAND(JLON,JB_ALBEDO) = PSWDIFFUSEBAND(JLON,JB_ALBEDO) &
1348           &  + flux%sw_dn_surf_band(JBAND,JLON) &
1349           &  - flux%sw_dn_direct_surf_band(JBAND,JLON)
1350      PSWDIRECTBAND(JLON,JB_ALBEDO)  = PSWDIRECTBAND(JLON,JB_ALBEDO) &
1351           &  + flux%sw_dn_direct_surf_band(JBAND,JLON)
1352    ENDDO
1353  ENDDO
1354ENDIF
1355endif
1356
1357print*,'********** DEALLOCATIONS ************************'
1358CALL single_level%deallocate
1359CALL thermodynamics%deallocate
1360CALL gas%deallocate
1361CALL cloud%deallocate
1362CALL aerosol%deallocate
1363CALL flux%deallocate
1364
1365IF (LHOOK) CALL DR_HOOK('RADIATION_SCHEME',1,ZHOOK_HANDLE)
1366
1367END SUBROUTINE RADIATION_SCHEME_S2
1368
1369end module interface_lmdz_ecrad
Note: See TracBrowser for help on using the repository browser.