Context Navigation

radiation_flux.F90 @ 5450

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

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

Download in other formats:

Original Format