source: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad.v1.5.1/radiation_flux.F90 @ 5218

Last change on this file since 5218 was 5185, checked in by abarral, 3 months ago

Replace REPROBUS CPP KEY by logical using handmade wonky wrapper

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