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

Last change on this file since 3981 was 3946, checked in by idelkadi, 3 years ago

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