source: LMDZ6/trunk/libf/phylmd/ecrad/lmdz/radiation_scheme_mod.f90 @ 5655

Last change on this file since 5655 was 5655, checked in by idelkadi, 3 weeks ago

Output of direct radiative fluxes when running LMDZ with Ecrad
File size: 54.4 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,KLEV+1)
163REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_DIR_CLEAR(KLON,KLEV+1)
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)
587! Direct component of flux into horizontal plane
588PFLUX_DIR(KIDIA:KFDIA,:) = flux%sw_dn_direct(KIDIA:KFDIA,:)
589PFLUX_DIR_CLEAR(KIDIA:KFDIA,:) = flux%sw_dn_direct_clear(KIDIA:KFDIA,:)
590PFLUX_DIR_INTO_SUN(KIDIA:KFDIA) = 0.0_JPRB
591WHERE (PMU0(KIDIA:KFDIA) > EPSILON(1.0_JPRB))
592! Direct Surface component of flux into a plane perpendicular to the sun       
593  PFLUX_DIR_INTO_SUN(KIDIA:KFDIA) = PFLUX_DIR(KIDIA:KFDIA,KLEV+1) / PMU0(KIDIA:KFDIA)
594END WHERE
595! Top-of-atmosphere downwelling flux
596!PFLUX_SW_DN_TOA(KIDIA:KFDIA) = flux%sw_dn(KIDIA:KFDIA,1)
597!PFLUX_SW_UP_TOA(KIDIA:KFDIA) = flux%sw_up(KIDIA:KFDIA,1)
598!PFLUX_LW_DN_TOA(KIDIA:KFDIA) = flux%lw_dn(KIDIA:KFDIA,1)
599!PFLUX_LW_UP_TOA(KIDIA:KFDIA) = flux%lw_up(KIDIA:KFDIA,1)
600!AI ATTENTION
601if (0.eq.1) then
602PFLUX_UV       (KIDIA:KFDIA) = 0.0_JPRB
603DO JBAND = 1,NWEIGHT_UV
604  PFLUX_UV(KIDIA:KFDIA) = PFLUX_UV(KIDIA:KFDIA) + WEIGHT_UV(JBAND) &
605       &  * flux%sw_dn_surf_band(IBAND_UV(JBAND),KIDIA:KFDIA)
606ENDDO
607! Compute photosynthetically active radiation similarly
608PFLUX_PAR      (KIDIA:KFDIA) = 0.0_JPRB
609PFLUX_PAR_CLEAR(KIDIA:KFDIA) = 0.0_JPRB
610DO JBAND = 1,NWEIGHT_PAR
611  PFLUX_PAR(KIDIA:KFDIA) = PFLUX_PAR(KIDIA:KFDIA) + WEIGHT_PAR(JBAND) &
612       &  * flux%sw_dn_surf_band(IBAND_PAR(JBAND),KIDIA:KFDIA)
613  PFLUX_PAR_CLEAR(KIDIA:KFDIA) = PFLUX_PAR_CLEAR(KIDIA:KFDIA) &
614       &  + WEIGHT_PAR(JBAND) &
615       &  * flux%sw_dn_surf_clear_band(IBAND_PAR(JBAND),KIDIA:KFDIA)
616ENDDO
617endif
618! Compute effective broadband emissivity
619ZBLACK_BODY_NET_LW = flux%lw_dn(KIDIA:KFDIA,KLEV+1) &
620     &  - RSIGMA*PTEMPERATURE_SKIN(KIDIA:KFDIA)**4
621PEMIS_OUT(KIDIA:KFDIA) = PEMIS(KIDIA:KFDIA)
622WHERE (ABS(ZBLACK_BODY_NET_LW) > 1.0E-5)
623  PEMIS_OUT(KIDIA:KFDIA) = PFLUX_LW(KIDIA:KFDIA,KLEV+1) / ZBLACK_BODY_NET_LW
624END WHERE
625! Copy longwave derivatives
626! AI ATTENTION
627!IF (YRERAD%LAPPROXLWUPDATE) THEN
628IF (rad_config%do_lw_derivatives) THEN
629  PLWDERIVATIVE(KIDIA:KFDIA,:) = flux%lw_derivatives(KIDIA:KFDIA,:)
630END IF
631! Store the shortwave downwelling fluxes in each albedo band
632!AI ATTENTION
633!IF (YRERAD%LAPPROXSWUPDATE) THEN
634if (0.eq.1) then
635IF (rad_config%do_surface_sw_spectral_flux) THEN
636  PSWDIFFUSEBAND(KIDIA:KFDIA,:) = 0.0_JPRB
637  PSWDIRECTBAND (KIDIA:KFDIA,:) = 0.0_JPRB
638  DO JBAND = 1,rad_config%n_bands_sw
639    JB_ALBEDO = rad_config%i_albedo_from_band_sw(JBAND)
640    DO JLON = KIDIA,KFDIA
641      PSWDIFFUSEBAND(JLON,JB_ALBEDO) = PSWDIFFUSEBAND(JLON,JB_ALBEDO) &
642           &  + flux%sw_dn_surf_band(JBAND,JLON) &
643           &  - flux%sw_dn_direct_surf_band(JBAND,JLON)
644      PSWDIRECTBAND(JLON,JB_ALBEDO)  = PSWDIRECTBAND(JLON,JB_ALBEDO) &
645           &  + flux%sw_dn_direct_surf_band(JBAND,JLON)
646    ENDDO
647  ENDDO
648ENDIF
649endif
650
651print*,'********** DEALLOCATIONS ************************'
652CALL single_level%deallocate
653CALL thermodynamics%deallocate
654CALL gas%deallocate
655CALL cloud%deallocate
656CALL aerosol%deallocate
657CALL flux%deallocate
658
659IF (LHOOK) CALL DR_HOOK('RADIATION_SCHEME',1,ZHOOK_HANDLE)
660
661END SUBROUTINE RADIATION_SCHEME
662
663SUBROUTINE RADIATION_SCHEME_S2 &
664! Inputs
665     & (KIDIA, KFDIA, KLON, KLEV, KAEROSOL, NSW, &
666     &  namelist_file, ok_3Deffect, &
667     &  debut, ok_volcan, flag_aerosol_strat, &
668     &  IDAY, TIME, &
669     &  PSOLAR_IRRADIANCE, &
670     &  PMU0, PTEMPERATURE_SKIN, &
671     &  PALBEDO_DIF, PALBEDO_DIR, &
672     &  PEMIS, PEMIS_WINDOW, &
673     &  PGELAM, PGEMU, &
674     &  PPRESSURE_H, PTEMPERATURE_H, PQ, PQSAT, &
675     &  PCO2, PCH4, PN2O, PNO2, PCFC11, PCFC12, PHCFC22, &
676     &  PCCL4, PO3, PO2, &
677     &  PCLOUD_FRAC, PQ_LIQUID, PQ_ICE, PQ_SNOW, &
678     &  ZRE_LIQUID_UM, ZRE_ICE_UM, &
679     &  PAEROSOL_OLD, PAEROSOL, &
680! Outputs
681     &  PFLUX_SW, PFLUX_LW, PFLUX_SW_CLEAR, PFLUX_LW_CLEAR, &
682     &  PFLUX_SW_DN, PFLUX_LW_DN, PFLUX_SW_DN_CLEAR, PFLUX_LW_DN_CLEAR, &
683     &  PFLUX_SW_UP, PFLUX_LW_UP, PFLUX_SW_UP_CLEAR, PFLUX_LW_UP_CLEAR, &
684     &  PFLUX_DIR, PFLUX_DIR_CLEAR, PFLUX_DIR_INTO_SUN, &
685     &  PFLUX_UV, PFLUX_PAR, PFLUX_PAR_CLEAR, &
686     &  PEMIS_OUT, PLWDERIVATIVE, &
687     &  PSWDIFFUSEBAND, PSWDIRECTBAND, &
688     &  ecrad_cloud_cover_sw)
689
690! RADIATION_SCHEME - Interface to modular radiation scheme
691!
692! (C) Copyright 2015- ECMWF.
693!
694! This software is licensed under the terms of the Apache Licence Version 2.0
695! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
696!
697! In applying this licence, ECMWF does not waive the privileges and immunities
698! granted to it by virtue of its status as an intergovernmental organisation
699! nor does it submit to any jurisdiction.
700!
701! PURPOSE
702! -------
703!   The modular radiation scheme is contained in a separate
704!   library. This routine puts the the IFS arrays into appropriate
705!   objects, computing the additional data that is required, and sends
706!   it to the radiation scheme.  It returns net fluxes and surface
707!   flux components needed by the rest of the model.
708!
709!   Lower case is used for variables and types taken from the
710!   radiation library
711!
712! INTERFACE
713! ---------
714!    RADIATION_SCHEME is called from RADLSWR. The
715!    SETUP_RADIATION_SCHEME routine (in the RADIATION_SETUP module)
716!    should have been run first.
717!
718! AUTHOR
719! ------
720!   Robin Hogan, ECMWF
721!   Original: 2015-09-16
722!
723! MODIFICATIONS
724! -------------
725!
726! TO DO
727! -----
728!
729!-----------------------------------------------------------------------
730
731! Modules from ifs or ifsaux libraries
732USE PARKIND1 , ONLY : JPIM, JPRB
733USE YOMHOOK  , ONLY : LHOOK, DR_HOOK
734USE RADIATION_SETUP
735USE YOMCST   , ONLY : RSIGMA ! Stefan-Boltzmann constant
736!USE RADIATION_SETUP, ONLY : SETUP_RADIATION_SCHEME, &
737!                         &  config_type, driver_config_type, &
738!                         &  NWEIGHT_UV,  IBAND_UV,  WEIGHT_UV, &
739!                         &  NWEIGHT_PAR, IBAND_PAR, WEIGHT_PAR, &
740!                         &  ITYPE_TROP_BG_AER,  TROP_BG_AER_MASS_EXT, &
741!                         &  ITYPE_STRAT_BG_AER, STRAT_BG_AER_MASS_EXT, &
742!                         &  ISolverSpartacus
743
744! Modules from radiation library
745USE radiation_single_level,   ONLY : single_level_type
746USE radiation_thermodynamics, ONLY : thermodynamics_type
747USE radiation_gas
748USE radiation_cloud,          ONLY : cloud_type
749USE radiation_aerosol,        ONLY : aerosol_type
750USE radiation_flux,           ONLY : flux_type
751USE radiation_interface,      ONLY : radiation, set_gas_units
752USE radiation_save,           ONLY : save_inputs
753
754USE mod_phys_lmdz_para
755
756IMPLICIT NONE
757
758! INPUT ARGUMENTS
759! *** Array dimensions and ranges
760INTEGER(KIND=JPIM),INTENT(IN) :: KIDIA    ! Start column to process
761INTEGER(KIND=JPIM),INTENT(IN) :: KFDIA    ! End column to process
762!INTEGER, INTENT(IN) :: KIDIA, KFDIA
763INTEGER(KIND=JPIM),INTENT(IN) :: KLON     ! Number of columns
764INTEGER(KIND=JPIM),INTENT(IN) :: KLEV     ! Number of levels
765!INTEGER, INTENT(IN) :: KLON, KLEV
766!INTEGER(KIND=JPIM),INTENT(IN) :: KAEROLMDZ ! Number of aerosol types
767INTEGER(KIND=JPIM),INTENT(IN) :: KAEROSOL
768INTEGER(KIND=JPIM),INTENT(IN) :: NSW ! Numbe of bands
769
770! AI ATTENTION
771!INTEGER, PARAMETER :: KAEROSOL = 12
772
773! *** Single-level fields
774REAL(KIND=JPRB),   INTENT(IN) :: PSOLAR_IRRADIANCE ! (W m-2)
775REAL(KIND=JPRB),   INTENT(IN) :: PMU0(KLON) ! Cosine of solar zenith ang
776REAL(KIND=JPRB),   INTENT(IN) :: PTEMPERATURE_SKIN(KLON) ! (K)
777! Diffuse and direct components of surface shortwave albedo
778!REAL(KIND=JPRB),   INTENT(IN) :: PALBEDO_DIF(KLON,YRERAD%NSW)
779!REAL(KIND=JPRB),   INTENT(IN) :: PALBEDO_DIR(KLON,YRERAD%NSW)
780REAL(KIND=JPRB),   INTENT(IN) :: PALBEDO_DIF(KLON,NSW)
781REAL(KIND=JPRB),   INTENT(IN) :: PALBEDO_DIR(KLON,NSW)
782! Longwave emissivity outside and inside the window region
783REAL(KIND=JPRB),   INTENT(IN) :: PEMIS(KLON)
784REAL(KIND=JPRB),   INTENT(IN) :: PEMIS_WINDOW(KLON)
785! Longitude (radians), sine of latitude
786REAL(KIND=JPRB),   INTENT(IN) :: PGELAM(KLON)
787REAL(KIND=JPRB),   INTENT(IN) :: PGEMU(KLON)
788! Land-sea mask
789!REAL(KIND=JPRB),   INTENT(IN) :: PLAND_SEA_MASK(KLON)
790
791! *** Variables on half levels
792REAL(KIND=JPRB),   INTENT(IN) :: PPRESSURE_H(KLON,KLEV+1)    ! (Pa)
793REAL(KIND=JPRB),   INTENT(IN) :: PTEMPERATURE_H(KLON,KLEV+1) ! (K)
794
795! *** Gas mass mixing ratios on full levels
796REAL(KIND=JPRB),   INTENT(IN) :: PQ(KLON,KLEV)
797! AI
798REAL(KIND=JPRB),   INTENT(IN) :: PQSAT(KLON,KLEV)
799REAL(KIND=JPRB),   INTENT(IN) :: PCO2
800REAL(KIND=JPRB),   INTENT(IN) :: PCH4
801REAL(KIND=JPRB),   INTENT(IN) :: PN2O
802REAL(KIND=JPRB),   INTENT(IN) :: PNO2
803REAL(KIND=JPRB),   INTENT(IN) :: PCFC11
804REAL(KIND=JPRB),   INTENT(IN) :: PCFC12
805REAL(KIND=JPRB),   INTENT(IN) :: PHCFC22
806REAL(KIND=JPRB),   INTENT(IN) :: PCCL4
807REAL(KIND=JPRB),   INTENT(IN) :: PO3(KLON,KLEV) ! AI (kg/kg) ATTENTION (Pa*kg/kg)
808REAL(KIND=JPRB),   INTENT(IN) :: PO2
809
810! *** Cloud fraction and hydrometeor mass mixing ratios
811REAL(KIND=JPRB),   INTENT(IN) :: PCLOUD_FRAC(KLON,KLEV)
812REAL(KIND=JPRB),   INTENT(IN) :: PQ_LIQUID(KLON,KLEV)
813REAL(KIND=JPRB),   INTENT(IN) :: PQ_ICE(KLON,KLEV)
814!REAL(KIND=JPRB),   INTENT(IN) :: PQ_RAIN(KLON,KLEV)
815REAL(KIND=JPRB),   INTENT(IN) :: PQ_SNOW(KLON,KLEV)
816
817! *** Aerosol mass mixing ratios
818REAL(KIND=JPRB),   INTENT(IN) :: PAEROSOL_OLD(KLON,6,KLEV)
819REAL(KIND=JPRB),   INTENT(IN) :: PAEROSOL(KLON,KLEV,KAEROSOL)
820
821!REAL(KIND=JPRB),   INTENT(IN) :: PCCN_LAND(KLON)
822!REAL(KIND=JPRB),   INTENT(IN) :: PCCN_SEA(KLON)
823
824!AI mars 2021
825INTEGER(KIND=JPIM), INTENT(IN) :: IDAY
826REAL(KIND=JPRB), INTENT(IN)    :: TIME
827
828! Name of file names specified on command line
829character(len=512), INTENT(IN) :: namelist_file
830logical, INTENT(IN)            :: ok_3Deffect, debut, ok_volcan
831INTEGER(KIND=JPIM), INTENT(IN) :: flag_aerosol_strat
832
833
834! OUTPUT ARGUMENTS
835
836! *** Net fluxes on half-levels (W m-2)
837REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW(KLON,KLEV+1)
838REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW(KLON,KLEV+1)
839REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_CLEAR(KLON,KLEV+1)
840REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_CLEAR(KLON,KLEV+1)
841
842!*** DN and UP flux on half-levels (W m-2)
843REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_DN(KLON,KLEV+1)
844REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_DN(KLON,KLEV+1)
845REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_DN_CLEAR(KLON,KLEV+1)
846REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_DN_CLEAR(KLON,KLEV+1)
847REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_UP(KLON,KLEV+1)
848REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_UP(KLON,KLEV+1)
849REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_SW_UP_CLEAR(KLON,KLEV+1)
850REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_LW_UP_CLEAR(KLON,KLEV+1)
851
852! Direct component of surface flux into horizontal plane
853REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_DIR(KLON,KLEV+1)
854REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_DIR_CLEAR(KLON,KLEV+1)
855! As PFLUX_DIR but into a plane perpendicular to the sun
856REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_DIR_INTO_SUN(KLON)
857
858! *** Ultraviolet and photosynthetically active radiation (W m-2)
859REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_UV(KLON)
860REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_PAR(KLON)
861REAL(KIND=JPRB),  INTENT(OUT) :: PFLUX_PAR_CLEAR(KLON)
862
863! Diagnosed longwave surface emissivity across the whole spectrum
864REAL(KIND=JPRB),  INTENT(OUT) :: PEMIS_OUT(KLON)   
865
866! Partial derivative of total-sky longwave upward flux at each level
867! with respect to upward flux at surface, used to correct heating
868! rates at gridpoints/timesteps between calls to the full radiation
869! scheme.  Note that this version uses the convention of level index
870! increasing downwards, unlike the local variable ZLwDerivative that
871! is returned from the LW radiation scheme.
872REAL(KIND=JPRB),  INTENT(OUT) :: PLWDERIVATIVE(KLON,KLEV+1)
873
874! Surface diffuse and direct downwelling shortwave flux in each
875! shortwave albedo band, used in RADINTG to update the surface fluxes
876! accounting for high-resolution albedo information
877REAL(KIND=JPRB),  INTENT(OUT) :: PSWDIFFUSEBAND(KLON,NSW)
878REAL(KIND=JPRB),  INTENT(OUT) :: PSWDIRECTBAND (KLON,NSW)
879
880!AI Nov 2023
881REAL(KIND=JPRB),  INTENT(OUT) :: ecrad_cloud_cover_sw(KLON)
882
883! LOCAL VARIABLES
884! AI ATTENTION
885type(config_type),save         :: rad_config
886!!$OMP THREADPRIVATE(rad_config)
887type(driver_config_type),save  :: driver_config
888!!$OMP THREADPRIVATE(driver_config)
889!type(config_type)        :: rad_config
890!type(driver_config_type)  :: driver_config
891TYPE(single_level_type)   :: single_level
892TYPE(thermodynamics_type) :: thermodynamics
893TYPE(gas_type)            :: gas
894TYPE(cloud_type)          :: cloud
895TYPE(aerosol_type)        :: aerosol
896TYPE(flux_type)           :: flux
897
898! Mass mixing ratio of ozone (kg/kg)
899REAL(KIND=JPRB)           :: ZO3(KLON,KLEV)
900
901! Cloud effective radii in microns
902REAL(KIND=JPRB)           :: ZRE_LIQUID_UM(KLON,KLEV)
903REAL(KIND=JPRB)           :: ZRE_ICE_UM(KLON,KLEV)
904
905! Cloud overlap decorrelation length for cloud boundaries in km
906REAL(KIND=JPRB)           :: ZDECORR_LEN_M
907REAL(KIND=JPRB)           :: ZDECORR_LEN_M_1D(KLON)
908REAL(KIND=JPRB)           :: ZDECORR_LEN_M_2D(KLON,KLEV)
909
910! Ratio of cloud overlap decorrelation length for cloud water
911! inhomogeneities to that for cloud boundaries (typically 0.5)
912!REAL(KIND=JPRB)           :: ZDECORR_LEN_RATIO = 0.5_jprb
913
914! The surface net longwave flux if the surface was a black body, used
915! to compute the effective broadband surface emissivity
916REAL(KIND=JPRB)           :: ZBLACK_BODY_NET_LW(KIDIA:KFDIA)
917
918! Layer mass in kg m-2
919REAL(KIND=JPRB)           :: ZLAYER_MASS(KIDIA:KFDIA,KLEV)
920
921! Time integers
922INTEGER :: ITIM
923
924! Loop indices
925INTEGER :: JLON, JLEV, JBAND, JB_ALBEDO, JAER
926
927REAL(KIND=JPRB) :: ZHOOK_HANDLE
928
929! AI ATTENTION traitement aerosols
930INTEGER, PARAMETER :: NAERMACC = 1
931
932logical :: loutput=.true.
933logical :: lprint_input=.false.
934logical :: lprint_config=.false.
935logical, save :: debut_ecrad=.true.
936!$OMP THREADPRIVATE(debut_ecrad)
937integer, save :: itap_ecrad=0
938!$OMP THREADPRIVATE(itap_ecrad)
939
940REAL(KIND=JPRB) ::  inv_cloud_effective_size(KLON,KLEV)
941REAL(KIND=JPRB) ::  inv_inhom_effective_size(KLON,KLEV)
942
943integer :: irang
944
945
946IF (LHOOK) CALL DR_HOOK('RADIATION_SCHEME',0,ZHOOK_HANDLE)
947
948   print*,'Entree radiation_scheme_s2, ok_3Deffect, namelist_file = ', &
949          ok_3Deffect, namelist_file
950! A.I juillet 2023 :
951! Initialisation dans radiation_setup au 1er passage dans Ecrad
952!$OMP MASTER
953!if (.not.ok_3Deffect) then
954  if (debut_ecrad) then
955   call SETUP_RADIATION_SCHEME(loutput,namelist_file,rad_config,driver_config)
956   debut_ecrad=.false.
957  endif
958!else
959!   call SETUP_RADIATION_SCHEME(loutput,namelist_file,rad_config,driver_config)
960!endif
961!$OMP END MASTER
962!$OMP BARRIER
963! Fin partie initialisation et configuration
964
965!AI print fichiers namelist utilise
966!if (is_omp_root) then
967! itap_ecrad=itap_ecrad+1
968! print*,'Dans radiation_scheme itap_ecrad, mpi_rank, omp_rank, namelist_file : ', &
969!         itap_ecrad, mpi_rank, omp_rank, namelist_file
970!else
971! print*,'mpi_rank omp_rank, namelist_file :', mpi_rank, omp_rank, namelist_file
972!endif
973
974! AI 11 23 Allocates depplaces au debut
975print*,'*********** ALLOCATES *******************************'
976! AI ATTENTION
977! Allocate memory in radiation objects
978! emissivite avec une seule bande
979CALL single_level%allocate(KLON, NSW, 1, &
980     &                     use_sw_albedo_direct=.TRUE.)
981CALL thermodynamics%allocate(KLON, KLEV, use_h2o_sat=.true.)
982CALL cloud%allocate(KLON, KLEV)
983CALL aerosol%allocate(KLON, 1, KLEV, KAEROSOL)
984CALL gas%allocate(KLON, KLEV)
985CALL flux%allocate(rad_config, 1, KLON, KLEV)
986
987print*,'************* THERMO (input) ************************************'
988! Set thermodynamic profiles: simply copy over the half-level
989! pressure and temperature
990! AI
991! pressure_hl > paprs
992! temperature_hl calculee dans radlsw de la meme facon que pour RRTM
993thermodynamics%pressure_hl   (KIDIA:KFDIA,:) = PPRESSURE_H   (KIDIA:KFDIA,:)
994thermodynamics%temperature_hl(KIDIA:KFDIA,:) = PTEMPERATURE_H(KIDIA:KFDIA,:)
995!print*,'Compute saturation specific humidity'
996! Compute saturation specific humidity, used to hydrate aerosols. The
997! "2" for the last argument indicates that the routine is not being
998! called from within the convection scheme.
999!CALL SATUR(KIDIA, KFDIA, KLON, 1, KLEV, &
1000!     &  PPRESSURE, PTEMPERATURE, thermodynamics%h2o_sat_liq, 2)
1001! Alternative approximate version using temperature and pressure from
1002! the thermodynamics structure
1003!CALL thermodynamics%calc_saturation_wrt_liquid(KIDIA, KFDIA)
1004!AI ATTENTION
1005thermodynamics%h2o_sat_liq = PQSAT
1006
1007print*,'********** SINGLE LEVEL VARS **********************************'
1008!AI ATTENTION
1009! Set single-level fileds
1010single_level%solar_irradiance              = PSOLAR_IRRADIANCE
1011single_level%cos_sza(KIDIA:KFDIA)          = PMU0(KIDIA:KFDIA)
1012single_level%skin_temperature(KIDIA:KFDIA) = PTEMPERATURE_SKIN(KIDIA:KFDIA)
1013single_level%sw_albedo(KIDIA:KFDIA,:)      = PALBEDO_DIF(KIDIA:KFDIA,:)
1014single_level%sw_albedo_direct(KIDIA:KFDIA,:)=PALBEDO_DIR(KIDIA:KFDIA,:)
1015single_level%lw_emissivity(KIDIA:KFDIA,1)  = PEMIS(KIDIA:KFDIA)
1016!single_level%lw_emissivity(KIDIA:KFDIA,2)  = PEMIS_WINDOW(KIDIA:KFDIA)
1017
1018! Create the relevant seed from date and time get the starting day
1019! and number of minutes since start
1020!IDAY = NDD(NINDAT)
1021!cur_day
1022!ITIM = NINT(NSTEP * YRRIP%TSTEP / 60.0_JPRB)
1023!ITIM = NINT(TIME / 60.0_JPRB)
1024!current_time
1025!allocate(single_level%iseed(KIDIA:KFDIA))
1026!DO JLON = KIDIA, KFDIA
1027  ! This method gives a unique value for roughly every 1-km square
1028  ! on the globe and every minute.  ASIN(PGEMU)*60 gives rough
1029  ! latitude in degrees, which we multiply by 100 to give a unique
1030  ! value for roughly every km. PGELAM*60*100 gives a unique number
1031  ! for roughly every km of longitude around the equator, which we
1032  ! multiply by 180*100 so there is no overlap with the latitude
1033  ! values.  The result can be contained in a 32-byte integer (but
1034  ! since random numbers are generated with the help of integer
1035  ! overflow, it should not matter if the number did overflow).
1036!  single_level%iseed(JLON) = ITIM + IDAY &
1037!       &  +  NINT(PGELAM(JLON)*108000000.0_JPRB &
1038!       &          + ASIN(PGEMU(JLON))*6000.0_JPRB)
1039!ENDDO
1040!AI Nov 23
1041! Simple initialization of the seeds for the Monte Carlo scheme
1042call single_level%init_seed_simple(kidia, kfdia)
1043
1044print*,'********** CLOUDS (allocate + input) *******************************************'
1045!print*,'Appel Allocate clouds'
1046! Set cloud fields
1047cloud%q_liq(KIDIA:KFDIA,:)    = PQ_LIQUID(KIDIA:KFDIA,:)
1048cloud%q_ice(KIDIA:KFDIA,:)    = PQ_ICE(KIDIA:KFDIA,:) + PQ_SNOW(KIDIA:KFDIA,:)
1049cloud%fraction(KIDIA:KFDIA,:) = PCLOUD_FRAC(KIDIA:KFDIA,:)
1050!!! ok AI ATTENTION a voir avec JL
1051! Compute effective radi and convert to metres
1052! AI. : on passe directement les champs de LMDZ
1053cloud%re_liq(KIDIA:KFDIA,:) = ZRE_LIQUID_UM(KIDIA:KFDIA,:)
1054cloud%re_ice(KIDIA:KFDIA,:) = ZRE_ICE_UM(KIDIA:KFDIA,:)
1055! Get the cloud overlap decorrelation length (for cloud boundaries),
1056! in km, according to the parameterization specified by NDECOLAT,
1057! and insert into the "cloud" object. Also get the ratio of
1058! decorrelation lengths for cloud water content inhomogeneities and
1059! cloud boundaries, and set it in the "rad_config" object.
1060! IFS :
1061!CALL CLOUD_OVERLAP_DECORR_LEN(KIDIA, KFDIA, KLON, PGEMU, YRERAD%NDECOLAT, &
1062!     &    ZDECORR_LEN_KM, PDECORR_LEN_RATIO=ZDECORR_LEN_RATIO)
1063! AI valeur dans namelist
1064! rad_config%cloud_inhom_decorr_scaling = ZDECORR_LEN_RATIO
1065!AI ATTENTION meme valeur que dans offline
1066! A mettre dans namelist
1067   CALL CALCUL_CLOUD_OVERLAP_DECORR_LEN &
1068     & (KIDIA, KFDIA, KLON, KLEV, &
1069     &  driver_config, &
1070     &  thermodynamics%pressure_hl , &
1071     &  ZDECORR_LEN_M_2D)
1072
1073 if (driver_config%kdecolat.eq.0) then
1074  ZDECORR_LEN_M = ZDECORR_LEN_M_2D(1,1)
1075  DO JLON = KIDIA,KFDIA
1076   CALL cloud%set_overlap_param(thermodynamics, &
1077       &                 ZDECORR_LEN_M, &
1078       &                 istartcol=JLON, iendcol=JLON)
1079  ENDDO
1080 else if (driver_config%kdecolat.eq.1.or.driver_config%kdecolat.eq.2) then
1081  ZDECORR_LEN_M_1D = ZDECORR_LEN_M_2D(:,1)
1082  DO JLON = KIDIA,KFDIA
1083   CALL cloud%set_overlap_param(thermodynamics, &
1084       &                 ZDECORR_LEN_M_1D, &
1085       &                 istartcol=JLON, iendcol=JLON)
1086  ENDDO
1087 else if (driver_config%kdecolat.eq.3) then
1088  DO JLON = KIDIA,KFDIA
1089   CALL cloud%set_overlap_param_var2D(thermodynamics, &
1090       &                 ZDECORR_LEN_M_2D, KLEV, &
1091       &                 istartcol=JLON, iendcol=JLON)
1092  ENDDO
1093 endif
1094   
1095! IFS :
1096! Cloud water content fractional standard deviation is configurable
1097! from namelist NAERAD but must be globally constant. Before it was
1098! hard coded at 1.0.
1099!CALL cloud%create_fractional_std(KLON, KLEV, YRERAD%RCLOUD_FRAC_STD)
1100! AI ATTENTION frac_std=0.75 meme valeur que dans la version offline
1101CALL cloud%create_fractional_std(KLON, KLEV, driver_config%frac_std)
1102
1103!if (ok_3Deffect) then
1104! if (driver_config%ok_effective_size) then
1105!     call cloud%create_inv_cloud_effective_size_eta(klon, klev, &
1106!               &  thermodynamics%pressure_hl, &
1107!               &  driver_config%low_inv_effective_size, &
1108!               &  driver_config%middle_inv_effective_size, &
1109!               &  driver_config%high_inv_effective_size, 0.8_jprb, 0.45_jprb, &
1110!               &  KIDIA, KFDIA)
1111!   else if (driver_config%ok_separation) then
1112!     call cloud%param_cloud_effective_separation_eta(klon, klev, &
1113!               &  thermodynamics%pressure_hl, &
1114!               &  driver_config%cloud_separation_scale_surface, &
1115!               &  driver_config%cloud_separation_scale_toa, &
1116!               &  driver_config%cloud_separation_scale_power, &
1117!               &  driver_config%cloud_inhom_separation_factor, &
1118!               &  KIDIA, KFDIA)
1119!   endif     
1120! else 
1121 if (rad_config%i_solver_sw == ISolverSPARTACUS &
1122      & .or.   rad_config%i_solver_lw == ISolverSPARTACUS) then
1123   ! AI ! Read cloud properties needed by SPARTACUS
1124   if (driver_config%ok_effective_size) then
1125     call cloud%create_inv_cloud_effective_size_eta(klon, klev, &
1126               &  thermodynamics%pressure_hl, &
1127               &  driver_config%low_inv_effective_size, &
1128               &  driver_config%middle_inv_effective_size, &
1129               &  driver_config%high_inv_effective_size, 0.8_jprb, 0.45_jprb, &
1130               &  KIDIA, KFDIA)
1131   else if (driver_config%ok_separation) then
1132     call cloud%param_cloud_effective_separation_eta(klon, klev, &
1133               &  thermodynamics%pressure_hl, &
1134               &  driver_config%cloud_separation_scale_surface, &
1135               &  driver_config%cloud_separation_scale_toa, &
1136               &  driver_config%cloud_separation_scale_power, &
1137               &  driver_config%cloud_inhom_separation_factor, &
1138               &  KIDIA, KFDIA)
1139   endif
1140 endif 
1141!endif
1142
1143print*,'******** AEROSOLS (input) **************************************'
1144!IF (NAERMACC > 0) THEN
1145!ELSE
1146!  CALL aerosol%allocate(KLON, 1, KLEV, 6) ! Tegen climatology
1147!ENDIF
1148! Compute the dry mass of each layer neglecting humidity effects, in
1149! kg m-2, needed to scale some of the aerosol inputs
1150! AI commente ATTENTION
1151!CALL thermodynamics%get_layer_mass(ZLAYER_MASS)
1152
1153! Copy over aerosol mass mixing ratio
1154!IF (NAERMACC > 0) THEN
1155
1156  ! MACC aerosol climatology - this is already in mass mixing ratio
1157  ! units with the required array orientation so we can copy it over
1158  ! directly
1159  aerosol%mixing_ratio(KIDIA:KFDIA,:,:) = PAEROSOL(KIDIA:KFDIA,:,:)
1160
1161  ! Add the tropospheric and stratospheric backgrounds contained in the
1162  ! old Tegen arrays - this is very ugly!
1163! AI ATTENTION
1164!  IF (TROP_BG_AER_MASS_EXT > 0.0_JPRB) THEN
1165!    aerosol%mixing_ratio(KIDIA:KFDIA,:,ITYPE_TROP_BG_AER) &
1166!         &  = aerosol%mixing_ratio(KIDIA:KFDIA,:,ITYPE_TROP_BG_AER) &
1167!         &  + PAEROSOL_OLD(KIDIA:KFDIA,1,:) &
1168!         &  / (ZLAYER_MASS * TROP_BG_AER_MASS_EXT)
1169!  ENDIF
1170!  IF (STRAT_BG_AER_MASS_EXT > 0.0_JPRB) THEN
1171!    aerosol%mixing_ratio(KIDIA:KFDIA,:,ITYPE_STRAT_BG_AER) &
1172!         &  = aerosol%mixing_ratio(KIDIA:KFDIA,:,ITYPE_STRAT_BG_AER) &
1173!         &  + PAEROSOL_OLD(KIDIA:KFDIA,6,:) &
1174!         &  / (ZLAYER_MASS * STRAT_BG_AER_MASS_EXT)
1175!  ENDIF
1176
1177!ELSE
1178
1179  ! Tegen aerosol climatology - the array PAEROSOL_OLD contains the
1180  ! 550-nm optical depth in each layer. The optics data file
1181  ! aerosol_ifs_rrtm_tegen.nc does not contain mass extinction
1182  ! coefficient, but a scaling factor that the 550-nm optical depth
1183  ! should be multiplied by to obtain the optical depth in each
1184  ! spectral band.  Therefore, in order for the units to work out, we
1185  ! need to divide by the layer mass (in kg m-2) to obtain the 550-nm
1186  ! cross-section per unit mass of dry air (so in m2 kg-1).  We also
1187  ! need to permute the array.
1188!  DO JLEV = 1,KLEV
1189!    DO JAER = 1,6
1190!      aerosol%mixing_ratio(KIDIA:KFDIA,JLEV,JAER) &
1191!         &  = PAEROSOL_OLD(KIDIA:KFDIA,JAER,JLEV) &
1192!         &  / ZLAYER_MASS(KIDIA:KFDIA,JLEV)
1193!    ENDDO
1194!  ENDDO
1195!ENDIF
1196
1197print*,'********** GAS (input) ************************************************'
1198!print*,'Appel Allocate gas'
1199! Convert ozone Pa*kg/kg to kg/kg
1200! AI ATTENTION
1201!DO JLEV = 1,KLEV
1202!  DO JLON = KIDIA,KFDIA
1203!    ZO3(JLON,JLEV) = PO3_DP(JLON,JLEV) &
1204!         &         / (PPRESSURE_H(JLON,JLEV+1)-PPRESSURE_H(JLON,JLEV))
1205!  ENDDO
1206!ENDDO
1207!  Insert gas mixing ratios
1208!print*,'Insert gas mixing ratios'
1209CALL gas%put(IH2O,    IMassMixingRatio, PQ)
1210CALL gas%put(IO3,     IMassMixingRatio, PO3)
1211CALL gas%put_well_mixed(ICO2,    IMAssMixingRatio, PCO2)
1212CALL gas%put_well_mixed(ICH4,    IMassMixingRatio, PCH4)
1213CALL gas%put_well_mixed(IN2O,    IMassMixingRatio, PN2O)
1214CALL gas%put_well_mixed(ICFC11,  IMassMixingRatio, PCFC11)
1215CALL gas%put_well_mixed(ICFC12,  IMassMixingRatio, PCFC12)
1216CALL gas%put_well_mixed(IHCFC22, IMassMixingRatio, PHCFC22)
1217CALL gas%put_well_mixed(ICCL4,   IMassMixingRatio, PCCL4)
1218CALL gas%put_well_mixed(IO2,     IMassMixingRatio, PO2)
1219! Ensure the units of the gas mixing ratios are what is required by
1220! the gas absorption model
1221call set_gas_units(rad_config, gas)
1222
1223! Call radiation scheme
1224!print*,'*** Appel radiation *** namelist **** omp_rank ****', &
1225!         omp_rank, namelist_file
1226!   if (rad_config%i_solver_sw == ISolverSPARTACUS) then
1227!    if (driver_config%ok_separation) then
1228!      print*,'Avant radiation, mpi_rank, omp_rank, size, chape inv_cloud = ',&
1229!              mpi_rank, omp_rank, &
1230!              shape(cloud%inv_cloud_effective_size), &
1231!              size(cloud%inv_cloud_effective_size)     
1232!      do jlon=KIDIA, KFDIA     
1233!        do jlev=1,klev
1234!         print*,' Avant radiation mpi_rank, omp_rank, jlon, jlev, &
1235!            &     cloud%inv_cloud_effective_size =', mpi_rank, &
1236!            &     omp_rank, jlon, jlev, &
1237!            &     cloud%inv_cloud_effective_size(jlon,jlev)
1238!        enddo
1239!      enddo
1240!       cloud%inv_cloud_effective_size=inv_cloud_effective_size
1241!       cloud%inv_inhom_effective_size=inv_inhom_effective_size
1242!    endif
1243!   endif
1244CALL radiation(KLON, KLEV, KIDIA, KFDIA, rad_config, &
1245     &  single_level, thermodynamics, gas, cloud, aerosol, flux)
1246
1247 if (rad_config%use_aerosols) then
1248    if (rad_config%i_gas_model_sw == IGasModelIFSRRTMG .or. &
1249     &    rad_config%i_gas_model_lw == IGasModelIFSRRTMG) then     
1250       CALL aeropt_5wv_ecrad(kidia, kfdia, 1, klev, &
1251                             rad_config,thermodynamics,aerosol)
1252    endif                 
1253    if (flag_aerosol_strat.eq.2) then
1254       CALL readaerosolstrato_ecrad(rad_config, debut, ok_volcan)
1255    endif
1256 endif               
1257
1258print*,'*********** Sortie flux ****************'
1259! Cloud cover
1260ecrad_cloud_cover_sw = flux%cloud_cover_sw
1261! Compute required output fluxes
1262! DN and UP flux
1263PFLUX_SW_DN(KIDIA:KFDIA,:) = flux%sw_dn(KIDIA:KFDIA,:)
1264PFLUX_SW_UP(KIDIA:KFDIA,:) = flux%sw_up(KIDIA:KFDIA,:)
1265PFLUX_LW_DN(KIDIA:KFDIA,:) = flux%lw_dn(KIDIA:KFDIA,:)
1266PFLUX_LW_UP(KIDIA:KFDIA,:) = flux%lw_up(KIDIA:KFDIA,:)
1267PFLUX_SW_DN_CLEAR(KIDIA:KFDIA,:) = flux%sw_dn_clear(KIDIA:KFDIA,:)
1268PFLUX_SW_UP_CLEAR(KIDIA:KFDIA,:) = flux%sw_up_clear(KIDIA:KFDIA,:)
1269PFLUX_LW_DN_CLEAR(KIDIA:KFDIA,:) = flux%lw_dn_clear(KIDIA:KFDIA,:)
1270PFLUX_LW_UP_CLEAR(KIDIA:KFDIA,:) = flux%lw_up_clear(KIDIA:KFDIA,:)
1271! First the net fluxes
1272PFLUX_SW(KIDIA:KFDIA,:) = flux%sw_dn(KIDIA:KFDIA,:) - flux%sw_up(KIDIA:KFDIA,:)
1273PFLUX_LW(KIDIA:KFDIA,:) = flux%lw_dn(KIDIA:KFDIA,:) - flux%lw_up(KIDIA:KFDIA,:)
1274PFLUX_SW_CLEAR(KIDIA:KFDIA,:) &
1275     &  = flux%sw_dn_clear(KIDIA:KFDIA,:) - flux%sw_up_clear(KIDIA:KFDIA,:)
1276PFLUX_LW_CLEAR(KIDIA:KFDIA,:) &
1277     &  = flux%lw_dn_clear(KIDIA:KFDIA,:) - flux%lw_up_clear(KIDIA:KFDIA,:)
1278! Now the surface fluxes
1279!PFLUX_SW_DN_SURF(KIDIA:KFDIA) = flux%sw_dn(KIDIA:KFDIA,KLEV+1)
1280!PFLUX_LW_DN_SURF(KIDIA:KFDIA) = flux%lw_dn(KIDIA:KFDIA,KLEV+1)
1281!PFLUX_SW_UP_SURF(KIDIA:KFDIA) = flux%sw_up(KIDIA:KFDIA,KLEV+1)
1282!PFLUX_LW_UP_SURF(KIDIA:KFDIA) = flux%lw_up(KIDIA:KFDIA,KLEV+1)
1283!PFLUX_SW_DN_CLEAR_SURF(KIDIA:KFDIA) = flux%sw_dn_clear(KIDIA:KFDIA,KLEV+1)
1284!PFLUX_LW_DN_CLEAR_SURF(KIDIA:KFDIA) = flux%lw_dn_clear(KIDIA:KFDIA,KLEV+1)
1285!PFLUX_SW_UP_CLEAR_SURF(KIDIA:KFDIA) = flux%sw_up_clear(KIDIA:KFDIA,KLEV+1)
1286!PFLUX_LW_UP_CLEAR_SURF(KIDIA:KFDIA) = flux%lw_up_clear(KIDIA:KFDIA,KLEV+1)
1287PFLUX_DIR(KIDIA:KFDIA,:) = flux%sw_dn_direct(KIDIA:KFDIA,:)
1288PFLUX_DIR_CLEAR(KIDIA:KFDIA,:) = flux%sw_dn_direct_clear(KIDIA:KFDIA,:)
1289PFLUX_DIR_INTO_SUN(KIDIA:KFDIA) = 0.0_JPRB
1290WHERE (PMU0(KIDIA:KFDIA) > EPSILON(1.0_JPRB))
1291  PFLUX_DIR_INTO_SUN(KIDIA:KFDIA) = PFLUX_DIR(KIDIA:KFDIA,KLEV+1) / PMU0(KIDIA:KFDIA)
1292END WHERE
1293! Top-of-atmosphere downwelling flux
1294!PFLUX_SW_DN_TOA(KIDIA:KFDIA) = flux%sw_dn(KIDIA:KFDIA,1)
1295!PFLUX_SW_UP_TOA(KIDIA:KFDIA) = flux%sw_up(KIDIA:KFDIA,1)
1296!PFLUX_LW_DN_TOA(KIDIA:KFDIA) = flux%lw_dn(KIDIA:KFDIA,1)
1297!PFLUX_LW_UP_TOA(KIDIA:KFDIA) = flux%lw_up(KIDIA:KFDIA,1)
1298!AI ATTENTION
1299if (0.eq.1) then
1300PFLUX_UV       (KIDIA:KFDIA) = 0.0_JPRB
1301DO JBAND = 1,NWEIGHT_UV
1302  PFLUX_UV(KIDIA:KFDIA) = PFLUX_UV(KIDIA:KFDIA) + WEIGHT_UV(JBAND) &
1303       &  * flux%sw_dn_surf_band(IBAND_UV(JBAND),KIDIA:KFDIA)
1304ENDDO
1305! Compute photosynthetically active radiation similarly
1306PFLUX_PAR      (KIDIA:KFDIA) = 0.0_JPRB
1307PFLUX_PAR_CLEAR(KIDIA:KFDIA) = 0.0_JPRB
1308DO JBAND = 1,NWEIGHT_PAR
1309  PFLUX_PAR(KIDIA:KFDIA) = PFLUX_PAR(KIDIA:KFDIA) + WEIGHT_PAR(JBAND) &
1310       &  * flux%sw_dn_surf_band(IBAND_PAR(JBAND),KIDIA:KFDIA)
1311  PFLUX_PAR_CLEAR(KIDIA:KFDIA) = PFLUX_PAR_CLEAR(KIDIA:KFDIA) &
1312       &  + WEIGHT_PAR(JBAND) &
1313       &  * flux%sw_dn_surf_clear_band(IBAND_PAR(JBAND),KIDIA:KFDIA)
1314ENDDO
1315endif
1316! Compute effective broadband emissivity
1317ZBLACK_BODY_NET_LW = flux%lw_dn(KIDIA:KFDIA,KLEV+1) &
1318     &  - RSIGMA*PTEMPERATURE_SKIN(KIDIA:KFDIA)**4
1319PEMIS_OUT(KIDIA:KFDIA) = PEMIS(KIDIA:KFDIA)
1320WHERE (ABS(ZBLACK_BODY_NET_LW) > 1.0E-5)
1321  PEMIS_OUT(KIDIA:KFDIA) = PFLUX_LW(KIDIA:KFDIA,KLEV+1) / ZBLACK_BODY_NET_LW
1322END WHERE
1323! Copy longwave derivatives
1324! AI ATTENTION
1325!IF (YRERAD%LAPPROXLWUPDATE) THEN
1326IF (rad_config%do_lw_derivatives) THEN
1327  PLWDERIVATIVE(KIDIA:KFDIA,:) = flux%lw_derivatives(KIDIA:KFDIA,:)
1328END IF
1329! Store the shortwave downwelling fluxes in each albedo band
1330!AI ATTENTION
1331!IF (YRERAD%LAPPROXSWUPDATE) THEN
1332if (0.eq.1) then
1333IF (rad_config%do_surface_sw_spectral_flux) THEN
1334  PSWDIFFUSEBAND(KIDIA:KFDIA,:) = 0.0_JPRB
1335  PSWDIRECTBAND (KIDIA:KFDIA,:) = 0.0_JPRB
1336  DO JBAND = 1,rad_config%n_bands_sw
1337    JB_ALBEDO = rad_config%i_albedo_from_band_sw(JBAND)
1338    DO JLON = KIDIA,KFDIA
1339      PSWDIFFUSEBAND(JLON,JB_ALBEDO) = PSWDIFFUSEBAND(JLON,JB_ALBEDO) &
1340           &  + flux%sw_dn_surf_band(JBAND,JLON) &
1341           &  - flux%sw_dn_direct_surf_band(JBAND,JLON)
1342      PSWDIRECTBAND(JLON,JB_ALBEDO)  = PSWDIRECTBAND(JLON,JB_ALBEDO) &
1343           &  + flux%sw_dn_direct_surf_band(JBAND,JLON)
1344    ENDDO
1345  ENDDO
1346ENDIF
1347endif
1348
1349print*,'********** DEALLOCATIONS ************************'
1350CALL single_level%deallocate
1351CALL thermodynamics%deallocate
1352CALL gas%deallocate
1353CALL cloud%deallocate
1354CALL aerosol%deallocate
1355CALL flux%deallocate
1356
1357IF (LHOOK) CALL DR_HOOK('RADIATION_SCHEME',1,ZHOOK_HANDLE)
1358
1359END SUBROUTINE RADIATION_SCHEME_S2
1360
1361end module interface_lmdz_ecrad
Note: See TracBrowser for help on using the repository browser.