source: LMDZ6/branches/LMDZ_ECRad/libf/phylmd/ecrad/radiation_scheme.F90 @ 4188

Last change on this file since 4188 was 4188, checked in by idelkadi, 2 years ago

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