source: LMDZ6/branches/LMDZ_ECRad/libf/phylmd/ecrad/radiation/radiation_ecckd.F90 @ 5450

Last change on this file since 5450 was 4728, checked in by idelkadi, 15 months ago

Update of ecrad in the LMDZ_ECRad branch of LMDZ:

  • version 1.6.1 of ecrad
  • files are no longer grouped in the same ecrad directory.
  • the structure of ecrad offline is preserved to facilitate updating in LMDZ
  • cfg.bld modified to take into account the new added subdirectories.
  • the interface routines and those added in ecrad are moved to the phylmd directory
File size: 39.2 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    ! Normalized amplitude of variations in the solar irradiance
72    ! through the solar cycle in each g point, dimensioned (ng).
73    ! Since the user always provides the solar irradiance SI
74    ! integrated across the spectrum, this variable must sum to zero:
75    ! this ensures that the solar irradiance in each g-point is
76    ! SSI=SI*(norm_solar_irradiance +
77    ! A*norm_amplitude_solar_irradiance) for any A, where A denotes
78    ! the amplitude of deviations from the mean solar spectrum,
79    ! typically between -1.0 and 1.0 and provided by
80    ! single_level%solar_spectral_multiplier.
81    real(jprb), allocatable :: norm_amplitude_solar_irradiance(:)
82   
83    ! Rayleigh molar scattering coefficient in m2 mol-1 in each g
84    ! point
85    real(jprb), allocatable :: rayleigh_molar_scat(:)
86
87    ! ! Spectral mapping of g points
88
89    ! ! Number of wavenumber intervals
90    ! integer :: nwav = 0
91
92    ! Number of k terms / g points
93    integer :: ng   = 0
94
95    ! Spectral definition describing bands and g points
96    type(spectral_definition_type) :: spectral_def
97
98    ! Shortwave: true, longwave: false
99    logical :: is_sw
100
101  contains
102
103    procedure :: read => read_ckd_model
104    procedure :: read_spectral_solar_cycle
105! Vectorized version of the optical depth look-up performs better on
106! NEC, but slower on x86
107#ifdef __SX__
108    procedure :: calc_optical_depth => calc_optical_depth_ckd_model_vec
109#else
110    procedure :: calc_optical_depth => calc_optical_depth_ckd_model
111#endif
112    procedure :: print => print_ckd_model
113    procedure :: calc_planck_function
114    procedure :: calc_incoming_sw
115!    procedure :: deallocate => deallocate_ckd_model
116
117  end type ckd_model_type
118
119
120contains
121
122  !---------------------------------------------------------------------
123  ! Read a complete ecCKD gas optics model from a NetCDF file
124  ! "filename"
125  subroutine read_ckd_model(this, filename, iverbose)
126
127    use easy_netcdf,  only : netcdf_file
128    !use radiation_io, only : nulerr, radiation_abort
129    use yomhook,      only : lhook, dr_hook, jphook
130
131    class(ckd_model_type), intent(inout) :: this
132    character(len=*),      intent(in)    :: filename
133    integer, optional,     intent(in)    :: iverbose
134
135    type(netcdf_file) :: file
136
137    real(jprb), allocatable :: pressure_lut(:)
138    real(jprb), allocatable :: temperature_full(:,:)
139    real(jprb), allocatable :: temperature_planck(:)
140
141    character(len=512) :: constituent_id
142
143    integer :: iverbose_local
144
145    ! Loop counters
146    integer :: jgas, jjgas
147
148    integer :: istart, inext, nchar, i_gas_code
149
150    real(jphook) :: hook_handle
151
152    if (lhook) call dr_hook('radiation_ecckd:read_ckd_model',0,hook_handle)
153
154    if (present(iverbose)) then
155      iverbose_local = iverbose
156    else
157      iverbose_local = 3
158    end if
159
160    call file%open(trim(filename), iverbose=iverbose_local)
161
162    ! Read temperature and pressure coordinate variables
163    call file%get('pressure', pressure_lut)
164    this%log_pressure1 = log(pressure_lut(1))
165    this%npress = size(pressure_lut)
166    this%d_log_pressure = log(pressure_lut(2)) - this%log_pressure1
167    call file%get('temperature', temperature_full)
168    ! AI oct 2023 ajout pour le double appel de ecrad   
169    if (allocated(this%temperature1)) deallocate(this%temperature1)
170    allocate(this%temperature1(this%npress));
171    this%temperature1 = temperature_full(:,1)
172    this%d_temperature = temperature_full(1,2)-temperature_full(1,1)
173    this%ntemp = size(temperature_full,2)
174    deallocate(temperature_full)
175   
176    ! Read Planck function, or solar irradiance and Rayleigh
177    ! scattering coefficient
178    if (file%exists('solar_irradiance')) then
179      this%is_sw = .true.
180      call file%get('solar_irradiance', this%norm_solar_irradiance)
181      this%norm_solar_irradiance = this%norm_solar_irradiance &
182           &  / sum(this%norm_solar_irradiance)
183      call file%get('rayleigh_molar_scattering_coeff', &
184           &  this%rayleigh_molar_scat)
185    else
186      this%is_sw = .false.
187      call file%get('temperature_planck', temperature_planck)
188      this%nplanck = size(temperature_planck)
189      this%temperature1_planck = temperature_planck(1)
190      this%d_temperature_planck = temperature_planck(2) - temperature_planck(1)
191      deallocate(temperature_planck)
192      call file%get('planck_function', this%planck_function)
193    end if
194
195    ! Read the spectral definition information into a separate
196    ! derived type
197    call this%spectral_def%read(file);
198    this%ng = this%spectral_def%ng
199
200    ! Read gases
201    call file%get('n_gases', this%ngas)
202    if (allocated(this%single_gas)) deallocate(this%single_gas)
203    allocate(this%single_gas(this%ngas))
204    call file%get_global_attribute('constituent_id', constituent_id)
205    nchar = len(trim(constituent_id))
206    istart = 1
207    this%i_gas_mapping = 0
208    do jgas = 1, this%ngas
209      if (jgas < this%ngas) then
210        inext = istart + scan(constituent_id(istart:nchar), ' ')
211      else
212        inext = nchar+2
213      end if
214      ! Find gas code
215      i_gas_code = 0
216      do jjgas = 1, NMaxGases
217        if (constituent_id(istart:inext-2) == trim(GasLowerCaseName(jjgas))) then
218          i_gas_code = jjgas
219          exit
220        end if
221      end do
222      ! if (i_gas_code == 0) then
223      !   write(nulerr,'(a,a,a)') '*** Error: Gas "', constituent_id(istart:inext-2), &
224      !        & '" not understood'
225      !   call radiation_abort('Radiation configuration error')
226      ! end if
227      this%i_gas_mapping(i_gas_code) = jgas;
228      call this%single_gas(jgas)%read(file, constituent_id(istart:inext-2), i_gas_code)
229      istart = inext
230    end do
231
232    call file%close()
233
234    if (lhook) call dr_hook('radiation_ecckd:read_ckd_model',1,hook_handle)
235
236  end subroutine read_ckd_model
237
238  !---------------------------------------------------------------------
239  ! Print a description of the correlated k-distribution model to the
240  ! "nulout" unit
241  subroutine print_ckd_model(this)
242
243    use radiation_io, only : nulout
244    use radiation_gas_constants
245
246    class(ckd_model_type), intent(in)  :: this
247
248    integer :: jgas
249   
250    if (this%is_sw) then
251      write(nulout,'(a)',advance='no') 'ecCKD shortwave gas optics model: '
252    else
253      write(nulout,'(a)',advance='no') 'ecCKD longwave gas optics model: '
254    end if
255
256    write(nulout,'(i0,a,i0,a,i0,a,i0,a)') &
257         &  nint(this%spectral_def%wavenumber1(1)), '-', &
258         &  nint(this%spectral_def%wavenumber2(size(this%spectral_def%wavenumber2))), &
259         &  ' cm-1, ', this%ng, ' g-points in ', this%spectral_def%nband, ' bands'
260    write(nulout,'(a,i0,a,i0,a,i0,a)') '  Look-up table sizes: ', this%npress, ' pressures, ', &
261         &  this%ntemp, ' temperatures, ', this%nplanck, ' planck-function entries'
262    write(nulout, '(a)') '  Gases:'
263    do jgas = 1,this%ngas
264      if (this%single_gas(jgas)%i_gas_code > 0) then
265        write(nulout, '(a,a,a)', advance='no') '    ', &
266             &  trim(GasName(this%single_gas(jgas)%i_gas_code)), ': '
267      else
268        write(nulout, '(a)', advance='no') '    Composite of well-mixed background gases: '
269      end if
270      select case (this%single_gas(jgas)%i_conc_dependence)
271        case (IConcDependenceNone)
272          write(nulout, '(a)') 'no concentration dependence'
273        case (IConcDependenceLinear)
274          write(nulout, '(a)') 'linear concentration dependence'
275        case (IConcDependenceRelativeLinear)
276          write(nulout, '(a,e14.6)') 'linear concentration dependence relative to a mole fraction of ', &
277               &  this%single_gas(jgas)%reference_mole_frac
278        case (IConcDependenceLUT)
279          write(nulout, '(a,i0,a,e14.6,a,e13.6)') 'look-up table with ', this%single_gas(jgas)%n_mole_frac, &
280               &  ' log-spaced mole fractions in range ', exp(this%single_gas(jgas)%log_mole_frac1), &
281               &  ' to ', exp(this%single_gas(jgas)%log_mole_frac1 &
282               &           + this%single_gas(jgas)%n_mole_frac*this%single_gas(jgas)%d_log_mole_frac)
283      end select
284    end do
285
286  end subroutine print_ckd_model
287
288
289  !---------------------------------------------------------------------
290  ! Read the amplitude of the spectral variations associated with the
291  ! solar cycle and map to g-points
292  subroutine read_spectral_solar_cycle(this, filename, iverbose, use_updated_solar_spectrum)
293
294    use easy_netcdf,  only : netcdf_file
295    use radiation_io, only : nulout, nulerr, radiation_abort
296    use yomhook,      only : lhook, dr_hook, jphook
297
298    ! Reference total solar irradiance (W m-2)
299    real(jprb), parameter :: ReferenceTSI = 1361.0_jprb
300   
301    class(ckd_model_type), intent(inout) :: this
302    character(len=*),      intent(in)    :: filename
303    integer, optional,     intent(in)    :: iverbose
304    ! Do we update the mean solar spectral irradiance for each g-point
305    ! based on the contents of the file?
306    logical, optional,     intent(in)    :: use_updated_solar_spectrum
307
308    type(netcdf_file) :: file
309
310    ! Solar spectral irradiance, its amplitude and wavenumber
311    ! coordinate variable, read from NetCDF file
312    real(jprb), allocatable :: wavenumber(:) ! cm-1
313    real(jprb), allocatable :: ssi(:) ! W m-2 cm
314    real(jprb), allocatable :: ssi_amplitude(:) ! W m-2 cm
315
316    ! As above but on the wavenumber grid delimited by
317    ! this%wavenumber1 and this%wavenumber2
318    real(jprb), allocatable :: ssi_grid(:)
319    real(jprb), allocatable :: ssi_amplitude_grid(:)
320    real(jprb), allocatable :: wavenumber_grid(:)
321   
322    ! Old normalized solar irradiance in case it gets changed and we
323    ! need to report the amplitude of the change
324    real(jprb), allocatable :: old_norm_solar_irradiance(:)
325   
326    real(jprb) :: dwav_grid
327   
328    ! Number of input wavenumbers, and number on ecCKD model's grid
329    integer :: nwav, nwav_grid
330    ! Corresponding loop indices
331    integer :: jwav, jwav_grid, jg
332
333    integer :: iband
334   
335    integer :: iverbose_local
336
337    real(jphook) :: hook_handle
338   
339    if (lhook) call dr_hook('radiation_ecckd:read_spectral_solar_cycle',0,hook_handle)
340
341    if (present(iverbose)) then
342      iverbose_local = iverbose
343    else
344      iverbose_local = 3
345    end if
346
347    call file%open(trim(filename), iverbose=iverbose_local)
348
349    call file%get('wavenumber', wavenumber)
350    call file%get('mean_solar_spectral_irradiance', ssi)
351    call file%get('ssi_solar_cycle_amplitude', ssi_amplitude)
352
353    call file%close()
354
355    nwav = size(wavenumber)
356   
357    nwav_grid = size(this%spectral_def%wavenumber1)
358    allocate(ssi_grid(nwav_grid))
359    allocate(ssi_amplitude_grid(nwav_grid))
360    allocate(wavenumber_grid(nwav_grid))
361    wavenumber_grid = 0.5_jprb * (this%spectral_def%wavenumber1+this%spectral_def%wavenumber2)
362    dwav_grid = this%spectral_def%wavenumber2(1)-this%spectral_def%wavenumber1(1)
363   
364    ssi_grid = 0.0_jprb
365    ssi_amplitude_grid = 0.0_jprb
366
367    ! Interpolate input SSI to regular wavenumber grid
368    do jwav_grid = 1,nwav_grid
369      do jwav = 1,nwav-1
370        if (wavenumber(jwav) < wavenumber_grid(jwav_grid) &
371             &  .and. wavenumber(jwav+1) >= wavenumber_grid(jwav_grid)) then
372          ! Linear interpolation - this is not perfect
373          ssi_grid(jwav_grid) = (ssi(jwav)*(wavenumber(jwav+1)-wavenumber_grid(jwav_grid)) &
374               &                +ssi(jwav+1)*(wavenumber_grid(jwav_grid)-wavenumber(jwav))) &
375               &                * dwav_grid / (wavenumber(jwav+1)-wavenumber(jwav))
376          ssi_amplitude_grid(jwav_grid) = (ssi_amplitude(jwav)*(wavenumber(jwav+1)-wavenumber_grid(jwav_grid)) &
377               &                +ssi_amplitude(jwav+1)*(wavenumber_grid(jwav_grid)-wavenumber(jwav))) &
378               &                * dwav_grid / (wavenumber(jwav+1)-wavenumber(jwav))
379          exit
380        end if
381      end do
382    end do
383
384    ! Optionally update the solar irradiances in each g-point, and the
385    ! spectral solar irradiance on the wavenumber grid corresponding
386    ! to gpoint_fraction
387    allocate(old_norm_solar_irradiance(nwav_grid))
388    old_norm_solar_irradiance = this%norm_solar_irradiance
389    if (present(use_updated_solar_spectrum)) then
390      if (use_updated_solar_spectrum) then
391        if (.not. allocated(this%spectral_def%solar_spectral_irradiance)) then
392          write(nulerr,'(a)') 'Cannot use_updated_solar_spectrum unless gas optics model is from ecCKD >= 1.4'
393          call radiation_abort()
394        end if
395        this%norm_solar_irradiance = old_norm_solar_irradiance &
396             &  * matmul(ssi_grid,this%spectral_def%gpoint_fraction) &
397             &  / matmul(this%spectral_def%solar_spectral_irradiance,this%spectral_def%gpoint_fraction)
398        this%norm_solar_irradiance = this%norm_solar_irradiance / sum(this%norm_solar_irradiance)
399        this%spectral_def%solar_spectral_irradiance = ssi_grid
400      end if
401    end if
402   
403    ! Map on to g-points
404    this%norm_amplitude_solar_irradiance &
405         &  = this%norm_solar_irradiance &
406         &  * matmul(ssi_amplitude_grid, this%spectral_def%gpoint_fraction) &
407         &  / matmul(ssi_grid,this%spectral_def%gpoint_fraction)
408   
409    ! Remove the mean from the solar-cycle fluctuations, since the
410    ! user will scale with total solar irradiance
411    this%norm_amplitude_solar_irradiance &
412         &  = (this%norm_solar_irradiance+this%norm_amplitude_solar_irradiance) &
413         &  / sum(this%norm_solar_irradiance+this%norm_amplitude_solar_irradiance) &
414         &  - this%norm_solar_irradiance
415
416    ! Print the spectral solar irradiance per g point, and solar cycle amplitude
417    if (iverbose_local >= 2) then
418      write(nulout,'(a,f6.1,a)') 'G-point, solar irradiance for nominal TSI = ', &
419           &  ReferenceTSI, ' W m-2, solar cycle amplitude (at solar maximum), update to original solar irradiance'
420      iband = 0
421      do jg = 1,this%ng
422        if (this%spectral_def%i_band_number(jg) > iband) then
423          iband = this%spectral_def%i_band_number(jg)
424          write(nulout, '(i2,f10.4,f7.3,a,f8.4,a,i2,a,f7.1,a,f7.1,a)') &
425               &  jg, ReferenceTSI*this%norm_solar_irradiance(jg), &
426               &  100.0_jprb * this%norm_amplitude_solar_irradiance(jg) &
427               &  / this%norm_solar_irradiance(jg), '% ', &
428               &  100.0_jprb * (this%norm_solar_irradiance(jg) &
429               &  / old_norm_solar_irradiance(jg) - 1.0_jprb), '% Band ', iband, ': ', &
430               &  this%spectral_def%wavenumber1_band(iband), '-', &
431               &  this%spectral_def%wavenumber2_band(iband), ' cm-1'
432        else
433          write(nulout, '(i2,f10.4,f7.3,a,f8.4,a)') jg, ReferenceTSI*this%norm_solar_irradiance(jg), &
434               &  100.0_jprb * this%norm_amplitude_solar_irradiance(jg) &
435               &  / this%norm_solar_irradiance(jg), '% ', &
436               &  100.0_jprb * (this%norm_solar_irradiance(jg) &
437               &  / old_norm_solar_irradiance(jg) - 1.0_jprb), '%'
438        end if
439      end do
440    end if
441   
442    if (lhook) call dr_hook('radiation_ecckd:read_spectral_solar_cycle',1,hook_handle)
443   
444  end subroutine read_spectral_solar_cycle
445 
446
447  !---------------------------------------------------------------------
448  ! Compute layerwise optical depth for each g point for ncol columns
449  ! at nlev layers
450  subroutine calc_optical_depth_ckd_model(this, ncol, nlev, istartcol, iendcol, nmaxgas, &
451       &  pressure_hl, temperature_fl, mole_fraction_fl, &
452       &  optical_depth_fl, rayleigh_od_fl)
453
454    use yomhook,             only : lhook, dr_hook, jphook
455    use radiation_constants, only : AccelDueToGravity
456
457    ! Input variables
458
459    class(ckd_model_type), intent(in), target  :: this
460    ! Number of columns, levels and input gases
461    integer,               intent(in)  :: ncol, nlev, nmaxgas, istartcol, iendcol
462    ! Pressure at half levels (Pa), dimensioned (ncol,nlev+1)
463    real(jprb),            intent(in)  :: pressure_hl(ncol,nlev+1)
464    ! Temperature at full levels (K), dimensioned (ncol,nlev)
465    real(jprb),            intent(in)  :: temperature_fl(istartcol:iendcol,nlev)
466    ! Gas mole fractions at full levels (mol mol-1), dimensioned (ncol,nlev,nmaxgas)
467    real(jprb),            intent(in)  :: mole_fraction_fl(ncol,nlev,nmaxgas)
468   
469    ! Output variables
470
471    ! Layer absorption optical depth for each g point
472    real(jprb),            intent(out) :: optical_depth_fl(this%ng,nlev,istartcol:iendcol)
473    ! In the shortwave only, the Rayleigh scattering optical depth
474    real(jprb),  optional, intent(out) :: rayleigh_od_fl(this%ng,nlev,istartcol:iendcol)
475
476    ! Local variables
477
478    real(jprb), pointer :: molar_abs(:,:,:), molar_abs_conc(:,:,:,:)
479
480    ! Natural logarithm of pressure at full levels
481    real(jprb) :: log_pressure_fl(nlev)
482
483    ! Optical depth of single gas at one point in space versus
484    ! spectral interval
485    !real(jprb) :: od_single_gas(this%ng)
486
487    real(jprb) :: multiplier(nlev), simple_multiplier(nlev), global_multiplier, temperature1
488
489    ! Indices and weights in temperature, pressure and concentration interpolation
490    real(jprb) :: pindex1, tindex1, cindex1
491    real(jprb) :: pw1(nlev), pw2(nlev), tw1(nlev), tw2(nlev), cw1(nlev), cw2(nlev)
492    integer    :: ip1(nlev), it1(nlev), ic1(nlev)
493
494    ! Natural logarithm of mole fraction at one point
495    real(jprb) :: log_conc
496
497    ! Minimum mole fraction in look-up-table
498    real(jprb) :: mole_frac1
499
500    integer :: jcol, jlev, jgas, igascode
501
502    real(jphook) :: hook_handle
503
504    if (lhook) call dr_hook('radiation_ecckd:calc_optical_depth',0,hook_handle)
505
506    global_multiplier = 1.0_jprb / (AccelDueToGravity * 0.001_jprb * AirMolarMass)
507
508    do jcol = istartcol,iendcol
509
510      log_pressure_fl = log(0.5_jprb * (pressure_hl(jcol,1:nlev)+pressure_hl(jcol,2:nlev+1)))
511
512      do jlev = 1,nlev
513        ! Find interpolation points in pressure
514        pindex1 = (log_pressure_fl(jlev)-this%log_pressure1) &
515             &    / this%d_log_pressure
516        pindex1 = 1.0_jprb + max(0.0_jprb, min(pindex1, this%npress-1.0001_jprb))
517        ip1(jlev) = int(pindex1)
518        pw2(jlev) = pindex1 - ip1(jlev)
519        pw1(jlev) = 1.0_jprb - pw2(jlev)
520
521        ! Find interpolation points in temperature
522        temperature1 = pw1(jlev)*this%temperature1(ip1(jlev)) &
523             &       + pw2(jlev)*this%temperature1(ip1(jlev)+1)
524        tindex1 = (temperature_fl(jcol,jlev) - temperature1) &
525             &    / this%d_temperature
526        tindex1 = 1.0_jprb + max(0.0_jprb, min(tindex1, this%ntemp-1.0001_jprb))
527        it1(jlev) = int(tindex1)
528        tw2(jlev) = tindex1 - it1(jlev)
529        tw1(jlev) = 1.0_jprb - tw2(jlev)
530
531        ! Concentration multiplier
532        simple_multiplier(jlev) = global_multiplier &
533             &  * (pressure_hl(jcol,jlev+1) - pressure_hl(jcol,jlev))
534      end do
535
536      optical_depth_fl(:,:,jcol) = 0.0_jprb
537     
538      do jgas = 1,this%ngas
539
540        associate (single_gas => this%single_gas(jgas))
541          igascode = this%single_gas(jgas)%i_gas_code
542         
543          select case (single_gas%i_conc_dependence)
544           
545          case (IConcDependenceLinear)
546            molar_abs => this%single_gas(jgas)%molar_abs
547            multiplier = simple_multiplier * mole_fraction_fl(jcol,:,igascode)
548
549            do jlev = 1,nlev
550              optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) &
551                   &        + (multiplier(jlev)*tw1(jlev)) * (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)) &
552                   &                +pw2(jlev) * molar_abs(:,ip1(jlev)+1,it1(jlev))) &
553                   &        + (multiplier(jlev)*tw2(jlev)) * (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)+1) &
554                   &                +pw2(jlev) * molar_abs(:,ip1(jlev)+1,it1(jlev)+1))
555            end do
556
557          case (IConcDependenceRelativeLinear)
558            molar_abs => this%single_gas(jgas)%molar_abs
559            multiplier = simple_multiplier  * (mole_fraction_fl(jcol,:,igascode) &
560                 &                            - single_gas%reference_mole_frac)
561            do jlev = 1,nlev
562              optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) &
563                   &        + (multiplier(jlev)*tw1(jlev)) * (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)) &
564                   &                +pw2(jlev) * molar_abs(:,ip1(jlev)+1,it1(jlev))) &
565                   &        + (multiplier(jlev)*tw2(jlev)) * (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)+1) &
566                   &                +pw2(jlev) * molar_abs(:,ip1(jlev)+1,it1(jlev)+1))
567            end do
568
569          case (IConcDependenceNone)
570            ! Composite gases
571            molar_abs => this%single_gas(jgas)%molar_abs
572            do jlev = 1,nlev
573              optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) &
574                   &  + (simple_multiplier(jlev)*tw1(jlev)) * (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)) &
575                   &                              +pw2(jlev) * molar_abs(:,ip1(jlev)+1,it1(jlev))) &
576                   &  + (simple_multiplier(jlev)*tw2(jlev)) * (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)+1) &
577                   &                              +pw2(jlev) * molar_abs(:,ip1(jlev)+1,it1(jlev)+1))
578            end do
579
580          case (IConcDependenceLUT)
581            ! Logarithmic interpolation in concentration space
582            molar_abs_conc => this%single_gas(jgas)%molar_abs_conc
583            mole_frac1 = exp(single_gas%log_mole_frac1)
584            do jlev = 1,nlev
585              ! Take care of mole_fraction == 0
586              log_conc = log(max(mole_fraction_fl(jcol,jlev,igascode), mole_frac1))
587              cindex1  = (log_conc - single_gas%log_mole_frac1) / single_gas%d_log_mole_frac
588              cindex1  = 1.0_jprb + max(0.0_jprb, min(cindex1, single_gas%n_mole_frac-1.0001_jprb))
589              ic1(jlev) = int(cindex1)
590              cw2(jlev) = cindex1 - ic1(jlev)
591              cw1(jlev) = 1.0_jprb - cw2(jlev)
592            end do
593              ! od_single_gas = cw1 * (tw1 * (pw1 * molar_abs_conc(:,ip1,it1,ic1) &
594              !      &                       +pw2 * molar_abs_conc(:,ip1+1,it1,ic1)) &
595              !      &                +tw2 * (pw1 * molar_abs_conc(:,ip1,it1+1,ic1) &
596              !      &                       +pw2 * molar_abs_conc(:,ip1+1,it1+1,ic1))) &
597              !      &         +cw2 * (tw1 * (pw1 * molar_abs_conc(:,ip1,it1,ic1+1) &
598              !      &                       +pw2 * molar_abs_conc(:,ip1+1,it1,ic1+1)) &
599              !      &                +tw2 * (pw1 * molar_abs_conc(:,ip1,it1+1,ic1+1) &
600              !      &                       +pw2 * molar_abs_conc(:,ip1+1,it1+1,ic1+1)))
601            do jlev = 1,nlev
602              optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) &
603                   &  + (simple_multiplier(jlev) * mole_fraction_fl(jcol,jlev,igascode)) * ( &
604                   &      (cw1(jlev) * tw1(jlev) * pw1(jlev)) * molar_abs_conc(:,ip1(jlev),it1(jlev),ic1(jlev)) &
605                   &     +(cw1(jlev) * tw1(jlev) * pw2(jlev)) * molar_abs_conc(:,ip1(jlev)+1,it1(jlev),ic1(jlev)) &
606                   &     +(cw1(jlev) * tw2(jlev) * pw1(jlev)) * molar_abs_conc(:,ip1(jlev),it1(jlev)+1,ic1(jlev)) &
607                   &     +(cw1(jlev) * tw2(jlev) * pw2(jlev)) * molar_abs_conc(:,ip1(jlev)+1,it1(jlev)+1,ic1(jlev)) &
608                   &     +(cw2(jlev) * tw1(jlev) * pw1(jlev)) * molar_abs_conc(:,ip1(jlev),it1(jlev),ic1(jlev)+1) &
609                   &     +(cw2(jlev) * tw1(jlev) * pw2(jlev)) * molar_abs_conc(:,ip1(jlev)+1,it1(jlev),ic1(jlev)+1) &
610                   &     +(cw2(jlev) * tw2(jlev) * pw1(jlev)) * molar_abs_conc(:,ip1(jlev),it1(jlev)+1,ic1(jlev)+1) &
611                   &     +(cw2(jlev) * tw2(jlev) * pw2(jlev)) * molar_abs_conc(:,ip1(jlev)+1,it1(jlev)+1,ic1(jlev)+1))
612            end do
613          end select
614           
615        end associate
616
617      end do
618
619      ! Ensure the optical depth is not negative
620      optical_depth_fl(:,:,jcol) = max(0.0_jprb, optical_depth_fl(:,:,jcol))
621
622      ! Rayleigh scattering
623      if (this%is_sw .and. present(rayleigh_od_fl)) then
624        do jlev = 1,nlev
625          rayleigh_od_fl(:,jlev,jcol) = global_multiplier &
626               &  * (pressure_hl(jcol,jlev+1) - pressure_hl(jcol,jlev)) * this%rayleigh_molar_scat
627        end do
628      end if
629
630    end do
631
632    if (lhook) call dr_hook('radiation_ecckd:calc_optical_depth',1,hook_handle)
633
634  end subroutine calc_optical_depth_ckd_model
635
636
637  !---------------------------------------------------------------------
638  ! Vectorized variant of above routine
639  subroutine calc_optical_depth_ckd_model_vec(this, ncol, nlev, istartcol, iendcol, nmaxgas, &
640       &  pressure_hl, temperature_fl, mole_fraction_fl, &
641       &  optical_depth_fl, rayleigh_od_fl)
642
643    use yomhook,             only : lhook, dr_hook, jphook
644    use radiation_constants, only : AccelDueToGravity
645
646    ! Input variables
647
648    class(ckd_model_type), intent(in), target  :: this
649    ! Number of columns, levels and input gases
650    integer,               intent(in)  :: ncol, nlev, nmaxgas, istartcol, iendcol
651    ! Pressure at half levels (Pa), dimensioned (ncol,nlev+1)
652    real(jprb),            intent(in)  :: pressure_hl(ncol,nlev+1)
653    ! Temperature at full levels (K), dimensioned (ncol,nlev)
654    real(jprb),            intent(in)  :: temperature_fl(istartcol:iendcol,nlev)
655    ! Gas mole fractions at full levels (mol mol-1), dimensioned (ncol,nlev,nmaxgas)
656    real(jprb),            intent(in)  :: mole_fraction_fl(ncol,nlev,nmaxgas)
657   
658    ! Output variables
659
660    ! Layer absorption optical depth for each g point
661    real(jprb),            intent(out) :: optical_depth_fl(this%ng,nlev,istartcol:iendcol)
662    ! In the shortwave only, the Rayleigh scattering optical depth
663    real(jprb),  optional, intent(out) :: rayleigh_od_fl(this%ng,nlev,istartcol:iendcol)
664
665    ! Local variables
666
667    real(jprb), pointer :: molar_abs(:,:,:), molar_abs_conc(:,:,:,:)
668
669    ! Natural logarithm of pressure at full levels
670    real(jprb) :: log_pressure_fl
671
672    ! Optical depth of single gas at one point in space versus
673    ! spectral interval
674    !real(jprb) :: od_single_gas(this%ng)
675
676    real(jprb) :: multiplier, simple_multiplier(ncol,nlev), global_multiplier, temperature1
677
678    ! Indices and weights in temperature, pressure and concentration interpolation
679    real(jprb) :: pindex1, tindex1, cindex1
680    real(jprb) :: pw1(ncol,nlev), pw2(ncol,nlev), tw1(ncol,nlev), tw2(ncol,nlev), cw1(ncol,nlev), cw2(ncol,nlev)
681    integer    :: ip1(ncol,nlev), it1(ncol,nlev), ic1(ncol,nlev)
682
683    ! Natural logarithm of mole fraction at one point
684    real(jprb) :: log_conc
685
686    ! Minimum mole fraction in look-up-table
687    real(jprb) :: mole_frac1
688
689    ! Layer absorption optical depth for each g point (memory layout adjusted to vectorization)
690    real(jprb) :: od_fl(ncol,this%ng,nlev)
691
692    integer :: jcol, jlev, jgas, igascode, jg
693
694    real(jphook) :: hook_handle
695
696    if (lhook) call dr_hook('radiation_ecckd:calc_optical_depth_vec',0,hook_handle)
697
698    global_multiplier = 1.0_jprb / (AccelDueToGravity * 0.001_jprb * AirMolarMass)
699
700    od_fl(:,:,:) = 0.0_jprb
701
702    do jlev = 1,nlev
703      do jcol = istartcol,iendcol
704
705        log_pressure_fl = log(0.5_jprb * (pressure_hl(jcol,jlev)+pressure_hl(jcol,jlev+1)))
706
707        ! Find interpolation points in pressure
708        pindex1 = (log_pressure_fl-this%log_pressure1) &
709             &    / this%d_log_pressure
710        pindex1 = 1.0_jprb + max(0.0_jprb, min(pindex1, this%npress-1.0001_jprb))
711        ip1(jcol,jlev) = int(pindex1)
712        pw2(jcol,jlev) = pindex1 - ip1(jcol,jlev)
713        pw1(jcol,jlev) = 1.0_jprb - pw2(jcol,jlev)
714
715        ! Find interpolation points in temperature
716        temperature1 = pw1(jcol,jlev)*this%temperature1(ip1(jcol,jlev)) &
717             &       + pw2(jcol,jlev)*this%temperature1(ip1(jcol,jlev)+1)
718        tindex1 = (temperature_fl(jcol,jlev) - temperature1) &
719             &    / this%d_temperature
720        tindex1 = 1.0_jprb + max(0.0_jprb, min(tindex1, this%ntemp-1.0001_jprb))
721        it1(jcol,jlev) = int(tindex1)
722        tw2(jcol,jlev) = tindex1 - it1(jcol,jlev)
723        tw1(jcol,jlev) = 1.0_jprb - tw2(jcol,jlev)
724
725        ! Concentration multiplier
726        simple_multiplier(jcol,jlev) = global_multiplier &
727             &  * (pressure_hl(jcol,jlev+1) - pressure_hl(jcol,jlev))
728      end do
729    end do
730
731    do jgas = 1,this%ngas
732
733      associate (single_gas => this%single_gas(jgas))
734        igascode = this%single_gas(jgas)%i_gas_code
735         
736        select case (single_gas%i_conc_dependence)
737           
738        case (IConcDependenceLinear)
739          molar_abs => this%single_gas(jgas)%molar_abs
740
741          do jlev = 1,nlev
742            do jg = 1, this%ng
743              do jcol = istartcol,iendcol
744                multiplier = simple_multiplier(jcol,jlev) * mole_fraction_fl(jcol,jlev,igascode)
745
746                od_fl(jcol,jg,jlev) = od_fl(jcol,jg,jlev) &
747                   &        + (multiplier*tw1(jcol,jlev)) * (pw1(jcol,jlev) * molar_abs(jg,ip1(jcol,jlev),it1(jcol,jlev)) &
748                   &                +pw2(jcol,jlev) * molar_abs(jg,ip1(jcol,jlev)+1,it1(jcol,jlev))) &
749                   &        + (multiplier*tw2(jcol,jlev)) * (pw1(jcol,jlev) * molar_abs(jg,ip1(jcol,jlev),it1(jcol,jlev)+1) &
750                   &                +pw2(jcol,jlev) * molar_abs(jg,ip1(jcol,jlev)+1,it1(jcol,jlev)+1))
751              end do
752            end do
753          end do
754
755        case (IConcDependenceRelativeLinear)
756          molar_abs => this%single_gas(jgas)%molar_abs
757
758          do jlev = 1,nlev
759            do jg = 1, this%ng
760              do jcol = istartcol,iendcol
761                multiplier = simple_multiplier(jcol,jlev)  * (mole_fraction_fl(jcol,jlev,igascode) &
762                   &                            - single_gas%reference_mole_frac)
763
764                od_fl(jcol,jg,jlev) = od_fl(jcol,jg,jlev) &
765                   &        + (multiplier*tw1(jcol,jlev)) * (pw1(jcol,jlev) * molar_abs(jg,ip1(jcol,jlev),it1(jcol,jlev)) &
766                   &                +pw2(jcol,jlev) * molar_abs(jg,ip1(jcol,jlev)+1,it1(jcol,jlev))) &
767                   &        + (multiplier*tw2(jcol,jlev)) * (pw1(jcol,jlev) * molar_abs(jg,ip1(jcol,jlev),it1(jcol,jlev)+1) &
768                   &                +pw2(jcol,jlev) * molar_abs(jg,ip1(jcol,jlev)+1,it1(jcol,jlev)+1))
769              end do
770            end do
771          end do
772
773        case (IConcDependenceNone)
774          ! Composite gases
775          molar_abs => this%single_gas(jgas)%molar_abs
776
777          do jlev = 1,nlev
778            do jg = 1, this%ng
779              do jcol = istartcol,iendcol
780                od_fl(jcol,jg,jlev) = od_fl(jcol,jg,jlev) &
781                   &        + (simple_multiplier(jcol,jlev)*tw1(jcol,jlev)) *   &
782                   &                (pw1(jcol,jlev) * molar_abs(jg,ip1(jcol,jlev),it1(jcol,jlev)) &
783                   &                +pw2(jcol,jlev) * molar_abs(jg,ip1(jcol,jlev)+1,it1(jcol,jlev))) &
784                   &        + (simple_multiplier(jcol,jlev)*tw2(jcol,jlev)) *   &
785                   &                (pw1(jcol,jlev) * molar_abs(jg,ip1(jcol,jlev),it1(jcol,jlev)+1) &
786                   &                +pw2(jcol,jlev) * molar_abs(jg,ip1(jcol,jlev)+1,it1(jcol,jlev)+1))
787              end do
788            end do
789          end do
790
791          case (IConcDependenceLUT)
792            ! Logarithmic interpolation in concentration space
793            molar_abs_conc => this%single_gas(jgas)%molar_abs_conc
794            mole_frac1 = exp(single_gas%log_mole_frac1)
795
796            do jlev = 1,nlev
797              do jcol = istartcol,iendcol
798                ! Take care of mole_fraction == 0
799                log_conc = log(max(mole_fraction_fl(jcol,jlev,igascode), mole_frac1))
800                cindex1  = (log_conc - single_gas%log_mole_frac1) / single_gas%d_log_mole_frac
801                cindex1  = 1.0_jprb + max(0.0_jprb, min(cindex1, single_gas%n_mole_frac-1.0001_jprb))
802                ic1(jcol,jlev) = int(cindex1)
803                cw2(jcol,jlev) = cindex1 - ic1(jcol,jlev)
804                cw1(jcol,jlev) = 1.0_jprb - cw2(jcol,jlev)
805              end do
806            end do
807
808          do jlev = 1,nlev
809            do jg = 1, this%ng
810!NEC$ select_vector
811              do jcol = istartcol,iendcol
812
813                od_fl(jcol,jg,jlev) = od_fl(jcol,jg,jlev) &
814                   &  + (simple_multiplier(jcol,jlev) * mole_fraction_fl(jcol,jlev,igascode)) * ( &
815                   &      (cw1(jcol,jlev) * tw1(jcol,jlev) * pw1(jcol,jlev)) * &
816                   &      molar_abs_conc(jg,ip1(jcol,jlev),it1(jcol,jlev),ic1(jcol,jlev)) &
817                   &     +(cw1(jcol,jlev) * tw1(jcol,jlev) * pw2(jcol,jlev)) * &
818                   &      molar_abs_conc(jg,ip1(jcol,jlev)+1,it1(jcol,jlev),ic1(jcol,jlev)) &
819                   &     +(cw1(jcol,jlev) * tw2(jcol,jlev) * pw1(jcol,jlev)) * &
820                   &      molar_abs_conc(jg,ip1(jcol,jlev),it1(jcol,jlev)+1,ic1(jcol,jlev)) &
821                   &     +(cw1(jcol,jlev) * tw2(jcol,jlev) * pw2(jcol,jlev)) * &
822                   &      molar_abs_conc(jg,ip1(jcol,jlev)+1,it1(jcol,jlev)+1,ic1(jcol,jlev)) &
823                   &     +(cw2(jcol,jlev) * tw1(jcol,jlev) * pw1(jcol,jlev)) * &
824                   &      molar_abs_conc(jg,ip1(jcol,jlev),it1(jcol,jlev),ic1(jcol,jlev)+1) &
825                   &     +(cw2(jcol,jlev) * tw1(jcol,jlev) * pw2(jcol,jlev)) * &
826                   &      molar_abs_conc(jg,ip1(jcol,jlev)+1,it1(jcol,jlev),ic1(jcol,jlev)+1) &
827                   &     +(cw2(jcol,jlev) * tw2(jcol,jlev) * pw1(jcol,jlev)) * &
828                   &      molar_abs_conc(jg,ip1(jcol,jlev),it1(jcol,jlev)+1,ic1(jcol,jlev)+1) &
829                   &     +(cw2(jcol,jlev) * tw2(jcol,jlev) * pw2(jcol,jlev)) * &
830                   &      molar_abs_conc(jg,ip1(jcol,jlev)+1,it1(jcol,jlev)+1,ic1(jcol,jlev)+1))
831              end do
832            end do
833          end do
834        end select
835           
836      end associate
837
838      ! Ensure the optical depth is not negative
839      do jcol = istartcol,iendcol
840        do jlev = 1,nlev
841          do jg = 1, this%ng
842            optical_depth_fl(jg,jlev,jcol) = max(0.0_jprb, od_fl(jcol,jg,jlev))
843          end do
844        end do
845      end do
846
847      ! Rayleigh scattering
848      if (this%is_sw .and. present(rayleigh_od_fl)) then
849        do jcol = istartcol,iendcol
850          do jlev = 1,nlev
851            do jg = 1, this%ng
852              rayleigh_od_fl(jg,jlev,jcol) = global_multiplier &
853               &  * (pressure_hl(jcol,jlev+1) - pressure_hl(jcol,jlev)) * this%rayleigh_molar_scat(jg)
854            end do
855          end do
856        end do
857      end if
858
859    end do
860
861    if (lhook) call dr_hook('radiation_ecckd:calc_optical_depth_vec',1,hook_handle)
862
863  end subroutine calc_optical_depth_ckd_model_vec
864
865 
866  !---------------------------------------------------------------------
867  ! Calculate the Planck function integrated across each of the g
868  ! points of this correlated k-distribution model, for a given
869  ! temperature, where Planck function is defined as the flux emitted
870  ! by a black body (rather than radiance)
871  subroutine calc_planck_function(this, nt, temperature, planck)
872
873    class(ckd_model_type), intent(in)  :: this
874    integer,    intent(in)  :: nt
875    real(jprb), intent(in)  :: temperature(:) ! K
876    real(jprb), intent(out) :: planck(this%ng,nt) ! W m-2
877
878    real(jprb) :: tindex1, tw1, tw2
879    integer    :: it1, jt
880
881    do jt = 1,nt
882      tindex1 = (temperature(jt) - this%temperature1_planck) &
883           &   * (1.0_jprb / this%d_temperature_planck)
884      if (tindex1 >= 0) then
885        ! Normal interpolation, and extrapolation for high temperatures
886        tindex1 = 1.0_jprb + tindex1
887        it1 = min(int(tindex1), this%nplanck-1)
888        tw2 = tindex1 - it1
889        tw1 = 1.0_jprb - tw2
890        planck(:,jt) = tw1 * this%planck_function(:,it1) &
891             &       + tw2 * this%planck_function(:,it1+1)
892      else
893        ! Interpolate linearly to zero
894        planck(:,jt) = this%planck_function(:,1) &
895             &       * (temperature(jt)/this%temperature1_planck)
896      end if
897    end do
898
899  end subroutine calc_planck_function
900 
901
902  !---------------------------------------------------------------------
903  ! Return the spectral solar irradiance integrated over each g point
904  ! of a solar correlated k-distribution model, given the
905  ! total_solar_irradiance
906  subroutine calc_incoming_sw(this, total_solar_irradiance, &
907       &                      spectral_solar_irradiance, &
908       &                      solar_spectral_multiplier)
909
910    use radiation_io, only : nulerr, radiation_abort
911
912    class(ckd_model_type), intent(in)    :: this
913    real(jprb),            intent(in)    :: total_solar_irradiance ! W m-2
914    real(jprb),            intent(inout) :: spectral_solar_irradiance(:,:) ! W m-2
915    real(jprb), optional,  intent(in)    :: solar_spectral_multiplier
916
917    if (.not. present(solar_spectral_multiplier)) then
918      spectral_solar_irradiance &
919           &  = spread(total_solar_irradiance * this%norm_solar_irradiance, &
920           &           2, size(spectral_solar_irradiance,2))
921    else if (allocated(this%norm_amplitude_solar_irradiance)) then
922      spectral_solar_irradiance &
923           &  = spread(total_solar_irradiance * (this%norm_solar_irradiance &
924           &   + solar_spectral_multiplier*this%norm_amplitude_solar_irradiance), &
925           &           2, size(spectral_solar_irradiance,2))
926    else if (solar_spectral_multiplier == 0.0_jprb) then
927      spectral_solar_irradiance &
928           &  = spread(total_solar_irradiance * this%norm_solar_irradiance, &
929           &           2, size(spectral_solar_irradiance,2))
930    else
931      write(nulerr, '(a)') '*** Error in calc_incoming_sw: no information present on solar cycle'
932      call radiation_abort()
933    end if
934
935  end subroutine calc_incoming_sw
936
937end module radiation_ecckd
Note: See TracBrowser for help on using the repository browser.