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 |
---|