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

Last change on this file since 5821 was 5821, checked in by idelkadi, 2 months ago

Added a formulation to prescribe effective cloud size as a hyperbolic tangent function of pressure for calculating radiative fluxes related to 3D cloud effects.
Activation is controlled in namelist_ecrad by the logical key ok_separation_tanh

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