Context Navigation

radiation_ifs_rrtm.F90 @ 4013

Last change on this file since 4013 was 3908, checked in by idelkadi, 4 years ago

Online implementation of the radiative transfer code ECRAD in the LMDZ model.

Inclusion of the ecrad directory containing the sources of the ECRAD code
- interface routine : radiation_scheme.F90
Adaptation of compilation scripts :
- compilation under CPP key CPP_ECRAD
- compilation with option "-rad ecard" or "-ecard true"
- The "-rad old/rtm/ecran" build option will need to replace the "-rrtm true" and "-ecrad true" options in the future.
Runing LMDZ simulations with ecrad, you need :
- logical key iflag_rrtm = 2 in physiq.def
- namelist_ecrad (DefLists?)
- the directory "data" containing the configuration files is temporarily placed in ../libfphylmd/ecrad/
Compilation and execution are tested in the 1D case. The repository under svn would allow to continue the implementation work: tests, verification of the results, ...

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

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

Download in other formats:

Original Format