Context Navigation

radiation_ifs_rrtm.F90 @ 4853

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

Ecrad update in LMDZ / Implementation of Ecrad double call in LMDZ

version 1.6.1 (commit 210d7911380f53a788c3cad73b3cf9b4e022ef87)
interface routines between lmdz and ecrad grouped in a new "lmdz" directory
double call method redesigned so as to go through the Ecrad initialization and configuration part only once for the entire simulation
clean-up in the read.F routine: delete unitules arguments
modification of default gas model in namelist (default: ECCKD)

File size: 37.3 KB

Line
1	! radiation_ifs_rrtm.F90 - Interface to IFS implementation of RRTM-G
2	!
3	! (C) Copyright 2015- ECMWF.
4	!
5	! This software is licensed under the terms of the Apache Licence Version 2.0
6	! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
7	!
8	! In applying this licence, ECMWF does not waive the privileges and immunities
9	! granted to it by virtue of its status as an intergovernmental organisation
10	! nor does it submit to any jurisdiction.
11	!
12	! Author: Robin Hogan
13	! Email: r.j.hogan@ecmwf.int
14	!
15	! Modifications
16	! 2017-04-11 R. Hogan Receive "surface" dummy argument
17	! 2017-09-08 R. Hogan Reverted some changes
18	! 2017-10-18 R. Hogan Added planck_function public function
19	! 2018-01-11 R. Hogan Added optional spectral scaling of incoming solar radiation
20	! 2018-02-22 R. Hogan Optimized reverse indexing of heights
21	! 2018-05-05 R. Hogan gas_optics can be called for reduced number of levels
22	! 2019-01-02 R. Hogan Initialize shortwave props to zero in case sun below horizon
23
24	module radiation_ifs_rrtm
25
26	implicit none
27
28	public :: setup_gas_optics, gas_optics, planck_function, set_gas_units
29
30	contains
31
32	!---------------------------------------------------------------------
33	! Setup the IFS implementation of RRTM-G gas absorption model
34	subroutine setup_gas_optics(config, directory)
35
36	use yoerrtm, only : jpglw
37	use yoesrtm, only : jpgsw
38	use yoerrtftr, only : ngb_lw => ngb
39	use yoesrtm, only : ngb_sw => ngbsw
40	use yomhook, only : lhook, dr_hook, jphook
41
42	use radiation_config
43	use radiation_spectral_definition, only &
44	& : SolarReferenceTemperature, TerrestrialReferenceTemperature
45
46	type(config_type), intent(inout), target :: config
47	character(len=*), intent(in) :: directory
48
49	integer :: irep ! For implied do
50
51	integer, parameter :: RRTM_GPOINT_REORDERING_LW(140) = (/ &
52	& 89, 90, 139, 77, 137, 69, 131, 97, 91, 70, 78, 71, 53, 72, 123, 54, 79, 98, &
53	& 92, 55, 80, 132, 124, 81, 73, 56, 99, 82, 57, 23, 125, 100, 24, 74, 93, 58, 25, &
54	& 83, 126, 75, 26, 11, 101, 133, 59, 27, 76, 140, 12, 84, 102, 94, 28, 127, 85, &
55	& 13, 39, 60, 86, 103, 87, 109, 14, 29, 115, 40, 95, 15, 61, 88, 41, 110, 104, 1, &
56	& 116, 42, 30, 134, 128, 138, 96, 62, 16, 43, 117, 63, 111, 44, 2, 64, 31, 65, &
57	& 105, 17, 45, 66, 118, 32, 3, 33, 67, 18, 129, 135, 46, 112, 34, 106, 68, 35, 4, &
58	& 119, 36, 47, 107, 19, 37, 38, 113, 48, 130, 5, 120, 49, 108, 20, 50, 51, 114, &
59	& 21, 121, 52, 136, 122, 6, 22, 7, 8, 9, 10 &
60	& /)
61	integer, parameter :: RRTM_GPOINT_REORDERING_SW(112) = (/ &
62	& 35, 45, 19, 27, 36, 57, 20, 46, 58, 21, 28, 67, 55, 68, 37, 1, 69, 22, 29, 59, &
63	& 78, 101, 79, 77, 70, 76, 47, 75, 30, 81, 60, 102, 80, 82, 23, 2, 83, 84, 85, &
64	& 86, 103, 61, 31, 87, 56, 38, 71, 48, 88, 3, 62, 89, 24, 7, 49, 32, 104, 72, 90, &
65	& 63, 39, 4, 8, 50, 91, 64, 40, 33, 25, 51, 95, 96, 73, 65, 9, 41, 97, 92, 105, &
66	& 52, 5, 98, 10, 42, 99, 100, 66, 11, 74, 34, 53, 26, 6, 106, 12, 43, 13, 54, 93, &
67	& 44, 107, 94, 14, 108, 15, 16, 109, 17, 18, 110, 111, 112 &
68	& /)
69
70	logical :: do_sw, do_lw
71
72	real(jphook) :: hook_handle
73
74	!#include "surdi.intfb.h"
75	#include "surrtab.intfb.h"
76	#include "surrtpk.intfb.h"
77	#include "surrtrf.intfb.h"
78	#include "rrtm_init_140gp.intfb.h"
79	#include "srtm_init.intfb.h"
80
81	if (lhook) call dr_hook('radiation_ifs_rrtm:setup_gas_optics',0,hook_handle)
82
83	do_sw = (config%do_sw .and. config%i_gas_model_sw == IGasModelIFSRRTMG)
84	do_lw = (config%do_lw .and. config%i_gas_model_lw == IGasModelIFSRRTMG)
85
86	! The IFS implementation of RRTMG uses many global variables. In
87	! the IFS these will have been set up already; otherwise set them
88	! up now.
89	if (config%do_setup_ifsrrtm) then
90	!call SURDI
91	call SURRTAB
92	call SURRTPK
93	call SURRTRF
94	if (do_lw) then
95	call RRTM_INIT_140GP(directory)
96	end if
97	if (do_sw) then
98	call SRTM_INIT(directory)
99	end if
100	end if
101
102	if (do_sw) then
103
104	! Cloud and aerosol properties can only be defined per band
105	config%do_cloud_aerosol_per_sw_g_point = .false.
106	config%n_g_sw = jpgsw
107	config%n_bands_sw = 14
108	! Wavenumber ranges of each band may be needed so that the user
109	! can compute UV and photosynthetically active radiation for a
110	! particular wavelength range
111	call config%gas_optics_sw%spectral_def%allocate_bands_only(SolarReferenceTemperature, &
112	& [2600.0_jprb, 3250.0_jprb, 4000.0_jprb, 4650.0_jprb, 5150.0_jprb, 6150.0_jprb, 7700.0_jprb, &
113	& 8050.0_jprb, 12850.0_jprb, 16000.0_jprb, 22650.0_jprb, 29000.0_jprb, 38000.0_jprb, 820.0_jprb], &
114	& [3250.0_jprb, 4000.0_jprb, 4650.0_jprb, 5150.0_jprb, 6150.0_jprb, 7700.0_jprb, 8050.0_jprb, &
115	& 12850.0_jprb, 16000.0_jprb, 22650.0_jprb, 29000.0_jprb, 38000.0_jprb, 50000.0_jprb, 2600.0_jprb])
116	allocate(config%i_band_from_g_sw (config%n_g_sw))
117	allocate(config%i_band_from_reordered_g_sw(config%n_g_sw))
118	allocate(config%i_g_from_reordered_g_sw (config%n_g_sw))
119	! Shortwave starts at 16: need to start at 1
120	config%i_band_from_g_sw = ngb_sw - ngb_sw(1)+1
121
122	if (config%i_solver_sw == ISolverSpartacus) then
123	! SPARTACUS requires g points ordered in approximately
124	! increasing order of optical depth
125	config%i_g_from_reordered_g_sw = RRTM_GPOINT_REORDERING_SW
126	else
127	! Implied-do for no reordering
128	! config%i_g_from_reordered_g_sw = RRTM_GPOINT_REORDERING_SW
129	config%i_g_from_reordered_g_sw = (/ (irep, irep=1,config%n_g_sw) /)
130	end if
131
132	config%i_band_from_reordered_g_sw &
133	= config%i_band_from_g_sw(config%i_g_from_reordered_g_sw)
134
135	! The i_spec_* variables are used solely for storing spectral
136	! data, and this can either be by band or by g-point
137	if (config%do_save_spectral_flux .or. config%do_toa_spectral_flux) then
138	if (config%do_save_gpoint_flux) then
139	config%n_spec_sw = config%n_g_sw
140	config%i_spec_from_reordered_g_sw => config%i_g_from_reordered_g_sw
141	else
142	config%n_spec_sw = config%n_bands_sw
143	config%i_spec_from_reordered_g_sw => config%i_band_from_reordered_g_sw
144	end if
145	else
146	config%n_spec_sw = 0
147	nullify(config%i_spec_from_reordered_g_sw)
148	end if
149
150	end if
151
152	if (do_lw) then
153	! Cloud and aerosol properties can only be defined per band
154	config%do_cloud_aerosol_per_lw_g_point = .false.
155	config%n_g_lw = jpglw
156	config%n_bands_lw = 16
157	call config%gas_optics_lw%spectral_def%allocate_bands_only(TerrestrialReferenceTemperature, &
158	& [10.0_jprb, 350.0_jprb, 500.0_jprb, 630.0_jprb, 700.0_jprb, 820.0_jprb, 980.0_jprb, 1080.0_jprb, &
159	& 1180.0_jprb, 1390.0_jprb, 1480.0_jprb, 1800.0_jprb, 2080.0_jprb, 2250.0_jprb, 2380.0_jprb, 2600.0_jprb], &
160	& [350.0_jprb, 500.0_jprb, 630.0_jprb, 700.0_jprb, 820.0_jprb, 980.0_jprb, 1080.0_jprb, 1180.0_jprb, &
161	& 1390.0_jprb, 1480.0_jprb, 1800.0_jprb, 2080.0_jprb, 2250.0_jprb, 2380.0_jprb, 2600.0_jprb, 3250.0_jprb])
162	allocate(config%i_band_from_g_lw (config%n_g_lw))
163	allocate(config%i_band_from_reordered_g_lw(config%n_g_lw))
164	allocate(config%i_g_from_reordered_g_lw (config%n_g_lw))
165	config%i_band_from_g_lw = ngb_lw
166
167	if (config%i_solver_lw == ISolverSpartacus) then
168	! SPARTACUS requires g points ordered in approximately
169	! increasing order of optical depth
170	config%i_g_from_reordered_g_lw = RRTM_GPOINT_REORDERING_LW
171	else
172	! Implied-do for no reordering
173	config%i_g_from_reordered_g_lw = (/ (irep, irep=1,config%n_g_lw) /)
174	end if
175
176	config%i_band_from_reordered_g_lw &
177	= config%i_band_from_g_lw(config%i_g_from_reordered_g_lw)
178
179	! The i_spec_* variables are used solely for storing spectral
180	! data, and this can either be by band or by g-point
181	if (config%do_save_spectral_flux .or. config%do_toa_spectral_flux) then
182	if (config%do_save_gpoint_flux) then
183	config%n_spec_lw = config%n_g_lw
184	config%i_spec_from_reordered_g_lw => config%i_g_from_reordered_g_lw
185	else
186	config%n_spec_lw = config%n_bands_lw
187	config%i_spec_from_reordered_g_lw => config%i_band_from_reordered_g_lw
188	end if
189	else
190	config%n_spec_lw = 0
191	nullify(config%i_spec_from_reordered_g_lw)
192	end if
193
194	end if
195
196	if (lhook) call dr_hook('radiation_ifs_rrtm:setup_gas_optics',1,hook_handle)
197
198	end subroutine setup_gas_optics
199
200
201	!---------------------------------------------------------------------
202	! Scale gas mixing ratios according to required units
203	subroutine set_gas_units(gas)
204
205	use radiation_gas, only : gas_type, IMassMixingRatio
206	type(gas_type), intent(inout) :: gas
207
208	call gas%set_units(IMassMixingRatio)
209
210	end subroutine set_gas_units
211
212
213	!---------------------------------------------------------------------
214	! Compute gas optical depths, shortwave scattering, Planck function
215	! and incoming shortwave radiation at top-of-atmosphere
216	subroutine gas_optics(ncol,nlev,istartcol,iendcol, &
217	& config, single_level, thermodynamics, gas, &
218	& od_lw, od_sw, ssa_sw, lw_albedo, planck_hl, lw_emission, &
219	& incoming_sw)
220
221	use parkind1, only : jprb, jpim
222
223	USE PARRRTM , ONLY : JPBAND, JPXSEC, JPINPX
224	USE YOERRTM , ONLY : JPGPT_LW => JPGPT
225	USE YOESRTM , ONLY : JPGPT_SW => JPGPT
226	use yomhook , only : lhook, dr_hook, jphook
227
228	use radiation_config, only : config_type, ISolverSpartacus, IGasModelIFSRRTMG
229	use radiation_thermodynamics, only : thermodynamics_type
230	use radiation_single_level, only : single_level_type
231	use radiation_gas
232
233	integer, intent(in) :: ncol ! number of columns
234	integer, intent(in) :: nlev ! number of model levels
235	integer, intent(in) :: istartcol, iendcol ! range of columns to process
236	type(config_type), intent(in) :: config
237	type(single_level_type), intent(in) :: single_level
238	type(thermodynamics_type),intent(in) :: thermodynamics
239	type(gas_type), intent(in) :: gas
240
241	! Longwave albedo of the surface
242	real(jprb), dimension(config%n_g_lw,istartcol:iendcol), &
243	& intent(in), optional :: lw_albedo
244
245	! Gaseous layer optical depth in longwave and shortwave, and
246	! shortwave single scattering albedo (i.e. fraction of extinction
247	! due to Rayleigh scattering) at each g-point
248	real(jprb), dimension(config%n_g_lw,nlev,istartcol:iendcol), intent(out) :: &
249	& od_lw
250	real(jprb), dimension(config%n_g_sw,nlev,istartcol:iendcol), intent(out) :: &
251	& od_sw, ssa_sw
252
253	! The Planck function (emitted flux from a black body) at half
254	! levels at each longwave g-point
255	real(jprb), dimension(config%n_g_lw,nlev+1,istartcol:iendcol), &
256	& intent(out), optional :: planck_hl
257	! Planck function for the surface (W m-2)
258	real(jprb), dimension(config%n_g_lw,istartcol:iendcol), &
259	& intent(out), optional :: lw_emission
260
261	! The incoming shortwave flux into a plane perpendicular to the
262	! incoming radiation at top-of-atmosphere in each of the shortwave
263	! g-points
264	real(jprb), dimension(config%n_g_sw,istartcol:iendcol), &
265	& intent(out), optional :: incoming_sw
266
267	real(jprb) :: incoming_sw_scale(istartcol:iendcol)
268
269	! The variables in capitals are used in the same way as the
270	! equivalent routine in the IFS
271
272	real(jprb) :: ZOD_LW(JPGPT_LW,nlev,istartcol:iendcol) ! Note ordering of dimensions
273	real(jprb) :: ZOD_SW(istartcol:iendcol,nlev,JPGPT_SW)
274	real(jprb) :: ZSSA_SW(istartcol:iendcol,nlev,JPGPT_SW)
275	real(jprb) :: ZINCSOL(istartcol:iendcol,JPGPT_SW)
276
277	real(jprb) :: ZCOLMOL(istartcol:iendcol,nlev)
278	real(jprb) :: ZCOLDRY(istartcol:iendcol,nlev)
279	real(jprb) :: ZWBRODL(istartcol:iendcol,nlev) !BROADENING GASES,column density (mol/cm2)
280	real(jprb) :: ZCOLBRD(istartcol:iendcol,nlev) !BROADENING GASES, column amount
281	real(jprb) :: ZWKL(istartcol:iendcol,JPINPX,nlev)
282
283	real(jprb) :: ZWX(istartcol:iendcol,JPXSEC,nlev) ! Amount of trace gases
284
285	real(jprb) :: ZFLUXFAC, ZPI
286
287	! - from AER
288	real(jprb) :: ZTAUAERL(istartcol:iendcol,nlev,JPBAND)
289
290	!- from INTFAC
291	real(jprb) :: ZFAC00(istartcol:iendcol,nlev)
292	real(jprb) :: ZFAC01(istartcol:iendcol,nlev)
293	real(jprb) :: ZFAC10(istartcol:iendcol,nlev)
294	real(jprb) :: ZFAC11(istartcol:iendcol,nlev)
295
296	!- from FOR
297	real(jprb) :: ZFORFAC(istartcol:iendcol,nlev)
298	real(jprb) :: ZFORFRAC(istartcol:iendcol,nlev)
299	integer :: INDFOR(istartcol:iendcol,nlev)
300
301	!- from MINOR
302	integer :: INDMINOR(istartcol:iendcol,nlev)
303	real(jprb) :: ZSCALEMINOR(istartcol:iendcol,nlev)
304	real(jprb) :: ZSCALEMINORN2(istartcol:iendcol,nlev)
305	real(jprb) :: ZMINORFRAC(istartcol:iendcol,nlev)
306
307	real(jprb) :: &
308	& ZRAT_H2OCO2(istartcol:iendcol,nlev),ZRAT_H2OCO2_1(istartcol:iendcol,nlev), &
309	& ZRAT_H2OO3(istartcol:iendcol,nlev) ,ZRAT_H2OO3_1(istartcol:iendcol,nlev), &
310	& ZRAT_H2ON2O(istartcol:iendcol,nlev),ZRAT_H2ON2O_1(istartcol:iendcol,nlev), &
311	& ZRAT_H2OCH4(istartcol:iendcol,nlev),ZRAT_H2OCH4_1(istartcol:iendcol,nlev), &
312	& ZRAT_N2OCO2(istartcol:iendcol,nlev),ZRAT_N2OCO2_1(istartcol:iendcol,nlev), &
313	& ZRAT_O3CO2(istartcol:iendcol,nlev) ,ZRAT_O3CO2_1(istartcol:iendcol,nlev)
314
315	!- from INTIND
316	integer :: JP(istartcol:iendcol,nlev)
317	integer :: JT(istartcol:iendcol,nlev)
318	integer :: JT1(istartcol:iendcol,nlev)
319
320	!- from PRECISE
321	real(jprb) :: ZONEMINUS, ZONEMINUS_ARRAY(istartcol:iendcol)
322
323	!- from PROFDATA
324	real(jprb) :: ZCOLH2O(istartcol:iendcol,nlev)
325	real(jprb) :: ZCOLCO2(istartcol:iendcol,nlev)
326	real(jprb) :: ZCOLO3(istartcol:iendcol,nlev)
327	real(jprb) :: ZCOLN2O(istartcol:iendcol,nlev)
328	real(jprb) :: ZCOLCH4(istartcol:iendcol,nlev)
329	real(jprb) :: ZCOLO2(istartcol:iendcol,nlev)
330	real(jprb) :: ZCO2MULT(istartcol:iendcol,nlev)
331	integer :: ILAYTROP(istartcol:iendcol)
332	integer :: ILAYSWTCH(istartcol:iendcol)
333	integer :: ILAYLOW(istartcol:iendcol)
334
335	!- from PROFILE
336	real(jprb) :: ZPAVEL(istartcol:iendcol,nlev)
337	real(jprb) :: ZTAVEL(istartcol:iendcol,nlev)
338	real(jprb) :: ZPZ(istartcol:iendcol,0:nlev)
339	real(jprb) :: ZTZ(istartcol:iendcol,0:nlev)
340
341	!- from SELF
342	real(jprb) :: ZSELFFAC(istartcol:iendcol,nlev)
343	real(jprb) :: ZSELFFRAC(istartcol:iendcol,nlev)
344	integer :: INDSELF(istartcol:iendcol,nlev)
345
346	!- from SP
347	real(jprb) :: ZPFRAC(istartcol:iendcol,JPGPT_LW,nlev)
348
349	!- from SURFACE
350	integer :: IREFLECT(istartcol:iendcol)
351
352	real(jprb) :: pressure_fl(ncol, nlev), temperature_fl(ncol, nlev)
353
354	! If nlev is less than the number of heights at which gas mixing
355	! ratios are stored, then we assume that the lower part of the
356	! atmosphere is required. This enables nlev=1 to be passed in to
357	! the routine, in which case the gas properties of the lowest
358	! layer are provided, useful for canopy radiative transfer.
359	integer :: istartlev, iendlev
360
361	logical :: do_sw, do_lw
362
363	integer :: jlev, jgreorder, jg, ig, iband, jcol
364
365	real(jphook) :: hook_handle
366
367	#include "rrtm_prepare_gases.intfb.h"
368	#include "rrtm_setcoef_140gp.intfb.h"
369	#include "rrtm_gas_optical_depth.intfb.h"
370	#include "srtm_setcoef.intfb.h"
371	#include "srtm_gas_optical_depth.intfb.h"
372
373	if (lhook) call dr_hook('radiation_ifs_rrtm:gas_optics',0,hook_handle)
374
375	do_sw = (config%do_sw .and. config%i_gas_model_sw == IGasModelIFSRRTMG)
376	do_lw = (config%do_lw .and. config%i_gas_model_lw == IGasModelIFSRRTMG)
377
378	! Compute start and end levels for indexing the gas mixing ratio
379	! and thermodynamics arrays
380	iendlev = ubound(gas%mixing_ratio,2)
381	istartlev = iendlev - nlev + 1
382
383	ZPI = 2.0_jprb*ASIN(1.0_jprb)
384	ZFLUXFAC = ZPI * 1.E+4
385	ZONEMINUS = 1.0_jprb - 1.0e-6_jprb
386	ZONEMINUS_ARRAY = ZONEMINUS
387
388	do jlev=1,nlev
389	do jcol= istartcol,iendcol
390	pressure_fl(jcol,jlev) &
391	& = 0.5_jprb * (thermodynamics%pressure_hl(jcol,jlev+istartlev-1) &
392	& +thermodynamics%pressure_hl(jcol,jlev+istartlev))
393	temperature_fl(jcol,jlev) &
394	& = 0.5_jprb * (thermodynamics%temperature_hl(jcol,jlev+istartlev-1) &
395	& +thermodynamics%temperature_hl(jcol,jlev+istartlev))
396	end do
397	end do
398
399	! Check we have gas mixing ratios in the right units
400	call gas%assert_units(IMassMixingRatio)
401
402	! Warning: O2 is hard-coded within the following function so the
403	! user-provided concentrations of this gas are ignored for both
404	! the longwave and shortwave
405	CALL RRTM_PREPARE_GASES &
406	& ( istartcol, iendcol, ncol, nlev, &
407	& thermodynamics%pressure_hl(:,istartlev:iendlev+1), &
408	& pressure_fl, &
409	& thermodynamics%temperature_hl(:,istartlev:iendlev+1), &
410	& temperature_fl, &
411	& gas%mixing_ratio(:,istartlev:iendlev,IH2O), &
412	& gas%mixing_ratio(:,istartlev:iendlev,ICO2), &
413	& gas%mixing_ratio(:,istartlev:iendlev,ICH4), &
414	& gas%mixing_ratio(:,istartlev:iendlev,IN2O), &
415	& gas%mixing_ratio(:,istartlev:iendlev,INO2), &
416	& gas%mixing_ratio(:,istartlev:iendlev,ICFC11), &
417	& gas%mixing_ratio(:,istartlev:iendlev,ICFC12), &
418	& gas%mixing_ratio(:,istartlev:iendlev,IHCFC22), &
419	& gas%mixing_ratio(:,istartlev:iendlev,ICCl4), &
420	& gas%mixing_ratio(:,istartlev:iendlev,IO3), &
421	& ZCOLDRY, ZWBRODL,ZWKL, ZWX, &
422	& ZPAVEL , ZTAVEL , ZPZ , ZTZ, IREFLECT)
423
424	if (do_lw) then
425
426	CALL RRTM_SETCOEF_140GP &
427	&( istartcol, iendcol, nlev , ZCOLDRY , ZWBRODL , ZWKL , &
428	& ZFAC00 , ZFAC01 , ZFAC10 , ZFAC11 , ZFORFAC,ZFORFRAC,INDFOR, JP, JT, JT1 , &
429	& ZCOLH2O, ZCOLCO2 , ZCOLO3 , ZCOLN2O, ZCOLCH4, ZCOLO2,ZCO2MULT , ZCOLBRD, &
430	& ILAYTROP,ILAYSWTCH, ILAYLOW, ZPAVEL , ZTAVEL , ZSELFFAC, ZSELFFRAC, INDSELF, &
431	& INDMINOR,ZSCALEMINOR,ZSCALEMINORN2,ZMINORFRAC,&
432	& ZRAT_H2OCO2, ZRAT_H2OCO2_1, ZRAT_H2OO3, ZRAT_H2OO3_1, &
433	& ZRAT_H2ON2O, ZRAT_H2ON2O_1, ZRAT_H2OCH4, ZRAT_H2OCH4_1, &
434	& ZRAT_N2OCO2, ZRAT_N2OCO2_1, ZRAT_O3CO2, ZRAT_O3CO2_1)
435
436	ZTAUAERL(istartcol:iendcol,:,:) = 0.0_jprb
437
438	CALL RRTM_GAS_OPTICAL_DEPTH &
439	&( istartcol, iendcol, nlev, ZOD_LW, ZPAVEL, ZCOLDRY, ZCOLBRD, ZWX ,&
440	& ZTAUAERL, ZFAC00 , ZFAC01, ZFAC10 , ZFAC11 , ZFORFAC,ZFORFRAC,INDFOR, &
441	& JP, JT, JT1, ZONEMINUS ,&
442	& ZCOLH2O , ZCOLCO2, ZCOLO3, ZCOLN2O, ZCOLCH4, ZCOLO2,ZCO2MULT ,&
443	& ILAYTROP, ILAYSWTCH,ILAYLOW, ZSELFFAC, ZSELFFRAC, INDSELF, ZPFRAC, &
444	& INDMINOR,ZSCALEMINOR,ZSCALEMINORN2,ZMINORFRAC,&
445	& ZRAT_H2OCO2, ZRAT_H2OCO2_1, ZRAT_H2OO3, ZRAT_H2OO3_1, &
446	& ZRAT_H2ON2O, ZRAT_H2ON2O_1, ZRAT_H2OCH4, ZRAT_H2OCH4_1, &
447	& ZRAT_N2OCO2, ZRAT_N2OCO2_1, ZRAT_O3CO2, ZRAT_O3CO2_1)
448
449	if (present(lw_albedo)) then
450
451	call planck_function_atmos(nlev, istartcol, iendcol, config, &
452	& thermodynamics, ZPFRAC, planck_hl)
453
454	if (single_level%is_simple_surface) then
455	call planck_function_surf(istartcol, iendcol, config, &
456	& single_level%skin_temperature, ZPFRAC(:,:,1), &
457	& lw_emission)
458
459	! The following can be used to extract the parameters defined at
460	! the top of the planck_function routine below:
461	!write(*,'(a,140(e12.5,","),a)') 'ZPFRAC_surf=[', &
462	!& sum(ZPFRAC(istartcol:iendcol,:,1),1) / (iendcol+1-istartcol), ']'
463
464	! lw_emission at this point is actually the planck function of
465	! the surface
466	lw_emission = lw_emission * (1.0_jprb - lw_albedo)
467	else
468	! Longwave emission has already been computed
469	if (config%use_canopy_full_spectrum_lw) then
470	lw_emission = transpose(single_level%lw_emission(istartcol:iendcol,:))
471	else
472	lw_emission = transpose(single_level%lw_emission(istartcol:iendcol, &
473	& config%i_emiss_from_band_lw(config%i_band_from_reordered_g_lw)))
474	end if
475	end if
476
477	end if
478
479	if (config%i_solver_lw == ISolverSpartacus) then
480	! We need to rearrange the gas optics info in memory: reordering
481	! the g points in order of approximately increasing optical
482	! depth (for efficient 3D processing on only the regions of the
483	! spectrum that are optically thin for gases) and reorder in
484	! pressure since the the functions above treat pressure
485	! decreasing with increasing index. Note that the output gas
486	! arrays have dimensions in a different order to the inputs,
487	! so there is some inefficiency here.
488	do jgreorder = 1,config%n_g_lw
489	iband = config%i_band_from_reordered_g_lw(jgreorder)
490	ig = config%i_g_from_reordered_g_lw(jgreorder)
491
492	! Top-of-atmosphere half level
493	do jlev = 1,nlev
494	do jcol = istartcol,iendcol
495	! Some g points can return negative optical depths;
496	! specifically original g points 54-56 which causes
497	! unphysical single-scattering albedo when combined with
498	! aerosol
499	od_lw(jgreorder,jlev,jcol) &
500	& = max(config%min_gas_od_lw, ZOD_LW(ig,nlev+1-jlev,jcol))
501	end do
502	end do
503	end do
504	else
505	! G points have not been reordered
506	do jcol = istartcol,iendcol
507	do jlev = 1,nlev
508	! Check for negative optical depth
509	od_lw(:,jlev,jcol) = max(config%min_gas_od_lw, ZOD_LW(:,nlev+1-jlev,jcol))
510	end do
511	end do
512	end if
513
514	end if
515
516	if (do_sw) then
517
518	CALL SRTM_SETCOEF &
519	& ( istartcol, iendcol, nlev,&
520	& ZPAVEL , ZTAVEL,&
521	& ZCOLDRY , ZWKL,&
522	& ILAYTROP,&
523	& ZCOLCH4 , ZCOLCO2 , ZCOLH2O , ZCOLMOL , ZCOLO2 , ZCOLO3,&
524	& ZFORFAC , ZFORFRAC , INDFOR , ZSELFFAC, ZSELFFRAC, INDSELF, &
525	& ZFAC00 , ZFAC01 , ZFAC10 , ZFAC11,&
526	& JP , JT , JT1 , single_level%cos_sza(istartcol:iendcol) &
527	& )
528
529	! SRTM_GAS_OPTICAL_DEPTH will not initialize profiles when the sun
530	! is below the horizon, so we do it here
531	ZOD_SW(istartcol:iendcol,:,:) = 0.0_jprb
532	ZSSA_SW(istartcol:iendcol,:,:) = 0.0_jprb
533	ZINCSOL(istartcol:iendcol,:) = 0.0_jprb
534
535	CALL SRTM_GAS_OPTICAL_DEPTH &
536	&( istartcol, iendcol , nlev , ZONEMINUS_ARRAY,&
537	& single_level%cos_sza(istartcol:iendcol), ILAYTROP,&
538	& ZCOLCH4 , ZCOLCO2 , ZCOLH2O, ZCOLMOL , ZCOLO2 , ZCOLO3,&
539	& ZFORFAC , ZFORFRAC , INDFOR , ZSELFFAC, ZSELFFRAC, INDSELF,&
540	& ZFAC00 , ZFAC01 , ZFAC10 , ZFAC11 ,&
541	& JP , JT , JT1 ,&
542	& ZOD_SW , ZSSA_SW , ZINCSOL )
543
544	! Scale the incoming solar per band, if requested
545	if (config%use_spectral_solar_scaling) then
546	do jg = 1,JPGPT_SW
547	do jcol = istartcol,iendcol
548	ZINCSOL(jcol,jg) = ZINCSOL(jcol,jg) * &
549	& single_level%spectral_solar_scaling(config%i_band_from_reordered_g_sw(jg))
550	end do
551	end do
552	end if
553
554	! Scaling factor to ensure that the total solar irradiance is as
555	! requested. Note that if the sun is below the horizon then
556	! ZINCSOL will be zero.
557	if (present(incoming_sw)) then
558	do jcol = istartcol,iendcol
559	if (single_level%cos_sza(jcol) > 0.0_jprb) then
560	! Added for DWD (2020)
561	!NEC$ nounroll
562	incoming_sw_scale(jcol) = single_level%solar_irradiance / sum(ZINCSOL(jcol,:))
563	else
564	incoming_sw_scale(jcol) = 1.0_jprb
565	end if
566	end do
567	end if
568
569	if (config%i_solver_sw == ISolverSpartacus) then
570	! if (.true.) then
571	! Account for reordered g points
572	do jgreorder = 1,config%n_g_sw
573	ig = config%i_g_from_reordered_g_sw(jgreorder)
574	do jlev = 1,nlev
575	do jcol = istartcol,iendcol
576	! Check for negative optical depth
577	od_sw (jgreorder,nlev+1-jlev,jcol) &
578	& = max(config%min_gas_od_sw, ZOD_SW (jcol,jlev,ig))
579	ssa_sw(jgreorder,nlev+1-jlev,jcol) = ZSSA_SW(jcol,jlev,ig)
580	end do
581	end do
582	if (present(incoming_sw)) then
583	incoming_sw(jgreorder,:) &
584	& = incoming_sw_scale(:) * ZINCSOL(:,ig)
585	end if
586	end do
587	else
588	! G points have not been reordered
589	do jcol = istartcol,iendcol
590	do jlev = 1,nlev
591	do jg = 1,config%n_g_sw
592	! Check for negative optical depth
593	od_sw (jg,nlev+1-jlev,jcol) = max(config%min_gas_od_sw, ZOD_SW(jcol,jlev,jg))
594	ssa_sw(jg,nlev+1-jlev,jcol) = ZSSA_SW(jcol,jlev,jg)
595	end do
596	end do
597	if (present(incoming_sw)) then
598	do jg = 1,config%n_g_sw
599	incoming_sw(jg,jcol) = incoming_sw_scale(jcol) * ZINCSOL(jcol,jg)
600	end do
601	end if
602	end do
603	end if
604
605	end if
606
607	if (lhook) call dr_hook('radiation_ifs_rrtm:gas_optics',1,hook_handle)
608
609	end subroutine gas_optics
610
611
612	!---------------------------------------------------------------------
613	! Compute Planck function of the atmosphere
614	subroutine planck_function_atmos(nlev,istartcol,iendcol, &
615	config, thermodynamics, PFRAC, &
616	planck_hl)
617
618	use parkind1, only : jprb, jpim
619
620	USE YOERRTM , ONLY : JPGPT_LW => JPGPT
621	use yoerrtwn, only : totplnk, delwave
622
623	use yomhook, only : lhook, dr_hook, jphook
624
625	use radiation_config, only : config_type, ISolverSpartacus
626	use radiation_thermodynamics, only : thermodynamics_type
627	!use radiation_gas
628
629	integer, intent(in) :: nlev ! number of model levels
630	integer, intent(in) :: istartcol, iendcol ! range of columns to process
631	type(config_type), intent(in) :: config
632	type(thermodynamics_type),intent(in) :: thermodynamics
633	real(jprb), intent(in) :: PFRAC(istartcol:iendcol,JPGPT_LW,nlev)
634
635	! The Planck function (emitted flux from a black body) at half
636	! levels at each longwave g-point
637	real(jprb), dimension(config%n_g_lw,nlev+1,istartcol:iendcol), intent(out) :: &
638	& planck_hl
639
640	! Planck function values per band
641	real(jprb), dimension(istartcol:iendcol, config%n_bands_lw) :: planck_store
642
643	! Look-up table variables for Planck function
644	real(jprb), dimension(istartcol:iendcol) :: frac
645	integer, dimension(istartcol:iendcol) :: ind
646
647	! Temperature (K) of a half-level
648	real(jprb) :: temperature
649
650	real(jprb) :: factor, planck_tmp(istartcol:iendcol,config%n_g_lw)
651	real(jprb) :: ZFLUXFAC
652
653	integer :: jlev, jgreorder, jg, ig, iband, jband, jcol, ilevoffset
654
655	real(jphook) :: hook_handle
656
657	if (lhook) call dr_hook('radiation_ifs_rrtm:planck_function_atmos',0,hook_handle)
658
659	ZFLUXFAC = 2.0_jprbASIN(1.0_jprb) 1.0e4_jprb
660
661	! nlev may be less than the number of original levels, in which
662	! case we assume that the user wants the lower part of the
663	! atmosphere
664	ilevoffset = ubound(thermodynamics%temperature_hl,2)-nlev-1
665
666	! Work out interpolations: for each half level, the index of the
667	! lowest interpolation bound, and the fraction into interpolation
668	! interval
669	do jlev = 1,nlev+1
670	do jcol = istartcol,iendcol
671	temperature = thermodynamics%temperature_hl(jcol,jlev+ilevoffset)
672	if (temperature < 339.0_jprb .and. temperature >= 160.0_jprb) then
673	! Linear interpolation between -113 and 66 degC
674	ind(jcol) = int(temperature - 159.0_jprb)
675	frac(jcol) = temperature - int(temperature)
676	else if(temperature >= 339.0_jprb) then
677	! Extrapolation above 66 degC
678	ind(jcol) = 180
679	frac(jcol) = temperature - 339.0_jprb
680	else
681	! Cap below -113 degC (to avoid possible negative Planck
682	! function values)
683	ind(jcol) = 1
684	frac(jcol) = 0.0_jprb
685	end if
686	end do
687
688	! Calculate Planck functions per band
689	do jband = 1,config%n_bands_lw
690	factor = zfluxfac * delwave(jband)
691	do jcol = istartcol,iendcol
692	planck_store(jcol,jband) = factor &
693	& * (totplnk(ind(jcol),jband) &
694	& + frac(jcol)*(totplnk(ind(jcol)+1,jband)-totplnk(ind(jcol),jband)))
695	end do
696	end do
697
698	if (config%i_solver_lw == ISolverSpartacus) then
699	! We need to rearrange the gas optics info in memory:
700	! reordering the g points in order of approximately increasing
701	! optical depth (for efficient 3D processing on only the
702	! regions of the spectrum that are optically thin for gases)
703	! and reorder in pressure since the the functions above treat
704	! pressure decreasing with increasing index.
705	if (jlev == 1) then
706	! Top-of-atmosphere half level - note that PFRAC is on model
707	! levels not half levels
708	do jgreorder = 1,config%n_g_lw
709	iband = config%i_band_from_reordered_g_lw(jgreorder)
710	ig = config%i_g_from_reordered_g_lw(jgreorder)
711	planck_hl(jgreorder,1,:) = planck_store(:,iband) &
712	& * PFRAC(:,ig,nlev)
713	end do
714	else
715	do jgreorder = 1,config%n_g_lw
716	iband = config%i_band_from_reordered_g_lw(jgreorder)
717	ig = config%i_g_from_reordered_g_lw(jgreorder)
718	planck_hl(jgreorder,jlev,:) &
719	& = planck_store(:,iband) &
720	& * PFRAC(:,ig,nlev+2-jlev)
721	end do
722	end if
723	else
724	! G points have not been reordered
725	if (jlev == 1) then
726	! Top-of-atmosphere half level - note that PFRAC is on model
727	! levels not half levels
728	do jg = 1,config%n_g_lw
729	iband = config%i_band_from_g_lw(jg)
730	planck_hl(jg,1,:) = planck_store(:,iband) * PFRAC(:,jg,nlev)
731	end do
732	else
733	do jg = 1,config%n_g_lw
734	iband = config%i_band_from_g_lw(jg)
735	planck_tmp(:,jg) = planck_store(:,iband) * PFRAC(:,jg,nlev+2-jlev)
736	end do
737	do jcol = istartcol,iendcol
738	planck_hl(:,jlev,jcol) = planck_tmp(jcol,:)
739	end do
740	end if
741	end if
742	end do
743
744	if (lhook) call dr_hook('radiation_ifs_rrtm:planck_function_atmos',1,hook_handle)
745
746	end subroutine planck_function_atmos
747
748
749	!---------------------------------------------------------------------
750	! Compute Planck function of the surface
751	subroutine planck_function_surf(istartcol, iendcol, config, temperature, PFRAC, &
752	& planck_surf)
753
754	use parkind1, only : jprb, jpim
755
756	USE YOERRTM , ONLY : JPGPT_LW => JPGPT
757	use yoerrtwn, only : totplnk, delwave
758
759	use yomhook, only : lhook, dr_hook, jphook
760
761	use radiation_config, only : config_type, ISolverSpartacus
762	! use radiation_gas
763
764	integer, intent(in) :: istartcol, iendcol ! range of columns to process
765	type(config_type), intent(in) :: config
766	real(jprb), intent(in) :: temperature(:)
767
768	real(jprb), intent(in) :: PFRAC(istartcol:iendcol,JPGPT_LW)
769
770	! Planck function of the surface (W m-2)
771	real(jprb), dimension(config%n_g_lw,istartcol:iendcol), &
772	& intent(out) :: planck_surf
773
774	! Planck function values per band
775	real(jprb), dimension(istartcol:iendcol, config%n_bands_lw) :: planck_store
776
777	! Look-up table variables for Planck function
778	real(jprb), dimension(istartcol:iendcol) :: frac
779	integer, dimension(istartcol:iendcol) :: ind
780
781	! Temperature (K)
782	real(jprb) :: Tsurf
783
784	real(jprb) :: factor
785	real(jprb) :: ZFLUXFAC
786
787	integer :: jgreorder, jg, ig, iband, jband, jcol
788
789	real(jphook) :: hook_handle
790
791	if (lhook) call dr_hook('radiation_ifs_rrtm:planck_function_surf',0,hook_handle)
792
793	ZFLUXFAC = 2.0_jprbASIN(1.0_jprb) 1.0e4_jprb
794
795	! Work out surface interpolations
796	do jcol = istartcol,iendcol
797	Tsurf = temperature(jcol)
798	if (Tsurf < 339.0_jprb .and. Tsurf >= 160.0_jprb) then
799	! Linear interpolation between -113 and 66 degC
800	ind(jcol) = int(Tsurf - 159.0_jprb)
801	frac(jcol) = Tsurf - int(Tsurf)
802	else if(Tsurf >= 339.0_jprb) then
803	! Extrapolation above 66 degC
804	ind(jcol) = 180
805	frac(jcol) = Tsurf - 339.0_jprb
806	else
807	! Cap below -113 degC (to avoid possible negative Planck
808	! function values)
809	ind(jcol) = 1
810	frac(jcol) = 0.0_jprb
811	end if
812	end do
813
814	! Calculate Planck functions per band
815	do jband = 1,config%n_bands_lw
816	factor = zfluxfac * delwave(jband)
817	do jcol = istartcol,iendcol
818	planck_store(jcol,jband) = factor &
819	& * (totplnk(ind(jcol),jband) &
820	& + frac(jcol)*(totplnk(ind(jcol)+1,jband)-totplnk(ind(jcol),jband)))
821	end do
822	end do
823
824	if (config%i_solver_lw == ISolverSpartacus) then
825	! We need to rearrange the gas optics info in memory: reordering
826	! the g points in order of approximately increasing optical
827	! depth (for efficient 3D processing on only the regions of
828	! the spectrum that are optically thin for gases) and reorder
829	! in pressure since the the functions above treat pressure
830	! decreasing with increasing index.
831	do jgreorder = 1,config%n_g_lw
832	iband = config%i_band_from_reordered_g_lw(jgreorder)
833	ig = config%i_g_from_reordered_g_lw(jgreorder)
834	planck_surf(jgreorder,:) = planck_store(:,iband) * PFRAC(:,ig)
835	end do
836	else
837	! G points have not been reordered
838	do jg = 1,config%n_g_lw
839	iband = config%i_band_from_g_lw(jg)
840	planck_surf(jg,:) = planck_store(:,iband) * PFRAC(:,jg)
841	end do
842	end if
843
844	if (lhook) call dr_hook('radiation_ifs_rrtm:planck_function_surf',1,hook_handle)
845
846	end subroutine planck_function_surf
847
848
849	!---------------------------------------------------------------------
850	! Externally facing function for computing the Planck function
851	! without reference to any gas profile; typically this would be used
852	! for computing the emission by facets of a complex surface. Note
853	! that this uses fixed "PFRAC" values, obtained by averaging over
854	! those derived from RRTM-G for near-surface conditions over a line
855	! of meridian from the ECMWF model.
856	subroutine planck_function(config, temperature, planck_surf)
857
858	use parkind1, only : jprb, jpim
859
860	use radiation_config, only : config_type
861
862	type(config_type), intent(in) :: config
863	real(jprb), intent(in) :: temperature
864
865	! Planck function of the surface (W m-2)
866	real(jprb), dimension(config%n_g_lw), &
867	& intent(out) :: planck_surf
868
869	! Fraction of each band contributed by each g-point within
870	! it. Since there are 16 bands, this array sums to 16
871	real(jprb), parameter, dimension(1,140) :: frac &
872	= reshape( (/ 0.21227E+00, 0.18897E+00, 0.25491E+00, 0.17864E+00, 0.11735E+00, 0.38298E-01, 0.57871E-02, &
873	& 0.31753E-02, 0.53169E-03, 0.76476E-04, 0.16388E+00, 0.15241E+00, 0.14290E+00, 0.12864E+00, &
874	& 0.11615E+00, 0.10047E+00, 0.80013E-01, 0.60445E-01, 0.44918E-01, 0.63395E-02, 0.32942E-02, &
875	& 0.54541E-03, 0.15380E+00, 0.15194E+00, 0.14339E+00, 0.13138E+00, 0.11701E+00, 0.10081E+00, &
876	& 0.82296E-01, 0.61735E-01, 0.41918E-01, 0.45918E-02, 0.37743E-02, 0.30121E-02, 0.22500E-02, &
877	& 0.14490E-02, 0.55410E-03, 0.78364E-04, 0.15938E+00, 0.15146E+00, 0.14213E+00, 0.13079E+00, &
878	& 0.11672E+00, 0.10053E+00, 0.81566E-01, 0.61126E-01, 0.41150E-01, 0.44488E-02, 0.36950E-02, &
879	& 0.29101E-02, 0.21357E-02, 0.19609E-02, 0.14134E+00, 0.14390E+00, 0.13913E+00, 0.13246E+00, &
880	& 0.12185E+00, 0.10596E+00, 0.87518E-01, 0.66164E-01, 0.44862E-01, 0.49402E-02, 0.40857E-02, &
881	& 0.32288E-02, 0.23613E-02, 0.15406E-02, 0.58258E-03, 0.82171E-04, 0.29127E+00, 0.28252E+00, &
882	& 0.22590E+00, 0.14314E+00, 0.45494E-01, 0.71792E-02, 0.38483E-02, 0.65712E-03, 0.29810E+00, &
883	& 0.27559E+00, 0.11997E+00, 0.10351E+00, 0.84515E-01, 0.62253E-01, 0.41050E-01, 0.44217E-02, &
884	& 0.36946E-02, 0.29113E-02, 0.34290E-02, 0.55993E-03, 0.31441E+00, 0.27586E+00, 0.21297E+00, &
885	& 0.14064E+00, 0.45588E-01, 0.65665E-02, 0.34232E-02, 0.53199E-03, 0.19811E+00, 0.16833E+00, &
886	& 0.13536E+00, 0.11549E+00, 0.10649E+00, 0.93264E-01, 0.75720E-01, 0.56405E-01, 0.41865E-01, &
887	& 0.59331E-02, 0.26510E-02, 0.40040E-03, 0.32328E+00, 0.26636E+00, 0.21397E+00, 0.14038E+00, &
888	& 0.52142E-01, 0.38852E-02, 0.14601E+00, 0.13824E+00, 0.27703E+00, 0.22388E+00, 0.15446E+00, &
889	& 0.48687E-01, 0.98054E-02, 0.18870E-02, 0.11961E+00, 0.12106E+00, 0.13215E+00, 0.13516E+00, &
890	& 0.25249E+00, 0.16542E+00, 0.68157E-01, 0.59725E-02, 0.49258E+00, 0.33651E+00, 0.16182E+00, &
891	& 0.90984E-02, 0.95202E+00, 0.47978E-01, 0.91716E+00, 0.82857E-01, 0.77464E+00, 0.22536E+00 /), (/ 1,140 /) )
892
893	call planck_function_surf(1, 1, config, spread(temperature,1,1), &
894	& frac, planck_surf)
895
896	end subroutine planck_function
897
898	end module radiation_ifs_rrtm

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format