source: LMDZ6/branches/Amaury_dev/libf/phylmd/ecrad.v1.5.1/radiation_spectral_definition.F90 @ 5185

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

Replace REPROBUS CPP KEY by logical using handmade wonky wrapper

File size: 33.7 KB
Line 
1! radiation_spectral_definition.F90 - Derived type to describe a spectral definition
2
3! (C) Copyright 2020- 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! License: see the COPYING file for details
15!
16
17module radiation_spectral_definition
18
19  use parkind1,    only : jprb
20
21  implicit none
22
23  public
24
25  real(jprb), parameter :: SolarReferenceTemperature = 5777.0_jprb ! K
26  real(jprb), parameter :: TerrestrialReferenceTemperature = 273.15_jprb ! K
27
28  !---------------------------------------------------------------------
29  ! A derived type describing the contribution of the g points of a
30  ! correlated k-distribution gas-optics model from each part of the
31  ! spectrum. This is used primarily to map the cloud and aerosol
32  ! optical properties on to the gas g points.
33  type spectral_definition_type
34   
35    ! Spectral mapping of g points
36
37    ! Number of wavenumber intervals
38    integer :: nwav = 0
39    ! Number of k terms / g points
40    integer :: ng   = 0
41    ! Start and end wavenumber (cm-1), dimensioned (nwav)
42    real(jprb), allocatable :: wavenumber1(:)
43    real(jprb), allocatable :: wavenumber2(:)
44    ! Fraction of each g point in each wavenumber interval,
45    ! dimensioned (nwav, ng)
46    real(jprb), allocatable :: gpoint_fraction(:,:)
47
48    ! Band information
49
50    ! Number of bands
51    integer :: nband = 0
52    ! Lower and upper bounds of wavenumber bands (cm-1), dimensioned
53    ! (nband)
54    real(jprb), allocatable :: wavenumber1_band(:)
55    real(jprb), allocatable :: wavenumber2_band(:)
56    ! Band (one based) to which each g point belongs
57    integer,    allocatable :: i_band_number(:)
58
59  contains
60    procedure :: read => read_spectral_definition
61    procedure :: allocate_bands_only
62    procedure :: deallocate
63    procedure :: find => find_wavenumber
64    procedure :: calc_mapping
65    procedure :: calc_mapping_from_bands
66    procedure :: calc_mapping_from_wavenumber_bands
67    procedure :: print_mapping_from_bands
68    procedure :: min_wavenumber, max_wavenumber
69
70  end type spectral_definition_type
71
72contains
73
74  !---------------------------------------------------------------------
75  ! Read the description of a spectral definition from a NetCDF
76  ! file of the type used to describe an ecCKD model
77  subroutine read_spectral_definition(this, file)
78
79    use easy_netcdf, only : netcdf_file
80    use yomhook,     only : lhook, dr_hook
81
82    class(spectral_definition_type), intent(inout) :: this
83    type(netcdf_file),               intent(inout) :: file
84
85    real(jprb) :: hook_handle
86
87    if (lhook) call dr_hook('radiation_spectral_definition:read',0,hook_handle)
88
89    ! Read spectral mapping of g points
90    call file%get('wavenumber1', this%wavenumber1)
91    call file%get('wavenumber2', this%wavenumber2)
92    call file%get('gpoint_fraction', this%gpoint_fraction)
93
94    ! Read band information
95    call file%get('wavenumber1_band', this%wavenumber1_band)
96    call file%get('wavenumber2_band', this%wavenumber2_band)
97    call file%get('band_number', this%i_band_number)
98
99    ! Band number is 0-based: add 1
100    this%i_band_number = this%i_band_number + 1
101
102    this%nwav  = size(this%wavenumber1)
103    this%ng    = size(this%gpoint_fraction, 2);
104    this%nband = size(this%wavenumber1_band)
105
106    if (lhook) call dr_hook('radiation_spectral_definition:read',1,hook_handle)
107
108  end subroutine read_spectral_definition
109
110
111  !---------------------------------------------------------------------
112  ! Store a simple band description by copying over the lower and
113  ! upper wavenumbers of each band
114  subroutine allocate_bands_only(this, wavenumber1, wavenumber2)
115
116    use yomhook,     only : lhook, dr_hook
117
118    class(spectral_definition_type), intent(inout) :: this
119    real(jprb),        dimension(:), intent(in)    :: wavenumber1, wavenumber2
120
121    real(jprb) :: hook_handle
122
123    if (lhook) call dr_hook('radiation_spectral_definition:allocate_bands_only',0,hook_handle)
124
125    call this%deallocate()
126
127    this%nband = size(wavenumber1)
128    allocate(this%wavenumber1_band(this%nband))
129    allocate(this%wavenumber2_band(this%nband))
130    this%wavenumber1_band = wavenumber1
131    this%wavenumber2_band = wavenumber2
132
133    if (lhook) call dr_hook('radiation_spectral_definition:allocate_bands_only',1,hook_handle)
134
135  end subroutine allocate_bands_only
136
137
138  !---------------------------------------------------------------------
139  ! Deallocate memory inside a spectral definition object
140  subroutine deallocate(this)
141
142    class(spectral_definition_type), intent(inout) :: this
143   
144    this%nwav  = 0
145    this%ng    = 0
146    this%nband = 0
147
148    if (allocated(this%wavenumber1))      deallocate(this%wavenumber1)
149    if (allocated(this%wavenumber2))      deallocate(this%wavenumber2)
150    if (allocated(this%wavenumber1_band)) deallocate(this%wavenumber1_band)
151    if (allocated(this%wavenumber2_band)) deallocate(this%wavenumber2_band)
152    if (allocated(this%gpoint_fraction))  deallocate(this%gpoint_fraction)
153    if (allocated(this%i_band_number))    deallocate(this%i_band_number)
154
155  end subroutine deallocate
156
157
158  !---------------------------------------------------------------------
159  ! Find the index to the highest wavenumber in the spectral
160  ! definition that is lower than or equal to "wavenumber", used for
161  ! implementing look-up tables
162  pure function find_wavenumber(this, wavenumber)
163    class(spectral_definition_type), intent(in) :: this
164    real(jprb),                      intent(in) :: wavenumber ! cm-1
165    integer                                     :: find_wavenumber
166
167    if (wavenumber < this%wavenumber1(1) .or. wavenumber > this%wavenumber2(this%nwav)) then
168      ! Wavenumber not present
169      find_wavenumber = 0
170    else
171      find_wavenumber = 1
172      DO while (wavenumber > this%wavenumber2(find_wavenumber) &
173           &    .AND. find_wavenumber < this%nwav)
174        find_wavenumber = find_wavenumber + 1
175      end do
176    end if
177  end function find_wavenumber
178
179
180  !---------------------------------------------------------------------
181  ! Compute a mapping matrix "mapping" that can be used in an
182  ! expression y=matmul(mapping,x) where x is a variable containing
183  ! optical properties at each input "wavenumber", and y is this
184  ! variable mapped on to the spectral intervals in the spectral
185  ! definition "this". Temperature (K) is used to generate a Planck
186  ! function to weight each wavenumber appropriately.
187  subroutine calc_mapping(this, temperature, wavenumber, mapping, use_bands)
188
189    use yomhook,      only : lhook, dr_hook
190    use radiation_io, only : nulerr, radiation_abort
191
192    class(spectral_definition_type), intent(in)    :: this
193    real(jprb),                      intent(in)    :: temperature   ! K
194    real(jprb),                      intent(in)    :: wavenumber(:) ! cm-1
195    real(jprb), allocatable,         intent(inout) :: mapping(:,:)
196    logical,    optional,            intent(in)    :: use_bands
197
198    ! Spectral weights to apply, same length as wavenumber above
199    real(jprb), dimension(:), allocatable :: weight, planck_weight
200
201    ! Wavenumbers (cm-1) marking triangle of influence of a cloud
202    ! spectral point
203    real(jprb) :: wavenum0, wavenum1, wavenum2
204
205    integer    :: nwav ! Number of wavenumbers describing cloud
206
207    ! Indices to wavenumber intervals in spectral definition structure
208    integer    :: isd, isd0, isd1, isd2
209
210    ! Wavenumber index
211    integer    :: iwav
212   
213    ! Loop indices
214    integer    :: jg, jwav, jband
215
216    logical    :: use_bands_local
217
218    real(jprb) :: hook_handle
219
220    if (lhook) call dr_hook('radiation_spectral_definition:calc_mapping',0,hook_handle)
221
222    if (present(use_bands)) then
223      use_bands_local = use_bands
224    else
225      use_bands_local = .false.
226    end if
227
228    nwav = size(wavenumber)
229
230    if (allocated(mapping)) then
231      deallocate(mapping)
232    end if
233   
234    ! Define the mapping matrix
235    if (use_bands_local) then
236      ! Cloud properties per band
237
238      allocate(mapping(this%nband, nwav))
239      allocate(weight(nwav))
240
241      ! Planck weight uses the wavenumbers of the cloud points
242      allocate(planck_weight(nwav))
243      planck_weight = calc_planck_function_wavenumber(wavenumber, &
244           &                                          temperature)
245
246      DO jband = 1,this%nband
247        weight = 0.0_jprb
248        DO jwav = 1,nwav
249          ! Work out wavenumber range for which this cloud wavenumber
250          ! will be applicable
251          if (wavenumber(jwav) >= this%wavenumber1_band(jband) &
252               & .AND. wavenumber(jwav) <= this%wavenumber2_band(jband)) then
253            if (jwav > 1) then
254              wavenum1 = max(this%wavenumber1_band(jband), &
255                   &  0.5_jprb*(wavenumber(jwav-1)+wavenumber(jwav)))
256            else
257              wavenum1 = this%wavenumber1_band(jband)
258            end if
259            if (jwav < nwav) then
260              wavenum2 = min(this%wavenumber2_band(jband), &
261                   &  0.5_jprb*(wavenumber(jwav)+wavenumber(jwav+1)))
262            else
263              wavenum2 = this%wavenumber2_band(jband)
264            end if
265            ! This cloud wavenumber is weighted by the wavenumber
266            ! range of its applicability multiplied by the Planck
267            ! function at an appropriate temperature
268            weight(jwav) = (wavenum2-wavenum1) * planck_weight(jwav)
269          end if
270        end do
271        if (sum(weight) <= 0.0_jprb) then
272          ! No cloud wavenumbers lie in the band; interpolate to
273          ! central wavenumber of band instead
274          if (wavenumber(1) >= this%wavenumber2_band(jband)) then
275            ! Band is entirely below first cloudy wavenumber
276            weight(1) = 1.0_jprb
277          else if (wavenumber(nwav) <= this%wavenumber1_band(jband)) then
278            ! Band is entirely above last cloudy wavenumber
279            weight(nwav) = 1.0_jprb
280          else
281            ! Find interpolating points
282            iwav = 2
283            DO while (wavenumber(iwav) < this%wavenumber2_band(jband))
284              iwav = iwav+1
285            end do
286            weight(iwav-1) = planck_weight(iwav-1) * (wavenumber(iwav) &
287                 &  - 0.5_jprb*(this%wavenumber2_band(jband)+this%wavenumber1_band(jband)))
288            weight(iwav) = planck_weight(iwav) * (-wavenumber(iwav-1) &
289                 &  + 0.5_jprb*(this%wavenumber2_band(jband)+this%wavenumber1_band(jband)))
290          end if
291        end if
292        mapping(jband,:) = weight / sum(weight)
293      end do
294
295      deallocate(weight)
296      deallocate(planck_weight)
297
298    else
299      ! Cloud properties per g-point
300
301      if (this%ng == 0) then
302        write(nulerr,'(a)') '*** Error: requested cloud/aerosol mapping per g-point but only available per band'
303        call radiation_abort('Radiation configuration error')
304      end if
305
306      allocate(mapping(this%ng, nwav))
307      allocate(weight(this%nwav))
308      allocate(planck_weight(this%nwav))
309
310      planck_weight = calc_planck_function_wavenumber(0.5_jprb &
311           &  * (this%wavenumber1 + this%wavenumber2), &
312           &  temperature)
313
314      mapping = 0.0_jprb
315      ! Loop over wavenumbers representing cloud
316      DO jwav = 1,nwav
317        ! Clear the weights. The weight says for one wavenumber in the
318        ! cloud file, what is its fractional contribution to each of
319        ! the spectral-definition intervals
320        weight = 0.0_jprb
321
322        ! Cloud properties are linearly interpolated between each of
323        ! the nwav cloud points; therefore, the influence of a
324        ! particular cloud point extends as a triangle between
325        ! wavenum0 and wavenum2, peaking at wavenum1
326        wavenum1 = wavenumber(jwav)
327        isd1 = this%find(wavenum1)
328        if (isd1 < 1) then
329          cycle
330        end if
331        if (jwav > 1) then
332          wavenum0 = wavenumber(jwav-1)
333
334          ! Map triangle under (wavenum0,0) to (wavenum1,1) to the
335          ! wavenumbers in this%gpoint_fraction
336          isd0 = this%find(wavenum0)
337          if (isd0 == isd1) then
338            ! Triangle completely within the range
339            ! this%wavenumber1(isd0)-this%wavenumber2(isd0)
340            weight(isd0) = 0.5_jprb*(wavenum1-wavenum0) &
341                 &       / (this%wavenumber2(isd0)-this%wavenumber1(isd0))
342          else
343            if (isd0 >= 1) then
344              ! Left part of triangle
345              weight(isd0) = 0.5_jprb * (this%wavenumber2(isd0)-wavenum0)**2 &
346                   &       / ((this%wavenumber2(isd0)-this%wavenumber1(isd0)) &
347                   &         *(wavenum1-wavenum0))
348            end if
349            ! Right part of triangle (trapezium)
350!            weight(isd1) = 0.5_jprb * (wavenum1-this%wavenumber1(isd1)) &
351!                 &       * (wavenum1 + this%wavenumber1(isd1) - 2.0_jprb*wavenum0) &
352!                 &       / (wavenum1-wavenum0)
353            weight(isd1) = 0.5_jprb * (1.0_jprb &
354                 &  + (this%wavenumber1(isd1)-wavenum1)/(wavenum1-wavenum0)) &
355                 &  * (wavenum1-this%wavenumber1(isd1)) &
356                 &  / (this%wavenumber2(isd1)-this%wavenumber1(isd1))
357            if (isd1-isd0 > 1) then
358              DO isd = isd0+1,isd1-1
359                ! Intermediate trapezia
360                weight(isd) = 0.5_jprb * (this%wavenumber1(isd)+this%wavenumber2(isd) &
361                     &                    - 2.0_jprb*wavenum0) &
362                     &      / (wavenum1-wavenum0)
363              end do
364            end if
365          end if
366
367        else
368          ! First cloud wavenumber: all wavenumbers in the spectral
369          ! definition below this will use the first one
370          if (isd1 >= 1) then
371            weight(1:isd1-1) = 1.0_jprb
372            weight(isd1) = (wavenum1-this%wavenumber1(isd1)) &
373                 &       / (this%wavenumber2(isd1)-this%wavenumber1(isd1))
374          end if
375        end if
376
377        if (jwav < nwav) then
378          wavenum2 = wavenumber(jwav+1)
379
380          ! Map triangle under (wavenum1,1) to (wavenum2,0) to the
381          ! wavenumbers in this%gpoint_fraction
382          isd2 = this%find(wavenum2)
383
384          if (isd1 == isd2) then
385            ! Triangle completely within the range
386            ! this%wavenumber1(isd1)-this%wavenumber2(isd1)
387            weight(isd1) = weight(isd1) + 0.5_jprb*(wavenum2-wavenum1) &
388                 &       / (this%wavenumber2(isd1)-this%wavenumber1(isd1))
389          else
390            if (isd2 >= 1 .AND. isd2 <= this%nwav) then
391              ! Right part of triangle
392              weight(isd2) = weight(isd2) + 0.5_jprb * (wavenum2-this%wavenumber1(isd2))**2 &
393                   &       / ((this%wavenumber2(isd2)-this%wavenumber1(isd2)) &
394                   &         *(wavenum2-wavenum1))
395            end if
396            ! Left part of triangle (trapezium)
397!            weight(isd1) = weight(isd1) + 0.5_jprb * (this%wavenumber2(isd1)-wavenum1) &
398!                 &       * (wavenum1 + this%wavenumber2(isd1) - 2.0_jprb*wavenum2) &
399!                 &       / (wavenum2-wavenum1)
400            weight(isd1) = weight(isd1) + 0.5_jprb * (1.0_jprb &
401                 &  + (wavenum2-this%wavenumber2(isd1)) / (wavenum2-wavenum1)) &
402                 &  * (this%wavenumber2(isd1)-wavenum1) &
403                 &  / (this%wavenumber2(isd1)-this%wavenumber1(isd1))
404            if (isd2-isd1 > 1) then
405              DO isd = isd1+1,isd2-1
406                ! Intermediate trapezia
407                weight(isd) = weight(isd) + 0.5_jprb * (2.0_jprb*wavenum2 &
408                     & - this%wavenumber1(isd) - this%wavenumber2(isd)) &
409                     &      / (wavenum2-wavenum1)
410              end do
411            end if
412          end if
413
414        else
415          ! Last cloud wavenumber: all wavenumbers in the spectral
416          ! definition above this will use the last one
417          if (isd1 <= this%nwav) then
418            weight(isd1+1:this%nwav) = 1.0_jprb
419            weight(isd1) = (this%wavenumber2(isd1)-wavenum1) &
420                 &       / (this%wavenumber2(isd1)-this%wavenumber1(isd1))
421          end if
422        end if
423
424        weight = weight * planck_weight
425
426        DO jg = 1,this%ng
427          mapping(jg, jwav) = sum(weight * this%gpoint_fraction(:,jg))
428        end do
429
430      end do
431
432      deallocate(weight)
433      deallocate(planck_weight)
434
435      ! Normalize mapping matrix
436      DO jg = 1,this%ng
437        mapping(jg,:) = mapping(jg,:) * (1.0_jprb/sum(mapping(jg,:)))
438      end do
439
440    end if
441
442    if (lhook) call dr_hook('radiation_spectral_definition:calc_mapping',1,hook_handle)
443
444  end subroutine calc_mapping
445
446
447  !---------------------------------------------------------------------
448  ! Under normal operation (if use_fluxes is .false. or not present),
449  ! compute a mapping matrix "mapping" that can be used in an
450  ! expression y=matmul(mapping^T,x) where x is a variable containing
451  ! optical properties in input bands (e.g. albedo in shortwave albedo
452  ! bands), and y is this variable mapped on to the spectral intervals
453  ! in the spectral definition "this". Temperature (K) is used to
454  ! generate a Planck function to weight each input band
455  ! appropriately. Note that "mapping" is here transposed from the
456  ! convention in the calc_mapping routine.  Under the alternative
457  ! operation (if use_fluxes is present and .true.), the mapping works
458  ! in the reverse sense: if y contains fluxes in each ecRad band or
459  ! g-point, then x=matmul(mapping,y) would return fluxes in x
460  ! averaged to user-supplied "input" bands. In this version, the
461  ! bands are described by their wavelength bounds (wavelength_bound,
462  ! which must be increasing and exclude the end points) and the index
463  ! of the mapping matrix that each band corresponds to (i_intervals,
464  ! which has one more element than wavelength_bound and can have
465  ! duplicated values if an albedo/emissivity value is to be
466  ! associated with more than one discontinuous ranges of the
467  ! spectrum).
468  subroutine calc_mapping_from_bands(this, temperature, &
469       &  wavelength_bound, i_intervals, mapping, use_bands, use_fluxes)
470
471    use yomhook,      only : lhook, dr_hook
472    use radiation_io, only : nulerr, radiation_abort
473
474    class(spectral_definition_type), intent(in)    :: this
475    real(jprb),                      intent(in)    :: temperature   ! K
476    ! Monotonically increasing wavelength bounds (m) between
477    ! intervals, not including the outer bounds (which are assumed to
478    ! be zero and infinity)
479    real(jprb), intent(in)    :: wavelength_bound(:)
480    ! The albedo band indices corresponding to each interval
481    integer,    intent(in)    :: i_intervals(:)
482    real(jprb), allocatable,         intent(inout) :: mapping(:,:)
483    logical,    optional,            intent(in)    :: use_bands
484    logical,    optional,            intent(in)    :: use_fluxes
485
486    ! Planck function and central wavenumber of each wavenumber
487    ! interval of the spectral definition
488    real(jprb) :: planck(this%nwav)         ! W m-2 (cm-1)-1
489    real(jprb) :: wavenumber_mid(this%nwav) ! cm-1
490
491    real(jprb), allocatable :: mapping_denom(:,:)
492
493    real(jprb) :: wavenumber1_bound, wavenumber2_bound
494
495    ! To work out weights we sample the Planck function at five points
496    ! in the interception between an input interval and a band, and
497    ! use the Trapezium Rule
498    integer, parameter :: nsample = 5
499    integer :: isamp
500    real(jprb), dimension(nsample) :: wavenumber_sample, planck_sample
501    real(jprb), parameter :: weight_sample(nsample) &
502         &        = [0.5_jprb, 1.0_jprb, 1.0_jprb, 1.0_jprb, 0.5_jprb]
503
504    ! Index of input value corresponding to each wavenumber interval
505    integer :: i_input(this%nwav)
506
507    ! Number of albedo/emissivity values that will be provided, some
508    ! of which may span discontinuous intervals in wavelength space
509    integer :: ninput
510
511    ! Number of albedo/emissivity intervals represented, where some
512    ! may be grouped to have the same value of albedo/emissivity (an
513    ! example is in the thermal infrared where classically the IFS has
514    ! ninput=2 and ninterval=3, since only two emissivities are
515    ! provided representing (1) the infrared window, and (2) the
516    ! intervals to each side of the infrared window.
517    integer :: ninterval
518
519    logical    :: use_bands_local, use_fluxes_local
520
521    ! Loop indices
522    integer    :: jg, jband, jin, jint
523
524    real(jprb) :: hook_handle
525
526    if (lhook) call dr_hook('radiation_spectral_definition:calc_mapping_from_bands',0,hook_handle)
527
528    if (present(use_bands)) then
529      use_bands_local = use_bands
530    else
531      use_bands_local = .false.
532    end if
533
534    if (present(use_fluxes)) then
535      use_fluxes_local = use_fluxes
536    else
537      use_fluxes_local = .false.
538    end if
539
540    ! Count the number of input intervals
541    ninterval = size(i_intervals)
542    ninput    = maxval(i_intervals)
543   
544    if (allocated(mapping)) then
545      deallocate(mapping)
546    end if
547   
548    ! Check wavelength is monotonically increasing
549    if (ninterval > 2) then
550      DO jint = 2,ninterval-1
551        if (wavelength_bound(jint) <= wavelength_bound(jint-1)) then
552          write(nulerr, '(a)') '*** Error: wavelength bounds must be monotonically increasing'
553          call radiation_abort()
554        end if
555      end do
556    end if
557
558    ! Define the mapping matrix
559    if (use_bands_local) then
560      ! Require properties per band
561
562      allocate(mapping(ninput, this%nband))
563      mapping = 0.0_jprb
564
565      if (use_fluxes_local) then
566        allocate(mapping_denom(ninput, this%nband))
567        mapping_denom = 0.0_jprb
568      end if
569
570      DO jband = 1,this%nband
571        DO jint = 1,ninterval
572          if (jint == 1) then
573            ! First input interval in wavelength space: lower
574            ! wavelength bound is 0 m, so infinity cm-1
575            wavenumber2_bound = this%wavenumber2_band(jband)
576          else
577            wavenumber2_bound = min(this%wavenumber2_band(jband), &
578                 &                  0.01_jprb/wavelength_bound(jint-1))
579          end if
580
581          if (jint == ninterval) then
582            ! Final input interval in wavelength space: upper
583            ! wavelength bound is infinity m, so 0 cm-1
584            wavenumber1_bound = this%wavenumber1_band(jband)
585          else
586            wavenumber1_bound = max(this%wavenumber1_band(jband), &
587                 &                  0.01_jprb/wavelength_bound(jint))
588
589          end if
590
591          if (wavenumber2_bound > wavenumber1_bound) then
592            ! Current input interval contributes to current band;
593            ! compute the weight of the contribution in proportion to
594            ! an approximate calculation of the integral of the Planck
595            ! function over the relevant part of the spectrum
596            wavenumber_sample = wavenumber1_bound + [(isamp,isamp=0,nsample-1)] &
597                 &  * (wavenumber2_bound-wavenumber1_bound) / real(nsample-1,jprb)
598            planck_sample = calc_planck_function_wavenumber(wavenumber_sample, temperature)
599            mapping(i_intervals(jint),jband) = mapping(i_intervals(jint),jband) &
600                 &  + sum(planck_sample*weight_sample) * (wavenumber2_bound-wavenumber1_bound)
601            if (use_fluxes_local) then
602              ! Compute an equivalent sample containing the entire ecRad band
603              wavenumber_sample = this%wavenumber1_band(jband) + [(isamp,isamp=0,nsample-1)] &
604                   &  * (this%wavenumber2_band(jband)-this%wavenumber1_band(jband)) &
605                   &  / real(nsample-1,jprb)
606              planck_sample = calc_planck_function_wavenumber(wavenumber_sample, temperature)
607              mapping_denom(i_intervals(jint),jband) = mapping_denom(i_intervals(jint),jband) &
608                 &  + sum(planck_sample*weight_sample) * (this%wavenumber2_band(jband)-this%wavenumber1_band(jband))
609            end if
610          end if
611
612        end do
613      end do
614
615      if (use_fluxes_local) then
616        mapping = mapping / max(1.0e-12_jprb, mapping_denom)
617        deallocate(mapping_denom)
618      end if
619
620    else
621      ! Require properties per g-point
622
623      if (this%ng == 0) then
624        write(nulerr,'(a)') '*** Error: requested surface mapping per g-point but only available per band'
625        call radiation_abort('Radiation configuration error')
626      end if
627
628      allocate(mapping(ninput,this%ng))
629      mapping = 0.0_jprb
630
631      wavenumber_mid = 0.5_jprb * (this%wavenumber1 + this%wavenumber2)
632      planck = calc_planck_function_wavenumber(wavenumber_mid, temperature)
633
634      ! In the processing that follows, we assume that the wavenumber
635      ! grid on which the g-points are defined in the spectral
636      ! definition is much finer than the albedo/emissivity intervals
637      ! that the user will provide.  This means that each wavenumber
638      ! is assigned to only one of the albedo/emissivity intervals.
639
640      ! By default set all wavenumbers to use first input
641      ! albedo/emissivity
642      i_input = 1
643     
644      ! All bounded intervals
645      DO jint = 2,ninterval-1
646        wavenumber1_bound = 0.01_jprb / wavelength_bound(jint)
647        wavenumber2_bound = 0.01_jprb / wavelength_bound(jint-1)
648        where (wavenumber_mid > wavenumber1_bound &
649             & .AND. wavenumber_mid <= wavenumber2_bound)
650          i_input = i_intervals(jint)
651        end where
652      end do
653
654      ! Final interval in wavelength space goes up to wavenumber of
655      ! infinity
656      if (ninterval > 1) then
657        wavenumber2_bound = 0.01_jprb / wavelength_bound(ninterval-1)
658        where (wavenumber_mid <= wavenumber2_bound)
659          i_input = i_intervals(ninterval)
660        end where
661      end if
662
663      DO jg = 1,this%ng
664        DO jin = 1,ninput
665          mapping(jin,jg) = sum(this%gpoint_fraction(:,jg) * planck, &
666               &                 mask=(i_input==jin))
667          if (use_fluxes_local) then
668            mapping(jin,jg) = mapping(jin,jg) / sum(this%gpoint_fraction(:,jg) * planck)
669          end if
670        end do
671      end do
672
673    end if
674
675    if (.not. use_fluxes_local) then
676      ! Normalize mapping matrix
677      DO jg = 1,size(mapping,dim=2)
678        mapping(:,jg) = mapping(:,jg) * (1.0_jprb/sum(mapping(:,jg)))
679      end do
680    end if
681
682    if (lhook) call dr_hook('radiation_spectral_definition:calc_mapping_from_bands',1,hook_handle)
683
684  end subroutine calc_mapping_from_bands
685
686
687  !---------------------------------------------------------------------
688  ! As calc_mapping_from_bands but in terms of wavenumber bounds from
689  ! wavenumber1 to wavenumber2
690  subroutine calc_mapping_from_wavenumber_bands(this, temperature, &
691       &  wavenumber1, wavenumber2, mapping, use_bands, use_fluxes)
692
693    use yomhook,      only : lhook, dr_hook
694    use radiation_io, only : nulerr, radiation_abort
695
696    class(spectral_definition_type), intent(in)    :: this
697    real(jprb),                      intent(in)    :: temperature   ! K
698    real(jprb), intent(in)    :: wavenumber1(:), wavenumber2(:)
699    real(jprb), allocatable,         intent(inout) :: mapping(:,:)
700    logical,    optional,            intent(in)    :: use_bands
701    logical,    optional,            intent(in)    :: use_fluxes
702
703    ! Monotonically increasing wavelength bounds (m) between
704    ! intervals, not including the outer bounds (which are assumed to
705    ! be zero and infinity)
706    real(jprb) :: wavelength_bound(size(wavenumber1)-1)
707    ! The albedo band indices corresponding to each interval
708    integer    :: i_intervals(size(wavenumber1))
709
710    ! Lower wavelength bound (m) of each band
711    real(jprb) :: wavelength1(size(wavenumber1))
712
713    logical    :: is_band_unassigned(size(wavenumber1))
714
715    ! Number of albedo/emissivity intervals represented, where some
716    ! may be grouped to have the same value of albedo/emissivity (an
717    ! example is in the thermal infrared where classically the IFS has
718    ! ninput=2 and ninterval=3, since only two emissivities are
719    ! provided representing (1) the infrared window, and (2) the
720    ! intervals to each side of the infrared window.
721    integer :: ninterval
722
723    ! Index to next band in order of increasing wavelength
724    integer :: inext
725
726    ! Loop indices
727    integer :: jint
728
729    real(jprb) :: hook_handle
730
731    if (lhook) call dr_hook('radiation_spectral_definition:calc_mapping_from_wavenumber_bands',0,hook_handle)
732
733    wavelength1 = 0.01_jprb / wavenumber2
734    ninterval = size(wavelength1)
735   
736    is_band_unassigned = .true.
737
738    DO jint = 1,ninterval
739      inext = minloc(wavelength1, dim=1, mask=is_band_unassigned)
740      if (jint > 1) then
741        wavelength_bound(jint-1) = wavelength1(inext)
742      end if
743      is_band_unassigned(inext) = .false.
744      i_intervals(jint) = inext
745    end do
746
747    call calc_mapping_from_bands(this, temperature, &
748         &  wavelength_bound, i_intervals, mapping, use_bands, use_fluxes)
749
750    if (lhook) call dr_hook('radiation_spectral_definition:calc_mapping_from_wavenumber_bands',1,hook_handle)
751
752  end subroutine calc_mapping_from_wavenumber_bands
753
754
755  !---------------------------------------------------------------------
756  ! Print out the mapping computed by calc_mapping_from_bands
757  subroutine print_mapping_from_bands(this, mapping, use_bands)
758
759    use radiation_io, only : nulout
760
761    class(spectral_definition_type), intent(in) :: this
762    real(jprb), allocatable,         intent(in) :: mapping(:,:) ! (ninput,nband/ng)
763    logical,    optional,            intent(in) :: use_bands
764
765    logical :: use_bands_local
766
767    integer :: nin, nout
768    integer :: jin, jout
769
770    if (present(use_bands)) then
771      use_bands_local = use_bands
772    else
773      use_bands_local = .false.
774    end if
775
776    nin = size(mapping,1)
777    nout = size(mapping,2)
778
779    if (nin <= 1) then
780      write(nulout, '(a)') '  All spectral intervals will use the same albedo/emissivity'
781    else if (use_bands_local) then
782      write(nulout, '(a,i0,a,i0,a)') '  Mapping from ', nin, ' values to ', nout, ' bands (wavenumber ranges in cm-1)'
783      if (nout <= 40) then
784        DO jout = 1,nout
785          write(nulout,'(i6,a,i6,a)',advance='no') nint(this%wavenumber1_band(jout)), ' to', &
786               &                        nint(this%wavenumber2_band(jout)), ':'
787          DO jin = 1,nin
788            write(nulout,'(f5.2)',advance='no') mapping(jin,jout)
789          end do
790          write(nulout, '()')
791        end do
792      else
793        DO jout = 1,30
794          write(nulout,'(i6,a,i6,a)',advance='no') nint(this%wavenumber1_band(jout)), ' to', &
795               &                        nint(this%wavenumber2_band(jout)), ':'
796          DO jin = 1,nin
797            write(nulout,'(f5.2)',advance='no') mapping(jin,jout)
798          end do
799          write(nulout, '()')
800        end do
801        write(nulout,'(a)') '  ...'
802        write(nulout,'(i6,a,i6,a)',advance='no') nint(this%wavenumber1_band(nout)), ' to', &
803             &                        nint(this%wavenumber2_band(nout)), ':'
804        DO jin = 1,nin
805          write(nulout,'(f5.2)',advance='no') mapping(jin,nout)
806        end do
807        write(nulout, '()')
808      end if
809    else
810      write(nulout, '(a,i0,a,i0,a)') '  Mapping from ', nin, ' values to ', nout, ' g-points'
811      if (nout <= 40) then
812        DO jout = 1,nout
813          write(nulout,'(i3,a)',advance='no') jout, ':'
814          DO jin = 1,nin
815            write(nulout,'(f5.2)',advance='no') mapping(jin,jout)
816          end do
817          write(nulout, '()')
818        end do
819      else
820        DO jout = 1,30
821          write(nulout,'(i3,a)',advance='no') jout, ':'
822          DO jin = 1,nin
823            write(nulout,'(f5.2)',advance='no') mapping(jin,jout)
824          end do
825          write(nulout, '()')
826        end do
827        write(nulout,'(a)') '  ...'
828        write(nulout,'(i3,a)',advance='no') nout, ':'
829        DO jin = 1,nin
830          write(nulout,'(f5.2)',advance='no') mapping(jin,nout)
831        end do
832        write(nulout, '()')
833      end if
834    end if
835
836  end subroutine print_mapping_from_bands
837
838
839  !---------------------------------------------------------------------
840  ! Return the minimum wavenumber of this object in cm-1
841  pure function min_wavenumber(this)
842
843    class(spectral_definition_type), intent(in)    :: this
844    real(jprb) :: min_wavenumber
845
846    if (this%nwav > 0) then
847      min_wavenumber = this%wavenumber1(1)
848    else
849      min_wavenumber = minval(this%wavenumber1_band)
850    end if
851
852  end function min_wavenumber
853
854
855  !---------------------------------------------------------------------
856  ! Return the maximum wavenumber of this object in cm-1
857  pure function max_wavenumber(this)
858
859    class(spectral_definition_type), intent(in)    :: this
860    real(jprb) :: max_wavenumber
861
862    if (this%nwav > 0) then
863      max_wavenumber = this%wavenumber1(this%nwav)
864    else
865      max_wavenumber = maxval(this%wavenumber2_band)
866    end if
867
868  end function max_wavenumber
869
870
871  !---------------------------------------------------------------------
872  ! Return the Planck function (in W m-2 (cm-1)-1) for a given
873  ! wavenumber (cm-1) and temperature (K), ensuring double precision
874  ! for internal calculation.  If temperature is 0 or less then unity
875  ! is returned; since this function is primarily used to weight an
876  ! integral by the Planck function, a temperature of 0 or less means
877  ! no weighting is to be applied.
878  elemental function calc_planck_function_wavenumber(wavenumber, temperature)
879
880    use parkind1,            only : jprb, jprd
881    use radiation_constants, only : SpeedOfLight, BoltzmannConstant, PlanckConstant
882
883    real(jprb), intent(in) :: wavenumber  ! cm-1
884    real(jprb), intent(in) :: temperature ! K
885    real(jprb) :: calc_planck_function_wavenumber
886
887    real(jprd) :: freq ! Hz
888    real(jprd) :: planck_fn_freq ! W m-2 Hz-1
889
890    if (temperature > 0.0_jprd) then
891      freq = 100.0_jprd * real(SpeedOfLight,jprd) * real(wavenumber,jprd)
892      planck_fn_freq = 2.0_jprd * real(PlanckConstant,jprd) * freq**3 &
893           &  / (real(SpeedOfLight,jprd)**2 * (exp(real(PlanckConstant,jprd)*freq &
894           &     /(real(BoltzmannConstant,jprd)*real(temperature,jprd))) - 1.0_jprd))
895      calc_planck_function_wavenumber = real(planck_fn_freq * 100.0_jprd * real(SpeedOfLight,jprd), jprb)
896    else
897      calc_planck_function_wavenumber = 1.0_jprb
898    end if
899
900  end function calc_planck_function_wavenumber
901
902end module radiation_spectral_definition
Note: See TracBrowser for help on using the repository browser.