source: LMDZ6/trunk/libf/phylmd/ecrad/radiation/radiation_flux.F90 @ 4773

Last change on this file since 4773 was 4773, checked in by idelkadi, 5 months ago
  • Update of Ecrad in LMDZ The same organization of the Ecrad offline version is retained in order to facilitate the updating of Ecrad in LMDZ and the comparison between online and offline results. version 1.6.1 of Ecrad (https://github.com/lguez/ecrad.git)
  • Implementation of the double call of Ecrad in LMDZ


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
23module 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
125contains
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 
855end module radiation_flux
Note: See TracBrowser for help on using the repository browser.