source: LMDZ6/trunk/libf/phylmd/ecrad/radiation_ecckd.F90 @ 4582

Last change on this file since 4582 was 4489, checked in by lguez, 21 months ago

Merge LMDZ_ECRad branch back into trunk!

File size: 20.3 KB
Line 
1! radiation_ecckd.F90 - ecCKD generalized gas optics model
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_ecckd
18
19  use parkind1, only : jprb
20  use radiation_gas_constants
21  use radiation_ecckd_gas
22  use radiation_spectral_definition, only : spectral_definition_type
23
24  implicit none
25
26  public
27
28  !---------------------------------------------------------------------
29  ! This derived type contains all the data needed to describe a
30  ! correlated k-distribution gas optics model created using the ecCKD
31  ! tool
32  type ckd_model_type
33
34    ! Gas information
35
36    ! Number of gases
37    integer :: ngas = 0
38    ! Array of individual-gas data objects
39    type(ckd_gas_type), allocatable :: single_gas(:)
40    ! Mapping from the "gas codes" in the radiation_gas_constants
41    ! module to an index to the single_gas array, where zero means gas
42    ! not present (or part of a composite gas)
43    integer :: i_gas_mapping(0:NMaxGases)
44
45    ! Coordinates of main look-up table for absorption coeffts
46
47    ! Number of pressure and temperature points
48    integer :: npress = 0
49    integer :: ntemp  = 0
50    ! Natural logarithm of first (lowest) pressure (Pa) and increment
51    real(jprb) :: log_pressure1, d_log_pressure
52    ! First temperature profile (K), dimensioned (npress)
53    real(jprb), allocatable :: temperature1(:)
54    ! Temperature increment (K)
55    real(jprb) :: d_temperature
56
57    ! Look-up table for Planck function
58
59    ! Number of entries
60    integer :: nplanck = 0
61    ! Temperature of first element of look-up table and increment (K)
62    real(jprb), allocatable :: temperature1_planck
63    real(jprb), allocatable :: d_temperature_planck
64    ! Planck function (black body flux into a horizontal plane) in W
65    ! m-2, dimensioned (ng,nplanck)
66    real(jprb), allocatable :: planck_function(:,:)
67
68    ! Normalized solar irradiance in each g point dimensioned (ng)
69    real(jprb), allocatable :: norm_solar_irradiance(:)
70
71    ! Rayleigh molar scattering coefficient in m2 mol-1 in each g
72    ! point
73    real(jprb), allocatable :: rayleigh_molar_scat(:)
74
75    ! ! Spectral mapping of g points
76
77    ! ! Number of wavenumber intervals
78    ! integer :: nwav = 0
79
80    ! Number of k terms / g points
81    integer :: ng   = 0
82
83    ! Spectral definition describing bands and g points
84    type(spectral_definition_type) :: spectral_def
85
86    ! Shortwave: true, longwave: false
87    logical :: is_sw
88
89  contains
90
91    procedure :: read => read_ckd_model
92    procedure :: calc_optical_depth => calc_optical_depth_ckd_model
93    procedure :: print => print_ckd_model
94    procedure :: calc_planck_function
95    procedure :: calc_incoming_sw
96!    procedure :: deallocate => deallocate_ckd_model
97
98  end type ckd_model_type
99
100
101contains
102
103  !---------------------------------------------------------------------
104  ! Read a complete ecCKD gas optics model from a NetCDF file
105  ! "filename"
106  subroutine read_ckd_model(this, filename, iverbose)
107
108    use easy_netcdf,  only : netcdf_file
109    !use radiation_io, only : nulerr, radiation_abort
110    use yomhook,      only : lhook, dr_hook
111
112    class(ckd_model_type), intent(inout) :: this
113    character(len=*),      intent(in)    :: filename
114    integer, optional,     intent(in)    :: iverbose
115
116    type(netcdf_file) :: file
117
118    real(jprb), allocatable :: pressure_lut(:)
119    real(jprb), allocatable :: temperature_full(:,:)
120    real(jprb), allocatable :: temperature_planck(:)
121
122    character(len=512) :: constituent_id
123
124    integer :: iverbose_local
125
126    ! Loop counters
127    integer :: jgas, jjgas
128
129    integer :: istart, inext, nchar, i_gas_code
130
131    real(jprb)         :: hook_handle
132
133    if (lhook) call dr_hook('radiation_ecckd:read_ckd_model',0,hook_handle)
134
135    if (present(iverbose)) then
136      iverbose_local = iverbose
137    else
138      iverbose_local = 3
139    end if
140
141    call file%open(trim(filename), iverbose=iverbose_local)
142
143    ! Read temperature and pressure coordinate variables
144    call file%get('pressure', pressure_lut)
145    this%log_pressure1 = log(pressure_lut(1))
146    this%npress = size(pressure_lut)
147    this%d_log_pressure = log(pressure_lut(2)) - this%log_pressure1
148    call file%get('temperature', temperature_full)
149    allocate(this%temperature1(this%npress));
150    this%temperature1 = temperature_full(:,1)
151    this%d_temperature = temperature_full(1,2)-temperature_full(1,1)
152    this%ntemp = size(temperature_full,2)
153    deallocate(temperature_full)
154   
155    ! Read Planck function, or solar irradiance and Rayleigh
156    ! scattering coefficient
157    if (file%exists('solar_irradiance')) then
158      this%is_sw = .true.
159      call file%get('solar_irradiance', this%norm_solar_irradiance)
160      this%norm_solar_irradiance = this%norm_solar_irradiance &
161           &  / sum(this%norm_solar_irradiance)
162      call file%get('rayleigh_molar_scattering_coeff', &
163           &  this%rayleigh_molar_scat)
164    else
165      this%is_sw = .false.
166      call file%get('temperature_planck', temperature_planck)
167      this%nplanck = size(temperature_planck)
168      this%temperature1_planck = temperature_planck(1)
169      this%d_temperature_planck = temperature_planck(2) - temperature_planck(1)
170      deallocate(temperature_planck)
171      call file%get('planck_function', this%planck_function)
172    end if
173
174    ! Read the spectral definition information into a separate
175    ! derived type
176    call this%spectral_def%read(file);
177    this%ng = this%spectral_def%ng
178
179    ! Read gases
180    call file%get('n_gases', this%ngas)
181    allocate(this%single_gas(this%ngas))
182    call file%get_global_attribute('constituent_id', constituent_id)
183    nchar = len(trim(constituent_id))
184    istart = 1
185    this%i_gas_mapping = 0
186    do jgas = 1, this%ngas
187      if (jgas < this%ngas) then
188        inext = istart + scan(constituent_id(istart:nchar), ' ')
189      else
190        inext = nchar+2
191      end if
192      ! Find gas code
193      i_gas_code = 0
194      do jjgas = 1, NMaxGases
195        if (constituent_id(istart:inext-2) == trim(GasLowerCaseName(jjgas))) then
196          i_gas_code = jjgas
197          exit
198        end if
199      end do
200      ! if (i_gas_code == 0) then
201      !   write(nulerr,'(a,a,a)') '*** Error: Gas "', constituent_id(istart:inext-2), &
202      !        & '" not understood'
203      !   call radiation_abort('Radiation configuration error')
204      ! end if
205      this%i_gas_mapping(i_gas_code) = jgas;
206      call this%single_gas(jgas)%read(file, constituent_id(istart:inext-2), i_gas_code)
207      istart = inext
208    end do
209   
210    if (lhook) call dr_hook('radiation_ecckd:read_ckd_model',1,hook_handle)
211
212  end subroutine read_ckd_model
213
214  !---------------------------------------------------------------------
215  ! Print a description of the correlated k-distribution model to the
216  ! "nulout" unit
217  subroutine print_ckd_model(this)
218
219    use radiation_io, only : nulout
220    use radiation_gas_constants
221
222    class(ckd_model_type), intent(in)  :: this
223
224    integer :: jgas
225   
226    if (this%is_sw) then
227      write(nulout,'(a)',advance='no') 'ecCKD shortwave gas optics model: '
228    else
229      write(nulout,'(a)',advance='no') 'ecCKD longwave gas optics model: '
230    end if
231
232    write(nulout,'(i0,a,i0,a,i0,a,i0,a)') &
233         &  nint(this%spectral_def%wavenumber1(1)), '-', &
234         &  nint(this%spectral_def%wavenumber2(size(this%spectral_def%wavenumber2))), &
235         &  ' cm-1, ', this%ng, ' g-points in ', this%spectral_def%nband, ' bands'
236    write(nulout,'(a,i0,a,i0,a,i0,a)') '  Look-up table sizes: ', this%npress, ' pressures, ', &
237         &  this%ntemp, ' temperatures, ', this%nplanck, ' planck-function entries'
238    write(nulout, '(a)') '  Gases:'
239    do jgas = 1,this%ngas
240      if (this%single_gas(jgas)%i_gas_code > 0) then
241        write(nulout, '(a,a,a)', advance='no') '    ', &
242             &  trim(GasName(this%single_gas(jgas)%i_gas_code)), ': '
243      else
244        write(nulout, '(a)', advance='no') '    Composite of well-mixed background gases: '
245      end if
246      select case (this%single_gas(jgas)%i_conc_dependence)
247        case (IConcDependenceNone)
248          write(nulout, '(a)') 'no concentration dependence'
249        case (IConcDependenceLinear)
250          write(nulout, '(a)') 'linear concentration dependence'
251        case (IConcDependenceRelativeLinear)
252          write(nulout, '(a,e14.6)') 'linear concentration dependence relative to a mole fraction of ', &
253               &  this%single_gas(jgas)%reference_mole_frac
254        case (IConcDependenceLUT)
255          write(nulout, '(a,i0,a,e14.6,a,e13.6)') 'look-up table with ', this%single_gas(jgas)%n_mole_frac, &
256               &  ' log-spaced mole fractions in range ', exp(this%single_gas(jgas)%log_mole_frac1), &
257               &  ' to ', exp(this%single_gas(jgas)%log_mole_frac1 &
258               &           + this%single_gas(jgas)%n_mole_frac*this%single_gas(jgas)%d_log_mole_frac)
259      end select
260    end do
261
262  end subroutine print_ckd_model
263
264
265  !---------------------------------------------------------------------
266  ! Compute layerwise optical depth for each g point for ncol columns
267  ! at nlev layers
268  subroutine calc_optical_depth_ckd_model(this, ncol, nlev, istartcol, iendcol, nmaxgas, &
269       &  pressure_hl, temperature_fl, mole_fraction_fl, &
270       &  optical_depth_fl, rayleigh_od_fl)
271
272    use yomhook,             only : lhook, dr_hook
273    use radiation_constants, only : AccelDueToGravity
274
275    ! Input variables
276
277    class(ckd_model_type), intent(in), target  :: this
278    ! Number of columns, levels and input gases
279    integer,               intent(in)  :: ncol, nlev, nmaxgas, istartcol, iendcol
280    ! Pressure at half levels (Pa), dimensioned (ncol,nlev+1)
281    real(jprb),            intent(in)  :: pressure_hl(ncol,nlev+1)
282    ! Temperature at full levels (K), dimensioned (ncol,nlev)
283    real(jprb),            intent(in)  :: temperature_fl(istartcol:iendcol,nlev)
284    ! Gas mole fractions at full levels (mol mol-1), dimensioned (ncol,nlev,nmaxgas)
285    real(jprb),            intent(in)  :: mole_fraction_fl(ncol,nlev,nmaxgas)
286   
287    ! Output variables
288
289    ! Layer absorption optical depth for each g point
290    real(jprb),            intent(out) :: optical_depth_fl(this%ng,nlev,istartcol:iendcol)
291    ! In the shortwave only, the Rayleigh scattering optical depth
292    real(jprb),  optional, intent(out) :: rayleigh_od_fl(this%ng,nlev,istartcol:iendcol)
293
294    ! Local variables
295
296    real(jprb), pointer :: molar_abs(:,:,:), molar_abs_conc(:,:,:,:)
297
298    ! Natural logarithm of pressure at full levels
299    real(jprb) :: log_pressure_fl(nlev)
300
301    ! Optical depth of single gas at one point in space versus
302    ! spectral interval
303    !real(jprb) :: od_single_gas(this%ng)
304
305    real(jprb) :: multiplier(nlev), simple_multiplier(nlev), global_multiplier, temperature1
306
307    ! Indices and weights in temperature, pressure and concentration interpolation
308    real(jprb) :: pindex1, tindex1, cindex1
309    real(jprb) :: pw1(nlev), pw2(nlev), tw1(nlev), tw2(nlev), cw1(nlev), cw2(nlev)
310    integer    :: ip1(nlev), it1(nlev), ic1(nlev)
311
312    ! Natural logarithm of mole fraction at one point
313    real(jprb) :: log_conc
314
315    ! Minimum mole fraction in look-up-table
316    real(jprb) :: mole_frac1
317
318    integer :: jcol, jlev, jgas, igascode
319
320    real(jprb) :: hook_handle
321
322    if (lhook) call dr_hook('radiation_ecckd:calc_optical_depth',0,hook_handle)
323
324    global_multiplier = 1.0_jprb / (AccelDueToGravity * 0.001_jprb * AirMolarMass)
325
326    do jcol = istartcol,iendcol
327
328      log_pressure_fl = log(0.5_jprb * (pressure_hl(jcol,1:nlev)+pressure_hl(jcol,2:nlev+1)))
329
330      do jlev = 1,nlev
331        ! Find interpolation points in pressure
332        pindex1 = (log_pressure_fl(jlev)-this%log_pressure1) &
333             &    / this%d_log_pressure
334        pindex1 = 1.0_jprb + max(0.0_jprb, min(pindex1, this%npress-1.0001_jprb))
335        ip1(jlev) = int(pindex1)
336        pw2(jlev) = pindex1 - ip1(jlev)
337        pw1(jlev) = 1.0_jprb - pw2(jlev)
338
339        ! Find interpolation points in temperature
340        temperature1 = pw1(jlev)*this%temperature1(ip1(jlev)) &
341             &       + pw2(jlev)*this%temperature1(ip1(jlev)+1)
342        tindex1 = (temperature_fl(jcol,jlev) - temperature1) &
343             &    / this%d_temperature
344        tindex1 = 1.0_jprb + max(0.0_jprb, min(tindex1, this%ntemp-1.0001_jprb))
345        it1(jlev) = int(tindex1)
346        tw2(jlev) = tindex1 - it1(jlev)
347        tw1(jlev) = 1.0_jprb - tw2(jlev)
348
349        ! Concentration multiplier
350        simple_multiplier(jlev) = global_multiplier &
351             &  * (pressure_hl(jcol,jlev+1) - pressure_hl(jcol,jlev))
352      end do
353
354      optical_depth_fl(:,:,jcol) = 0.0_jprb
355     
356      do jgas = 1,this%ngas
357
358        associate (single_gas => this%single_gas(jgas))
359          igascode = this%single_gas(jgas)%i_gas_code
360         
361          select case (single_gas%i_conc_dependence)
362           
363          case (IConcDependenceLinear)
364            molar_abs => this%single_gas(jgas)%molar_abs
365            multiplier = simple_multiplier * mole_fraction_fl(jcol,:,igascode)
366
367            do jlev = 1,nlev
368              optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) &
369                   &        + (multiplier(jlev)*tw1(jlev)) * (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)) &
370                   &                +pw2(jlev) * molar_abs(:,ip1(jlev)+1,it1(jlev))) &
371                   &        + (multiplier(jlev)*tw2(jlev)) * (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)+1) &
372                   &                +pw2(jlev) * molar_abs(:,ip1(jlev)+1,it1(jlev)+1))
373            end do
374
375          case (IConcDependenceRelativeLinear)
376            molar_abs => this%single_gas(jgas)%molar_abs
377            multiplier = simple_multiplier  * (mole_fraction_fl(jcol,:,igascode) &
378                 &                            - single_gas%reference_mole_frac)
379            do jlev = 1,nlev
380              optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) &
381                   &        + (multiplier(jlev)*tw1(jlev)) * (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)) &
382                   &                +pw2(jlev) * molar_abs(:,ip1(jlev)+1,it1(jlev))) &
383                   &        + (multiplier(jlev)*tw2(jlev)) * (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)+1) &
384                   &                +pw2(jlev) * molar_abs(:,ip1(jlev)+1,it1(jlev)+1))
385            end do
386
387          case (IConcDependenceNone)
388            ! Composite gases
389            molar_abs => this%single_gas(jgas)%molar_abs
390            do jlev = 1,nlev
391              optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) &
392                   &  + (simple_multiplier(jlev)*tw1(jlev)) * (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)) &
393                   &                              +pw2(jlev) * molar_abs(:,ip1(jlev)+1,it1(jlev))) &
394                   &  + (simple_multiplier(jlev)*tw2(jlev)) * (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)+1) &
395                   &                              +pw2(jlev) * molar_abs(:,ip1(jlev)+1,it1(jlev)+1))
396            end do
397
398          case (IConcDependenceLUT)
399            ! Logarithmic interpolation in concentration space
400            molar_abs_conc => this%single_gas(jgas)%molar_abs_conc
401            mole_frac1 = exp(single_gas%log_mole_frac1)
402            do jlev = 1,nlev
403              ! Take care of mole_fraction == 0
404              log_conc = log(max(mole_fraction_fl(jcol,jlev,igascode), mole_frac1))
405              cindex1  = (log_conc - single_gas%log_mole_frac1) / single_gas%d_log_mole_frac
406              cindex1  = 1.0_jprb + max(0.0_jprb, min(cindex1, single_gas%n_mole_frac-1.0001_jprb))
407              ic1(jlev) = int(cindex1)
408              cw2(jlev) = cindex1 - ic1(jlev)
409              cw1(jlev) = 1.0_jprb - cw2(jlev)
410            end do
411              ! od_single_gas = cw1 * (tw1 * (pw1 * molar_abs_conc(:,ip1,it1,ic1) &
412              !      &                       +pw2 * molar_abs_conc(:,ip1+1,it1,ic1)) &
413              !      &                +tw2 * (pw1 * molar_abs_conc(:,ip1,it1+1,ic1) &
414              !      &                       +pw2 * molar_abs_conc(:,ip1+1,it1+1,ic1))) &
415              !      &         +cw2 * (tw1 * (pw1 * molar_abs_conc(:,ip1,it1,ic1+1) &
416              !      &                       +pw2 * molar_abs_conc(:,ip1+1,it1,ic1+1)) &
417              !      &                +tw2 * (pw1 * molar_abs_conc(:,ip1,it1+1,ic1+1) &
418              !      &                       +pw2 * molar_abs_conc(:,ip1+1,it1+1,ic1+1)))
419            do jlev = 1,nlev
420              optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) &
421                   &  + (simple_multiplier(jlev) * mole_fraction_fl(jcol,jlev,igascode)) * ( &
422                   &      (cw1(jlev) * tw1(jlev) * pw1(jlev)) * molar_abs_conc(:,ip1(jlev),it1(jlev),ic1(jlev)) &
423                   &     +(cw1(jlev) * tw1(jlev) * pw2(jlev)) * molar_abs_conc(:,ip1(jlev)+1,it1(jlev),ic1(jlev)) &
424                   &     +(cw1(jlev) * tw2(jlev) * pw1(jlev)) * molar_abs_conc(:,ip1(jlev),it1(jlev)+1,ic1(jlev)) &
425                   &     +(cw1(jlev) * tw2(jlev) * pw2(jlev)) * molar_abs_conc(:,ip1(jlev)+1,it1(jlev)+1,ic1(jlev)) &
426                   &     +(cw2(jlev) * tw1(jlev) * pw1(jlev)) * molar_abs_conc(:,ip1(jlev),it1(jlev),ic1(jlev)+1) &
427                   &     +(cw2(jlev) * tw1(jlev) * pw2(jlev)) * molar_abs_conc(:,ip1(jlev)+1,it1(jlev),ic1(jlev)+1) &
428                   &     +(cw2(jlev) * tw2(jlev) * pw1(jlev)) * molar_abs_conc(:,ip1(jlev),it1(jlev)+1,ic1(jlev)+1) &
429                   &     +(cw2(jlev) * tw2(jlev) * pw2(jlev)) * molar_abs_conc(:,ip1(jlev)+1,it1(jlev)+1,ic1(jlev)+1))
430            end do
431          end select
432           
433        end associate
434
435      end do
436
437      ! Ensure the optical depth is not negative
438      optical_depth_fl(:,:,jcol) = max(0.0_jprb, optical_depth_fl(:,:,jcol))
439
440      ! Rayleigh scattering
441      if (this%is_sw .and. present(rayleigh_od_fl)) then
442        do jlev = 1,nlev
443          rayleigh_od_fl(:,jlev,jcol) = global_multiplier &
444               &  * (pressure_hl(jcol,jlev+1) - pressure_hl(jcol,jlev)) * this%rayleigh_molar_scat
445        end do
446      end if
447
448    end do
449
450    if (lhook) call dr_hook('radiation_ecckd:calc_optical_depth',1,hook_handle)
451
452  end subroutine calc_optical_depth_ckd_model
453
454  !---------------------------------------------------------------------
455  ! Calculate the Planck function integrated across each of the g
456  ! points of this correlated k-distribution model, for a given
457  ! temperature, where Planck function is defined as the flux emitted
458  ! by a black body (rather than radiance)
459  subroutine calc_planck_function(this, nt, temperature, planck)
460
461    class(ckd_model_type), intent(in)  :: this
462    integer,    intent(in)  :: nt
463    real(jprb), intent(in)  :: temperature(:) ! K
464    real(jprb), intent(out) :: planck(this%ng,nt) ! W m-2
465
466    real(jprb) :: tindex1, tw1, tw2
467    integer    :: it1, jt
468
469    do jt = 1,nt
470      tindex1 = (temperature(jt) - this%temperature1_planck) &
471           &   * (1.0_jprb / this%d_temperature_planck)
472      if (tindex1 >= 0) then
473        ! Normal interpolation, and extrapolation for high temperatures
474        tindex1 = 1.0_jprb + tindex1
475        it1 = min(int(tindex1), this%nplanck-1)
476        tw2 = tindex1 - it1
477        tw1 = 1.0_jprb - tw2
478        planck(:,jt) = tw1 * this%planck_function(:,it1) &
479             &       + tw2 * this%planck_function(:,it1+1)
480      else
481        ! Interpolate linearly to zero
482        planck(:,jt) = this%planck_function(:,1) &
483             &       * (temperature(jt)/this%temperature1_planck)
484      end if
485    end do
486
487  end subroutine calc_planck_function
488 
489
490  !---------------------------------------------------------------------
491  ! Return the spectral solar irradiance integrated over each g point
492  ! of a solar correlated k-distribution model, given the
493  ! total_solar_irradiance
494  subroutine calc_incoming_sw(this, total_solar_irradiance, &
495       &                      spectral_solar_irradiance)
496
497    class(ckd_model_type), intent(in)    :: this
498    real(jprb),            intent(in)    :: total_solar_irradiance ! W m-2
499    real(jprb),            intent(inout) :: spectral_solar_irradiance(:,:) ! W m-2
500 
501    spectral_solar_irradiance &
502         &  = spread(total_solar_irradiance * this%norm_solar_irradiance, &
503         &           2, size(spectral_solar_irradiance,2))
504
505  end subroutine calc_incoming_sw
506
507end module radiation_ecckd
Note: See TracBrowser for help on using the repository browser.