Context Navigation

radiation_flux.F90 @ 5422

Last change on this file since 5422 was 4853, checked in by idelkadi, 9 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: 36.1 KB

Line
1	! radiation_flux.F90 - Derived type to store the output fluxes
2	!
3	! (C) Copyright 2014- 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-09-08 R. Hogan Store g-point fluxes
17	! 2017-10-23 R. Hogan Renamed single-character variables
18	! 2019-01-08 R. Hogan Added "indexed_sum_profile"
19	! 2019-01-14 R. Hogan out_of_physical_bounds calls routine in radiation_config
20	! 2021-01-20 R. Hogan Added heating_rate_out_of_physical_bounds function
21	! 2022-12-07 R. Hogan Added top-of-atmosphere spectral output
22
23	#include "ecrad_config.h"
24
25	module radiation_flux
26
27	use parkind1, only : jprb
28
29	implicit none
30	public
31
32	!---------------------------------------------------------------------
33	! This derived type contains the output from the radiation
34	! calculation. Currently this is solely flux profiles, but in
35	! future surface fluxes in each band may be stored in order that the
36	! calling program can compute surface-radiation such as
37	! photosynthetically active radiation and UV index.
38	type flux_type
39	! All the following are broad-band fluxes in W m-2 with
40	! dimensions (ncol,nlev+1). Note that only those fluxes that are
41	! requested will be used, so clear-sky and direct-beam arrays may
42	! not be allocated
43	real(jprb), allocatable, dimension(:,:) :: &
44	& lw_up, lw_dn, & ! Upwelling and downwelling longwave
45	& sw_up, sw_dn, & ! Upwelling and downwelling shortwave
46	& sw_dn_direct, & ! Direct-beam shortwave into a horizontal plane
47	& lw_up_clear, lw_dn_clear, & ! Clear-sky quantities...
48	& sw_up_clear, sw_dn_clear, &
49	& sw_dn_direct_clear
50	! As above but fluxes in each spectral band in W m-2 with
51	! dimensions (nband,ncol,nlev+1). These are only allocated if
52	! config%do_save_spectral_flux==.true.
53	real(jprb), allocatable, dimension(:,:,:) :: &
54	& lw_up_band, lw_dn_band, & ! Upwelling and downwelling longwave
55	& sw_up_band, sw_dn_band, & ! Upwelling and downwelling shortwave
56	& sw_dn_direct_band, & ! Direct-beam shortwave
57	& lw_up_clear_band, lw_dn_clear_band, & ! Clear-sky quantities...
58	& sw_up_clear_band, sw_dn_clear_band, &
59	& sw_dn_direct_clear_band
60	! Surface downwelling quantities at each g point, dimensioned
61	! (ng,ncol), that are always saved by the solver, except for the
62	! clear-sky ones that are only produced if
63	! config%do_clear==.true.
64	real(jprb), allocatable, dimension(:,:) :: &
65	& lw_dn_surf_g, lw_dn_surf_clear_g, &
66	& sw_dn_diffuse_surf_g, sw_dn_direct_surf_g, &
67	& sw_dn_diffuse_surf_clear_g, sw_dn_direct_surf_clear_g
68	! Top-of-atmosphere quantities at each g point, dimensioned
69	! (ng,ncol), that are always saved by the solver, except for the
70	! clear-sky ones that are only produced if
71	! config%do_clear==.true.
72	real(jprb), allocatable, dimension(:,:) :: &
73	& lw_up_toa_g, lw_up_toa_clear_g, &
74	& sw_dn_toa_g, sw_up_toa_g, sw_up_toa_clear_g
75	! Shortwave downwelling spectral fluxes in W m-2 at the surface,
76	! from which quantities such as photosynthetically active and UV
77	! radiation can be computed. Only allocated if
78	! config%do_surface_sw_spectral_flux==.true. Note that the
79	! clear-sky quantities are only computed if
80	! config%do_clear==.true., but direct fluxes are computed whether
81	! or not do_direct==.true.. The dimensions are (nband,ncol).
82	real(jprb), allocatable, dimension(:,:) :: &
83	& sw_dn_surf_band, sw_dn_direct_surf_band, &
84	& sw_dn_surf_clear_band, sw_dn_direct_surf_clear_band
85	! Top-of-atmosphere spectral fluxes in W m-2. Only allocated if
86	! config%do_toa_spectral_flux=.true.. Note that the clear-sky
87	! quantities are only computed if config%do_clear==.true.. The
88	! dimensions are (nband,ncol).
89	real(jprb), allocatable, dimension(:,:) :: &
90	& lw_up_toa_band, lw_up_toa_clear_band, &
91	& sw_dn_toa_band, sw_up_toa_band, sw_up_toa_clear_band
92	! Surface downwelling fluxes in W m-2 at the spectral resolution
93	! needed by any subsequent canopy radiative transfer. If
94	! config%use_canopy_full_spectrum_[sw\|lw] then these will be at
95	! g-point resolution; otherwise they will be at
96	! config%n_albedo_bands and config%n_emiss_bands resolution.
97	real(jprb), allocatable, dimension(:,:) :: &
98	& lw_dn_surf_canopy, &
99	& sw_dn_diffuse_surf_canopy, sw_dn_direct_surf_canopy
100
101	! Diagnosed cloud cover from the short- and long-wave solvers
102	real(jprb), allocatable, dimension(:) :: &
103	& cloud_cover_lw, cloud_cover_sw
104	! Longwave derivatives needed by Hogan and Bozzo (2015) method
105	! for approximate longwave updates in between the full radiation
106	! calls: rate of change of upwelling broad-band flux with respect
107	! to surface value, dimensioned (ncol,nlev+1)
108	real(jprb), allocatable, dimension(:,:) :: &
109	& lw_derivatives
110
111	contains
112	procedure :: allocate => allocate_flux_type
113	procedure :: deallocate => deallocate_flux_type
114	procedure :: calc_surface_spectral
115	procedure :: calc_toa_spectral
116	procedure :: out_of_physical_bounds
117	procedure :: heating_rate_out_of_physical_bounds
118	end type flux_type
119
120	! Added for DWD (2020)
121	#ifdef DWD_VECTOR_OPTIMIZATIONS
122	logical, parameter :: use_indexed_sum_vec = .true.
123	#else
124	logical, parameter :: use_indexed_sum_vec = .false.
125	#endif
126
127	contains
128
129	!---------------------------------------------------------------------
130	! Allocate arrays for flux profiles, using config to define which
131	! fluxes are needed. The arrays are dimensioned for columns between
132	! istartcol, iendcol and levels from 1 to nlev+1
133	subroutine allocate_flux_type(this, config, istartcol, iendcol, nlev)
134
135	use yomhook, only : lhook, dr_hook, jphook
136	use radiation_io, only : nulerr, radiation_abort
137	use radiation_config, only : config_type
138
139	integer, intent(in) :: istartcol, iendcol, nlev
140	class(flux_type), intent(inout) :: this
141	type(config_type), intent(in) :: config
142
143	real(jphook) :: hook_handle
144
145	if (lhook) call dr_hook('radiation_flux:allocate',0,hook_handle)
146
147	! Allocate longwave arrays
148	if (config%do_lw) then
149	allocate(this%lw_up(istartcol:iendcol,nlev+1))
150	allocate(this%lw_dn(istartcol:iendcol,nlev+1))
151	if (config%do_clear) then
152	allocate(this%lw_up_clear(istartcol:iendcol,nlev+1))
153	allocate(this%lw_dn_clear(istartcol:iendcol,nlev+1))
154	end if
155
156	if (config%do_save_spectral_flux) then
157	if (config%n_spec_lw == 0) then
158	write(nulerr,'(a)') '*** Error: number of LW spectral points to save not yet defined ' &
159	& // 'so cannot allocate spectral flux arrays'
160	call radiation_abort()
161	end if
162
163	allocate(this%lw_up_band(config%n_spec_lw,istartcol:iendcol,nlev+1))
164	allocate(this%lw_dn_band(config%n_spec_lw,istartcol:iendcol,nlev+1))
165	if (config%do_clear) then
166	allocate(this%lw_up_clear_band(config%n_spec_lw, &
167	& istartcol:iendcol,nlev+1))
168	allocate(this%lw_dn_clear_band(config%n_spec_lw, &
169	& istartcol:iendcol,nlev+1))
170	end if
171	end if
172
173	if (config%do_lw_derivatives) then
174	allocate(this%lw_derivatives(istartcol:iendcol,nlev+1))
175	end if
176
177	if (config%do_toa_spectral_flux) then
178	if (config%n_bands_lw == 0) then
179	write(nulerr,'(a)') '*** Error: number of LW bands not yet defined ' &
180	& // 'so cannot allocate TOA spectral flux arrays'
181	call radiation_abort()
182	end if
183	allocate(this%lw_up_toa_band(config%n_bands_lw, istartcol:iendcol))
184	if (config%do_clear) then
185	allocate(this%lw_up_toa_clear_band(config%n_bands_lw, istartcol:iendcol))
186	end if
187	end if
188
189	! Allocate g-point downwelling fluxes at surface, and TOA fluxes
190	allocate(this%lw_dn_surf_g(config%n_g_lw,istartcol:iendcol))
191	allocate(this%lw_up_toa_g (config%n_g_lw,istartcol:iendcol))
192	if (config%do_clear) then
193	allocate(this%lw_dn_surf_clear_g(config%n_g_lw,istartcol:iendcol))
194	allocate(this%lw_up_toa_clear_g (config%n_g_lw,istartcol:iendcol))
195	end if
196
197	if (config%do_canopy_fluxes_lw) then
198	! Downward fluxes at top of canopy at the spectral resolution
199	! used in the canopy radiative transfer scheme
200	allocate(this%lw_dn_surf_canopy(config%n_canopy_bands_lw,istartcol:iendcol))
201	end if
202	end if
203
204	! Allocate shortwave arrays
205	if (config%do_sw) then
206	allocate(this%sw_up(istartcol:iendcol,nlev+1))
207	allocate(this%sw_dn(istartcol:iendcol,nlev+1))
208	if (config%do_sw_direct) then
209	allocate(this%sw_dn_direct(istartcol:iendcol,nlev+1))
210	end if
211	if (config%do_clear) then
212	allocate(this%sw_up_clear(istartcol:iendcol,nlev+1))
213	allocate(this%sw_dn_clear(istartcol:iendcol,nlev+1))
214	if (config%do_sw_direct) then
215	allocate(this%sw_dn_direct_clear(istartcol:iendcol,nlev+1))
216	end if
217	end if
218
219	if (config%do_save_spectral_flux) then
220	if (config%n_spec_sw == 0) then
221	write(nulerr,'(a)') '*** Error: number of SW spectral points to save not yet defined ' &
222	& // 'so cannot allocate spectral flux arrays'
223	call radiation_abort()
224	end if
225
226	allocate(this%sw_up_band(config%n_spec_sw,istartcol:iendcol,nlev+1))
227	allocate(this%sw_dn_band(config%n_spec_sw,istartcol:iendcol,nlev+1))
228
229	if (config%do_sw_direct) then
230	allocate(this%sw_dn_direct_band(config%n_spec_sw, &
231	& istartcol:iendcol,nlev+1))
232	end if
233	if (config%do_clear) then
234	allocate(this%sw_up_clear_band(config%n_spec_sw, &
235	& istartcol:iendcol,nlev+1))
236	allocate(this%sw_dn_clear_band(config%n_spec_sw, &
237	& istartcol:iendcol,nlev+1))
238	if (config%do_sw_direct) then
239	allocate(this%sw_dn_direct_clear_band(config%n_spec_sw, &
240	& istartcol:iendcol, nlev+1))
241	end if
242	end if
243	end if
244
245	if (config%do_surface_sw_spectral_flux) then
246	if (config%n_bands_sw == 0) then
247	write(nulerr,'(a)') '*** Error: number of SW bands not yet defined ' &
248	& // 'so cannot allocate TOA spectral flux arrays'
249	call radiation_abort()
250	end if
251	allocate(this%sw_dn_surf_band(config%n_bands_sw,istartcol:iendcol))
252	allocate(this%sw_dn_direct_surf_band(config%n_bands_sw,istartcol:iendcol))
253	if (config%do_clear) then
254	allocate(this%sw_dn_surf_clear_band(config%n_bands_sw, &
255	& istartcol:iendcol))
256	allocate(this%sw_dn_direct_surf_clear_band(config%n_bands_sw, &
257	& istartcol:iendcol))
258	end if
259	end if
260
261	if (config%do_toa_spectral_flux) then
262	if (config%n_bands_sw == 0) then
263	write(nulerr,'(a)') '*** Error: number of SW bands not yet defined ' &
264	& // 'so cannot allocate surface spectral flux arrays'
265	call radiation_abort()
266	end if
267	allocate(this%sw_dn_toa_band(config%n_bands_sw, istartcol:iendcol))
268	allocate(this%sw_up_toa_band(config%n_bands_sw, istartcol:iendcol))
269	if (config%do_clear) then
270	allocate(this%sw_up_toa_clear_band(config%n_bands_sw, istartcol:iendcol))
271	end if
272	end if
273
274	! Allocate g-point downwelling fluxes at surface, and TOA fluxes
275	allocate(this%sw_dn_diffuse_surf_g(config%n_g_sw,istartcol:iendcol))
276	allocate(this%sw_dn_direct_surf_g (config%n_g_sw,istartcol:iendcol))
277	allocate(this%sw_dn_toa_g (config%n_g_sw,istartcol:iendcol))
278	allocate(this%sw_up_toa_g (config%n_g_sw,istartcol:iendcol))
279	if (config%do_clear) then
280	allocate(this%sw_dn_diffuse_surf_clear_g(config%n_g_sw,istartcol:iendcol))
281	allocate(this%sw_dn_direct_surf_clear_g (config%n_g_sw,istartcol:iendcol))
282	allocate(this%sw_up_toa_clear_g (config%n_g_sw,istartcol:iendcol))
283	end if
284
285	if (config%do_canopy_fluxes_sw) then
286	! Downward fluxes at top of canopy at the spectral resolution
287	! used in the canopy radiative transfer scheme
288	allocate(this%sw_dn_diffuse_surf_canopy(config%n_canopy_bands_sw,istartcol:iendcol))
289	allocate(this%sw_dn_direct_surf_canopy (config%n_canopy_bands_sw,istartcol:iendcol))
290	end if
291	end if
292
293	! Allocate cloud cover arrays
294	allocate(this%cloud_cover_lw(istartcol:iendcol))
295	allocate(this%cloud_cover_sw(istartcol:iendcol))
296
297	! Some solvers may not write to cloud cover, so we initialize to
298	! an unphysical value
299	this%cloud_cover_lw = -1.0_jprb
300	this%cloud_cover_sw = -1.0_jprb
301
302	if (lhook) call dr_hook('radiation_flux:allocate',1,hook_handle)
303
304	end subroutine allocate_flux_type
305
306
307	!---------------------------------------------------------------------
308	! Deallocate flux arrays
309	subroutine deallocate_flux_type(this)
310
311	use yomhook, only : lhook, dr_hook, jphook
312
313	class(flux_type), intent(inout) :: this
314	real(jphook) :: hook_handle
315
316	if (lhook) call dr_hook('radiation_flux:deallocate',0,hook_handle)
317
318	if (allocated(this%lw_up)) then
319	deallocate(this%lw_up)
320	if (allocated(this%lw_dn)) deallocate(this%lw_dn)
321	if (allocated(this%lw_up_clear)) deallocate(this%lw_up_clear)
322	if (allocated(this%lw_dn_clear)) deallocate(this%lw_dn_clear)
323	end if
324
325	if (allocated(this%sw_up)) then
326	deallocate(this%sw_up)
327	if (allocated(this%sw_dn)) deallocate(this%sw_dn)
328	if (allocated(this%sw_up_clear)) deallocate(this%sw_up_clear)
329	if (allocated(this%sw_dn_clear)) deallocate(this%sw_dn_clear)
330	if (allocated(this%sw_dn_direct)) deallocate(this%sw_dn_direct)
331	if (allocated(this%sw_dn_direct_clear)) &
332	& deallocate(this%sw_dn_direct_clear)
333	end if
334
335	if (allocated(this%lw_up_band)) then
336	deallocate(this%lw_up_band)
337	if (allocated(this%lw_dn_band)) deallocate(this%lw_dn_band)
338	if (allocated(this%lw_up_clear_band)) deallocate(this%lw_up_clear_band)
339	if (allocated(this%lw_dn_clear_band)) deallocate(this%lw_dn_clear_band)
340	end if
341
342	if (allocated(this%sw_up_band)) then
343	deallocate(this%sw_up_band)
344	if (allocated(this%sw_dn_band)) deallocate(this%sw_dn_band)
345	if (allocated(this%sw_up_clear_band)) deallocate(this%sw_up_clear_band)
346	if (allocated(this%sw_dn_clear_band)) deallocate(this%sw_dn_clear_band)
347	if (allocated(this%sw_dn_direct_band)) deallocate(this%sw_dn_direct_band)
348	if (allocated(this%sw_dn_direct_clear_band)) &
349	& deallocate(this%sw_dn_direct_clear_band)
350	end if
351
352	if (allocated(this%sw_dn_surf_band)) then
353	deallocate(this%sw_dn_surf_band)
354	deallocate(this%sw_dn_direct_surf_band)
355	end if
356	if (allocated(this%sw_dn_surf_clear_band)) then
357	deallocate(this%sw_dn_surf_clear_band)
358	deallocate(this%sw_dn_direct_surf_clear_band)
359	end if
360
361	if (allocated(this%lw_dn_surf_canopy)) deallocate(this%lw_dn_surf_canopy)
362	if (allocated(this%sw_dn_diffuse_surf_canopy)) deallocate(this%sw_dn_diffuse_surf_canopy)
363	if (allocated(this%sw_dn_direct_surf_canopy)) deallocate(this%sw_dn_direct_surf_canopy)
364
365	if (allocated(this%cloud_cover_sw)) then
366	deallocate(this%cloud_cover_sw)
367	end if
368	if (allocated(this%cloud_cover_lw)) then
369	deallocate(this%cloud_cover_lw)
370	end if
371
372	if (allocated(this%lw_derivatives)) then
373	deallocate(this%lw_derivatives)
374	end if
375
376	if (allocated(this%lw_dn_surf_g)) deallocate(this%lw_dn_surf_g)
377	if (allocated(this%lw_dn_surf_clear_g)) deallocate(this%lw_dn_surf_clear_g)
378	if (allocated(this%sw_dn_diffuse_surf_g)) deallocate(this%sw_dn_diffuse_surf_g)
379	if (allocated(this%sw_dn_direct_surf_g)) deallocate(this%sw_dn_direct_surf_g)
380	if (allocated(this%sw_dn_diffuse_surf_clear_g)) deallocate(this%sw_dn_diffuse_surf_clear_g)
381	if (allocated(this%sw_dn_direct_surf_clear_g)) deallocate(this%sw_dn_direct_surf_clear_g)
382
383	if (allocated(this%lw_up_toa_g)) deallocate(this%lw_up_toa_g)
384	if (allocated(this%sw_up_toa_g)) deallocate(this%sw_up_toa_g)
385	if (allocated(this%sw_dn_toa_g)) deallocate(this%sw_dn_toa_g)
386	if (allocated(this%lw_up_toa_clear_g)) deallocate(this%lw_up_toa_clear_g)
387	if (allocated(this%sw_up_toa_clear_g)) deallocate(this%sw_up_toa_clear_g)
388
389	if (lhook) call dr_hook('radiation_flux:deallocate',1,hook_handle)
390
391	end subroutine deallocate_flux_type
392
393
394	!---------------------------------------------------------------------
395	! Calculate surface downwelling fluxes in each band using the
396	! downwelling surface fluxes at each g point
397	subroutine calc_surface_spectral(this, config, istartcol, iendcol)
398
399	use yomhook, only : lhook, dr_hook, jphook
400	use radiation_config, only : config_type
401
402	class(flux_type), intent(inout) :: this
403	type(config_type), intent(in) :: config
404	integer, intent(in) :: istartcol, iendcol
405
406	integer :: jcol, jband, jalbedoband, nalbedoband
407
408	! Longwave surface downwelling in each band needed to compute
409	! canopy fluxes
410	real(jprb) :: lw_dn_surf_band(config%n_bands_lw,istartcol:iendcol)
411
412	real(jphook) :: hook_handle
413
414	if (lhook) call dr_hook('radiation_flux:calc_surface_spectral',0,hook_handle)
415
416	if (config%do_sw .and. config%do_surface_sw_spectral_flux) then
417
418	if (use_indexed_sum_vec) then
419	call indexed_sum_vec(this%sw_dn_direct_surf_g, &
420	& config%i_band_from_reordered_g_sw, &
421	& this%sw_dn_direct_surf_band, istartcol, iendcol)
422	call indexed_sum_vec(this%sw_dn_diffuse_surf_g, &
423	& config%i_band_from_reordered_g_sw, &
424	& this%sw_dn_surf_band, istartcol, iendcol)
425	do jcol = istartcol,iendcol
426	this%sw_dn_surf_band(:,jcol) &
427	& = this%sw_dn_surf_band(:,jcol) &
428	& + this%sw_dn_direct_surf_band(:,jcol)
429	end do
430	else
431	do jcol = istartcol,iendcol
432	call indexed_sum(this%sw_dn_direct_surf_g(:,jcol), &
433	& config%i_band_from_reordered_g_sw, &
434	& this%sw_dn_direct_surf_band(:,jcol))
435	call indexed_sum(this%sw_dn_diffuse_surf_g(:,jcol), &
436	& config%i_band_from_reordered_g_sw, &
437	& this%sw_dn_surf_band(:,jcol))
438	this%sw_dn_surf_band(:,jcol) &
439	& = this%sw_dn_surf_band(:,jcol) &
440	& + this%sw_dn_direct_surf_band(:,jcol)
441	end do
442	end if
443
444	if (config%do_clear) then
445	if (use_indexed_sum_vec) then
446	call indexed_sum_vec(this%sw_dn_direct_surf_clear_g, &
447	& config%i_band_from_reordered_g_sw, &
448	& this%sw_dn_direct_surf_clear_band, istartcol, iendcol)
449	call indexed_sum_vec(this%sw_dn_diffuse_surf_clear_g, &
450	& config%i_band_from_reordered_g_sw, &
451	& this%sw_dn_surf_clear_band, istartcol, iendcol)
452	do jcol = istartcol,iendcol
453	this%sw_dn_surf_clear_band(:,jcol) &
454	& = this%sw_dn_surf_clear_band(:,jcol) &
455	& + this%sw_dn_direct_surf_clear_band(:,jcol)
456	end do
457	else
458	do jcol = istartcol,iendcol
459	call indexed_sum(this%sw_dn_direct_surf_clear_g(:,jcol), &
460	& config%i_band_from_reordered_g_sw, &
461	& this%sw_dn_direct_surf_clear_band(:,jcol))
462	call indexed_sum(this%sw_dn_diffuse_surf_clear_g(:,jcol), &
463	& config%i_band_from_reordered_g_sw, &
464	& this%sw_dn_surf_clear_band(:,jcol))
465	this%sw_dn_surf_clear_band(:,jcol) &
466	& = this%sw_dn_surf_clear_band(:,jcol) &
467	& + this%sw_dn_direct_surf_clear_band(:,jcol)
468	end do
469	end if
470	end if
471
472	end if ! do_surface_sw_spectral_flux
473
474	! Fluxes in bands required for canopy radiative transfer
475	if (config%do_sw .and. config%do_canopy_fluxes_sw) then
476	if (config%use_canopy_full_spectrum_sw) then
477	this%sw_dn_diffuse_surf_canopy(:,istartcol:iendcol) = this%sw_dn_diffuse_surf_g(:,istartcol:iendcol)
478	this%sw_dn_direct_surf_canopy (:,istartcol:iendcol) = this%sw_dn_direct_surf_g (:,istartcol:iendcol)
479	else if (config%do_nearest_spectral_sw_albedo) then
480	if (use_indexed_sum_vec) then
481	call indexed_sum_vec(this%sw_dn_direct_surf_g, &
482	& config%i_albedo_from_band_sw(config%i_band_from_reordered_g_sw), &
483	& this%sw_dn_direct_surf_canopy, istartcol, iendcol)
484	call indexed_sum_vec(this%sw_dn_diffuse_surf_g, &
485	& config%i_albedo_from_band_sw(config%i_band_from_reordered_g_sw), &
486	& this%sw_dn_diffuse_surf_canopy, istartcol, iendcol)
487	else
488	do jcol = istartcol,iendcol
489	call indexed_sum(this%sw_dn_direct_surf_g(:,jcol), &
490	& config%i_albedo_from_band_sw(config%i_band_from_reordered_g_sw), &
491	& this%sw_dn_direct_surf_canopy(:,jcol))
492	call indexed_sum(this%sw_dn_diffuse_surf_g(:,jcol), &
493	& config%i_albedo_from_band_sw(config%i_band_from_reordered_g_sw), &
494	& this%sw_dn_diffuse_surf_canopy(:,jcol))
495	end do
496	end if
497	else
498	! More accurate calculations using weights, but requires
499	! this%sw_dn_[direct_]surf_band to be defined, i.e.
500	! config%do_surface_sw_spectral_flux == .true.
501	nalbedoband = size(config%sw_albedo_weights,1)
502	this%sw_dn_diffuse_surf_canopy(:,istartcol:iendcol) = 0.0_jprb
503	this%sw_dn_direct_surf_canopy (:,istartcol:iendcol) = 0.0_jprb
504	do jband = 1,config%n_bands_sw
505	do jalbedoband = 1,nalbedoband
506	if (config%sw_albedo_weights(jalbedoband,jband) /= 0.0_jprb) then
507	! Initially, "diffuse" is actually "total"
508	this%sw_dn_diffuse_surf_canopy(jalbedoband,istartcol:iendcol) &
509	& = this%sw_dn_diffuse_surf_canopy(jalbedoband,istartcol:iendcol) &
510	& + config%sw_albedo_weights(jalbedoband,jband) &
511	& * this%sw_dn_surf_band(jband,istartcol:iendcol)
512	this%sw_dn_direct_surf_canopy(jalbedoband,istartcol:iendcol) &
513	& = this%sw_dn_direct_surf_canopy(jalbedoband,istartcol:iendcol) &
514	& + config%sw_albedo_weights(jalbedoband,jband) &
515	& * this%sw_dn_direct_surf_band(jband,istartcol:iendcol)
516	end if
517	end do
518	end do
519	! Subtract the direct from total to get diffuse
520	this%sw_dn_diffuse_surf_canopy(:,istartcol:iendcol) &
521	& = this%sw_dn_diffuse_surf_canopy(:,istartcol:iendcol) &
522	& - this%sw_dn_direct_surf_canopy(:,istartcol:iendcol)
523	end if
524
525	end if ! do_canopy_fluxes_sw
526
527	if (config%do_lw .and. config%do_canopy_fluxes_lw) then
528	if (config%use_canopy_full_spectrum_lw) then
529	this%lw_dn_surf_canopy(:,istartcol:iendcol) = this%lw_dn_surf_g(:,istartcol:iendcol)
530	else if (config%do_nearest_spectral_lw_emiss) then
531	if (use_indexed_sum_vec) then
532	call indexed_sum_vec(this%lw_dn_surf_g, &
533	& config%i_emiss_from_band_lw(config%i_band_from_reordered_g_lw), &
534	& this%lw_dn_surf_canopy, istartcol, iendcol)
535	else
536	do jcol = istartcol,iendcol
537	call indexed_sum(this%lw_dn_surf_g(:,jcol), &
538	& config%i_emiss_from_band_lw(config%i_band_from_reordered_g_lw), &
539	& this%lw_dn_surf_canopy(:,jcol))
540	end do
541	end if
542	else
543	! Compute fluxes in each longwave emissivity interval using
544	! weights; first sum over g points to get the values in bands
545	if (use_indexed_sum_vec) then
546	call indexed_sum_vec(this%lw_dn_surf_g, &
547	& config%i_band_from_reordered_g_lw, &
548	& lw_dn_surf_band, istartcol, iendcol)
549	else
550	do jcol = istartcol,iendcol
551	call indexed_sum(this%lw_dn_surf_g(:,jcol), &
552	& config%i_band_from_reordered_g_lw, &
553	& lw_dn_surf_band(:,jcol))
554	end do
555	end if
556	nalbedoband = size(config%lw_emiss_weights,1)
557	this%lw_dn_surf_canopy(:,istartcol:iendcol) = 0.0_jprb
558	do jband = 1,config%n_bands_lw
559	do jalbedoband = 1,nalbedoband
560	if (config%lw_emiss_weights(jalbedoband,jband) /= 0.0_jprb) then
561	this%lw_dn_surf_canopy(jalbedoband,istartcol:iendcol) &
562	& = this%lw_dn_surf_canopy(jalbedoband,istartcol:iendcol) &
563	& + config%lw_emiss_weights(jalbedoband,jband) &
564	& * lw_dn_surf_band(jband,istartcol:iendcol)
565	end if
566	end do
567	end do
568	end if
569	end if
570
571	if (lhook) call dr_hook('radiation_flux:calc_surface_spectral',1,hook_handle)
572
573	end subroutine calc_surface_spectral
574
575
576	!---------------------------------------------------------------------
577	! Calculate top-of-atmosphere fluxes in each band using the fluxes
578	! at each g point
579	subroutine calc_toa_spectral(this, config, istartcol, iendcol)
580
581	use yomhook, only : lhook, dr_hook, jphook
582	use radiation_config, only : config_type
583
584	class(flux_type), intent(inout) :: this
585	type(config_type), intent(in) :: config
586	integer, intent(in) :: istartcol, iendcol
587
588	integer :: jcol, jband
589
590	real(jphook) :: hook_handle
591
592	if (lhook) call dr_hook('radiation_flux:calc_toa_spectral',0,hook_handle)
593
594	if (config%do_sw .and. config%do_toa_spectral_flux) then
595
596	if (use_indexed_sum_vec) then
597	call indexed_sum_vec(this%sw_dn_toa_g, &
598	& config%i_band_from_reordered_g_sw, &
599	& this%sw_dn_toa_band, istartcol, iendcol)
600	call indexed_sum_vec(this%sw_up_toa_g, &
601	& config%i_band_from_reordered_g_sw, &
602	& this%sw_up_toa_band, istartcol, iendcol)
603	else
604	do jcol = istartcol,iendcol
605	call indexed_sum(this%sw_dn_toa_g(:,jcol), &
606	& config%i_band_from_reordered_g_sw, &
607	& this%sw_dn_toa_band(:,jcol))
608	call indexed_sum(this%sw_up_toa_g(:,jcol), &
609	& config%i_band_from_reordered_g_sw, &
610	& this%sw_up_toa_band(:,jcol))
611	end do
612	end if
613
614	if (config%do_clear) then
615	if (use_indexed_sum_vec) then
616	call indexed_sum_vec(this%sw_up_toa_clear_g, &
617	& config%i_band_from_reordered_g_sw, &
618	& this%sw_up_toa_clear_band, istartcol, iendcol)
619	else
620	do jcol = istartcol,iendcol
621	call indexed_sum(this%sw_up_toa_clear_g(:,jcol), &
622	& config%i_band_from_reordered_g_sw, &
623	& this%sw_up_toa_clear_band(:,jcol))
624	end do
625	end if
626	end if
627	end if
628
629	if (config%do_lw .and. config%do_toa_spectral_flux) then
630
631	if (use_indexed_sum_vec) then
632	call indexed_sum_vec(this%lw_up_toa_g, &
633	& config%i_band_from_reordered_g_lw, &
634	& this%lw_up_toa_band, istartcol, iendcol)
635	else
636	do jcol = istartcol,iendcol
637	call indexed_sum(this%lw_up_toa_g(:,jcol), &
638	& config%i_band_from_reordered_g_lw, &
639	& this%lw_up_toa_band(:,jcol))
640	end do
641	end if
642
643	if (config%do_clear) then
644	if (use_indexed_sum_vec) then
645	call indexed_sum_vec(this%lw_up_toa_clear_g, &
646	& config%i_band_from_reordered_g_lw, &
647	& this%lw_up_toa_clear_band, istartcol, iendcol)
648	else
649	do jcol = istartcol,iendcol
650	call indexed_sum(this%lw_up_toa_clear_g(:,jcol), &
651	& config%i_band_from_reordered_g_lw, &
652	& this%lw_up_toa_clear_band(:,jcol))
653	end do
654	end if
655	end if
656	end if
657
658	if (lhook) call dr_hook('radiation_flux:calc_toa_spectral',1,hook_handle)
659
660	end subroutine calc_toa_spectral
661
662
663	!---------------------------------------------------------------------
664	! Return .true. if the most important flux variables are out of a
665	! physically sensible range, optionally only considering columns
666	! between istartcol and iendcol
667	function out_of_physical_bounds(this, istartcol, iendcol) result(is_bad)
668
669	use yomhook, only : lhook, dr_hook, jphook
670	use radiation_check, only : out_of_bounds_2d
671
672	class(flux_type), intent(inout) :: this
673	integer, optional,intent(in) :: istartcol, iendcol
674	logical :: is_bad
675
676	real(jphook) :: hook_handle
677
678	if (lhook) call dr_hook('radiation_flux:out_of_physical_bounds',0,hook_handle)
679
680	is_bad = out_of_bounds_2d(this%lw_up, 'lw_up', 10.0_jprb, 900.0_jprb, .false., i1=istartcol, i2=iendcol) &
681	& .or. out_of_bounds_2d(this%lw_dn, 'lw_dn', 0.0_jprb, 800.0_jprb, .false., i1=istartcol, i2=iendcol) &
682	& .or. out_of_bounds_2d(this%sw_up, 'sw_up', 0.0_jprb, 1500.0_jprb, .false., i1=istartcol, i2=iendcol) &
683	& .or. out_of_bounds_2d(this%sw_dn, 'sw_dn', 0.0_jprb, 1500.0_jprb, .false., i1=istartcol, i2=iendcol) &
684	& .or. out_of_bounds_2d(this%sw_dn_direct, 'sw_dn_direct', 0.0_jprb, 1500.0_jprb, .false., i1=istartcol, i2=iendcol) &
685	& .or. out_of_bounds_2d(this%lw_derivatives, 'lw_derivatives', 0.0_jprb, 1.0_jprb, .false., i1=istartcol, i2=iendcol) &
686	& .or. out_of_bounds_2d(this%sw_dn_surf_band, 'sw_dn_surf_band', 0.0_jprb, 1500.0_jprb, &
687	& .false., j1=istartcol, j2=iendcol) &
688	& .or. out_of_bounds_2d(this%sw_dn_surf_clear_band, 'sw_dn_surf_clear_band', 0.0_jprb, 1500.0_jprb, &
689	& .false., j1=istartcol, j2=iendcol)
690
691	if (lhook) call dr_hook('radiation_flux:out_of_physical_bounds',1,hook_handle)
692
693	end function out_of_physical_bounds
694
695	!---------------------------------------------------------------------
696	! Return .true. if the heating rates are out of a physically
697	! sensible range, optionally only considering columns between
698	! istartcol and iendcol. This function allocates and deallocates
699	! memory due to the requirements for inputs of out_of_bounds_2d.
700	function heating_rate_out_of_physical_bounds(this, nlev, istartcol, iendcol, pressure_hl) result(is_bad)
701
702	use radiation_check, only : out_of_bounds_2d
703	use radiation_constants, only : AccelDueToGravity
704
705	! "Cp" (J kg-1 K-1)
706	real(jprb), parameter :: SpecificHeatDryAir = 1004.0
707
708	class(flux_type), intent(inout) :: this
709	integer, intent(in) :: istartcol, iendcol, nlev
710	logical :: is_bad
711
712	real(jprb), intent(in) :: pressure_hl(:,:)
713
714	real(jprb), allocatable :: hr_K_day(:,:)
715
716	real(jprb) :: scaling(istartcol:iendcol,nlev)
717
718	allocate(hr_K_day(istartcol:iendcol,nlev))
719
720	scaling = -(24.0_jprb * 3600.0_jprb * AccelDueToGravity / SpecificHeatDryAir) &
721	& / (pressure_hl(istartcol:iendcol,2:nlev+1) - pressure_hl(istartcol:iendcol,1:nlev))
722	! Shortwave
723	hr_K_day = scaling * (this%sw_dn(istartcol:iendcol,2:nlev+1) - this%sw_up(istartcol:iendcol,2:nlev+1) &
724	& -this%sw_dn(istartcol:iendcol,1:nlev) + this%sw_up(istartcol:iendcol,1:nlev))
725	is_bad = out_of_bounds_2d(hr_K_day, 'sw_heating_rate_K_day', 0.0_jprb, 200.0_jprb, &
726	& .false., i1=istartcol, i2=iendcol)
727
728	! Longwave
729	hr_K_day = scaling * (this%lw_dn(istartcol:iendcol,2:nlev+1) - this%lw_up(istartcol:iendcol,2:nlev+1) &
730	& -this%lw_dn(istartcol:iendcol,1:nlev) + this%lw_up(istartcol:iendcol,1:nlev))
731	is_bad = is_bad .or. out_of_bounds_2d(hr_K_day, 'lw_heating_rate_K_day', -250.0_jprb, 150.0_jprb, &
732	& .false., i1=istartcol, i2=iendcol)
733
734	deallocate(hr_K_day)
735
736	end function heating_rate_out_of_physical_bounds
737
738
739	!---------------------------------------------------------------------
740	! Sum elements of "source" into "dest" according to index "ind".
741	! "source" and "ind" should have the same size and bounds, and no
742	! element of "ind" should refer outside the bounds of "dest". This
743	! version increments existing contents of "dest".
744	pure subroutine add_indexed_sum(source, ind, dest)
745
746	real(jprb), intent(in) :: source(:)
747	integer, intent(in) :: ind(:)
748	real(jprb), intent(inout) :: dest(:)
749
750	integer :: ig, jg, istart, iend
751
752	istart = lbound(source,1)
753	iend = ubound(source,1)
754
755	do jg = istart, iend
756	ig = ind(jg)
757	dest(ig) = dest(ig) + source(jg)
758	end do
759
760	end subroutine add_indexed_sum
761
762
763	!---------------------------------------------------------------------
764	! As "add_indexed_sum" but this version overwrites existing contents
765	! of "dest"
766	pure subroutine indexed_sum(source, ind, dest)
767
768	real(jprb), intent(in) :: source(:)
769	integer, intent(in) :: ind(:)
770	real(jprb), intent(out) :: dest(:)
771
772	integer :: ig, jg, istart, iend
773
774	dest = 0.0
775
776	istart = lbound(source,1)
777	iend = ubound(source,1)
778
779	do jg = istart, iend
780	ig = ind(jg)
781	dest(ig) = dest(ig) + source(jg)
782	end do
783
784	end subroutine indexed_sum
785
786	!---------------------------------------------------------------------
787	! Vectorized version of "add_indexed_sum"
788	subroutine indexed_sum_vec(source, ind, dest, ist, iend)
789
790	real(jprb), intent(in) :: source(:,:)
791	integer, intent(in) :: ind(:)
792	real(jprb), intent(out) :: dest(:,:)
793	integer, intent(in) :: ist, iend
794
795	integer :: ig, jg, jc
796
797	dest = 0.0
798
799	do jg = lbound(source,1), ubound(source,1)
800	ig = ind(jg)
801	do jc = ist, iend
802	dest(ig,jc) = dest(ig,jc) + source(jg,jc)
803	end do
804	end do
805
806	end subroutine indexed_sum_vec
807
808	!---------------------------------------------------------------------
809	! As "add_indexed_sum" but a whole vertical profiles
810	pure subroutine add_indexed_sum_profile(source, ind, dest)
811
812	real(jprb), intent(in) :: source(:,:)
813	integer, intent(in) :: ind(:)
814	real(jprb), intent(out) :: dest(:,:)
815
816	integer :: ig, jg, istart, iend, jlev, nlev
817
818	istart = lbound(source,1)
819	iend = ubound(source,1)
820	nlev = size(source,2)
821
822	do jlev = 1,nlev
823	do jg = istart, iend
824	ig = ind(jg)
825	dest(ig,jlev) = dest(ig,jlev) + source(jg,jlev)
826	end do
827	end do
828
829	end subroutine add_indexed_sum_profile
830
831
832	!---------------------------------------------------------------------
833	! As "indexed_sum" but a whole vertical profiles
834	pure subroutine indexed_sum_profile(source, ind, dest)
835
836	real(jprb), intent(in) :: source(:,:)
837	integer, intent(in) :: ind(:)
838	real(jprb), intent(out) :: dest(:,:)
839
840	integer :: ig, jg, istart, iend, jlev, nlev
841
842	dest = 0.0
843
844	istart = lbound(source,1)
845	iend = ubound(source,1)
846	nlev = size(source,2)
847
848	do jlev = 1,nlev
849	do jg = istart, iend
850	ig = ind(jg)
851	dest(ig,jlev) = dest(ig,jlev) + source(jg,jlev)
852	end do
853	end do
854
855	end subroutine indexed_sum_profile
856
857	end module radiation_flux

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

Download in other formats:

Original Format