source: LMDZ6/branches/LMDZ-ECRAD/libf/phylmd/ecrad/radiation_scheme.F90 @ 4171

Last change on this file since 4171 was 4171, checked in by lguez, 2 years ago

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