source: LMDZ6/branches/contrails/libf/phylmd/ecrad/lmdz/radiation_scheme_mod.f90 @ 5456

Last change on this file since 5456 was 5268, checked in by abarral, 2 months ago

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