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

Last change on this file since 4791 was 4791, checked in by idelkadi, 4 months ago

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