Context Navigation

radiation_ifs_rrtm.F90 @ 4773

Last change on this file since 4773 was 4773, checked in by idelkadi, 5 months ago

Update of Ecrad in LMDZ The same organization of the Ecrad offline version is retained in order to facilitate the updating of Ecrad in LMDZ and the comparison between online and offline results. version 1.6.1 of Ecrad (https://github.com/lguez/ecrad.git)

Implementation of the double call of Ecrad in LMDZ

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

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

Download in other formats:

Original Format