Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

radiation_flux.F90 @ 5451

Last change on this file since 5451 was 4728, checked in by idelkadi, 14 months ago

Update of ecrad in the LMDZ_ECRad branch of LMDZ:

version 1.6.1 of ecrad
files are no longer grouped in the same ecrad directory.
the structure of ecrad offline is preserved to facilitate updating in LMDZ
cfg.bld modified to take into account the new added subdirectories.
the interface routines and those added in ecrad are moved to the phylmd directory

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

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

Download in other formats:

Original Format