source: LMDZ6/trunk/libf/phylmd/ecrad/radiation_aerosol_optics_data.F90 @ 3908

Last change on this file since 3908 was 3908, checked in by idelkadi, 3 years ago

Online implementation of the radiative transfer code ECRAD in the LMDZ model.

  • Inclusion of the ecrad directory containing the sources of the ECRAD code
    • interface routine : radiation_scheme.F90
  • Adaptation of compilation scripts :
    • compilation under CPP key CPP_ECRAD
    • compilation with option "-rad ecard" or "-ecard true"
    • The "-rad old/rtm/ecran" build option will need to replace the "-rrtm true" and "-ecrad true" options in the future.
  • Runing LMDZ simulations with ecrad, you need :
    • logical key iflag_rrtm = 2 in physiq.def
    • namelist_ecrad (DefLists?)
    • the directory "data" containing the configuration files is temporarily placed in ../libfphylmd/ecrad/
  • Compilation and execution are tested in the 1D case. The repository under svn would allow to continue the implementation work: tests, verification of the results, ...
File size: 18.5 KB
Line 
1! radiation_aerosol_optics_data.F90 - Type to store aerosol optical properties
2!
3! (C) Copyright 2015- ECMWF.
4!
5! This software is licensed under the terms of the Apache Licence Version 2.0
6! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
7!
8! In applying this licence, ECMWF does not waive the privileges and immunities
9! granted to it by virtue of its status as an intergovernmental organisation
10! nor does it submit to any jurisdiction.
11!
12! Author:  Robin Hogan
13! Email:   r.j.hogan@ecmwf.int
14!
15! Modifications
16!   2017-10-23  R. Hogan  Renamed single-character variables
17!   2018-04-20  A. Bozzo  Read optical properties at selected wavelengths
18
19
20module radiation_aerosol_optics_data
21
22  use parkind1,      only : jprb
23  use radiation_io,  only : nulerr, radiation_abort
24
25  implicit none
26  public
27
28  private :: get_line
29
30  ! The following provide possible values for
31  ! aerosol_optics_type%iclass, which is used to map the user's type
32  ! index to the hydrophobic or hydrophilic types loaded from the
33  ! aerosol optics file. Initially iclass is equal to
34  ! AerosolClassUndefined, which will throw an error if ever the user
35  ! tries to use this aerosol type. The user may specify that an
36  ! aerosol type is to be ignored in the radiation calculation, in
37  ! which case iclass will be set equal to AerosolClassIgnored.
38  enum, bind(c)
39     enumerator IAerosolClassUndefined,   IAerosolClassIgnored, &
40          &     IAerosolClassHydrophobic, IAerosolClassHydrophilic
41  end enum
42
43  integer, parameter :: NMaxStringLength = 2000
44  integer, parameter :: NMaxLineLength   = 200
45
46  !---------------------------------------------------------------------
47  ! This type holds the configuration information to compute
48  ! aerosol optical properties
49  type aerosol_optics_type
50     ! A vector of length ntype, iclass maps user-defined types on to
51     ! the hydrophilic or hydrophobic aerosol classes using the
52     ! enumerators above
53     integer, allocatable, dimension(:) :: iclass
54
55     ! Also a vector of length ntype, itype maps user-defined types on
56     ! to specific hydrophilic or hydrophobic aerosol types
57     integer, allocatable, dimension(:) :: itype
58
59     ! Scattering properties are provided separately in the shortwave
60     ! and longwave for hydrophobic and hydrophilic aerosols.
61     ! Hydrophobic aerosols are dimensioned (nband,n_type_phobic):
62     real(jprb), allocatable, dimension(:,:) :: &
63          &  mass_ext_sw_phobic, & ! Mass-extinction coefficient (m2 kg-1)
64          &  ssa_sw_phobic,      & ! Single scattering albedo
65          &  g_sw_phobic,        & ! Asymmetry factor
66          &  mass_ext_lw_phobic, & ! Mass-extinction coefficient (m2 kg-1)
67          &  ssa_lw_phobic,      & ! Single scattering albedo
68          &  g_lw_phobic           ! Asymmetry factor
69
70     ! Hydrophilic aerosols are dimensioned (nband, nrh, n_type_philic):
71     real(jprb), allocatable, dimension(:,:,:) :: &
72          &  mass_ext_sw_philic, & ! Mass-extinction coefficient (m2 kg-1)
73          &  ssa_sw_philic,      & ! Single scattering albedo
74          &  g_sw_philic,        & ! Asymmetry factor
75          &  mass_ext_lw_philic, & ! Mass-extinction coefficient (m2 kg-1)
76          &  ssa_lw_philic,      & ! Single scattering albedo
77          &  g_lw_philic           ! Asymmetry factor
78
79     ! Scattering properties at selected wavelengths
80     ! (n_mono_wl,n_type_phobic/philic)
81     real(jprb), allocatable, dimension(:,:) :: &
82          &  mass_ext_mono_phobic, & ! Mass-extinction coefficient (m2 kg-1)
83          &  ssa_mono_phobic,      & ! Single scattering albedo
84          &  g_mono_phobic,        & ! Asymmetry factor
85          &  lidar_ratio_mono_phobic ! Lidar Ratio
86     real(jprb), allocatable, dimension(:,:,:) :: &
87          &  mass_ext_mono_philic, & ! Mass-extinction coefficient (m2 kg-1)
88          &  ssa_mono_philic,      & ! Single scattering albedo
89          &  g_mono_philic,        & ! Asymmetry factor
90          &  lidar_ratio_mono_philic ! Lidar Ratio
91
92     ! For hydrophilic aerosols, the lower bounds of the relative
93     ! humidity bins is a vector of length nrh:
94     real(jprb), allocatable, dimension(:) :: &
95          &  rh_lower    ! Dimensionless (1.0 = 100% humidity)
96
97     ! Strings describing the aerosol types
98     character(len=NMaxStringLength) :: description_phobic_str = ' '
99     character(len=NMaxStringLength) :: description_philic_str = ' '
100
101     ! The number of user-defined aerosol types
102     integer :: ntype
103
104     ! The number of hydrophobic and hydrophilic types read from the
105     ! aerosol optics file
106     integer :: n_type_phobic, n_type_philic
107
108     ! Number of relative humidity bins
109     integer :: nrh
110
111     ! Number of longwave and shortwave bands of the data in the file,
112     ! and monochromatic wavelengths
113     integer :: n_bands_lw = 0, n_bands_sw = 0, n_mono_wl = 0
114
115     ! Do we have any hydrophilic types?
116     logical :: use_hydrophilic = .true.
117
118     ! Do we have monochromatic optical properties
119     logical :: use_monochromatic = .false.
120
121   contains
122     procedure :: setup => setup_aerosol_optics
123     procedure :: set_hydrophobic_type
124     procedure :: set_hydrophilic_type
125     procedure :: set_empty_type
126     procedure :: set_types
127     procedure :: calc_rh_index
128     procedure :: print_description
129
130  end type aerosol_optics_type
131
132contains
133
134
135  !---------------------------------------------------------------------
136  ! Setup aerosol optics coefficients by reading them from a file
137  subroutine setup_aerosol_optics(this, file_name, ntype, iverbose)
138
139    use yomhook,              only : lhook, dr_hook
140    use easy_netcdf,          only : netcdf_file
141    use radiation_io,         only : nulerr, radiation_abort
142
143    class(aerosol_optics_type), intent(inout) :: this
144    character(len=*), intent(in)              :: file_name
145    integer, intent(in)                       :: ntype
146    integer, intent(in), optional             :: iverbose
147
148    ! The NetCDF file containing the aerosol optics data
149    type(netcdf_file)  :: file
150    integer            :: iverb
151    real(jprb)         :: hook_handle
152
153    if (lhook) call dr_hook('radiation_aerosol_optics_data:setup',0,hook_handle)
154
155    if (present(iverbose)) then
156       iverb = iverbose
157    else
158       iverb = 2
159    end if
160
161    ! Open the aerosol scattering file and configure the way it is
162    ! read
163    call file%open(trim(file_name), iverbose=iverb)
164
165    if (file%exists('mass_ext_sw_hydrophilic')) then
166      this%use_hydrophilic = .true.
167    else
168      this%use_hydrophilic = .false.
169    end if
170
171    ! Read the raw scattering data
172    call file%get('mass_ext_sw_hydrophobic', this%mass_ext_sw_phobic)
173    call file%get('ssa_sw_hydrophobic',      this%ssa_sw_phobic)
174    call file%get('asymmetry_sw_hydrophobic',this%g_sw_phobic)
175    call file%get('mass_ext_lw_hydrophobic', this%mass_ext_lw_phobic)
176    call file%get('ssa_lw_hydrophobic',      this%ssa_lw_phobic)
177    call file%get('asymmetry_lw_hydrophobic',this%g_lw_phobic)
178
179    call file%get_global_attribute('description_hydrophobic', &
180         &                         this%description_phobic_str)
181
182    if (this%use_hydrophilic) then
183      call file%get('mass_ext_sw_hydrophilic', this%mass_ext_sw_philic)
184      call file%get('ssa_sw_hydrophilic',      this%ssa_sw_philic)
185      call file%get('asymmetry_sw_hydrophilic',this%g_sw_philic)
186      call file%get('mass_ext_lw_hydrophilic', this%mass_ext_lw_philic)
187      call file%get('ssa_lw_hydrophilic',      this%ssa_lw_philic)
188      call file%get('asymmetry_lw_hydrophilic',this%g_lw_philic)
189
190      call file%get('relative_humidity1',      this%rh_lower)
191
192      call file%get_global_attribute('description_hydrophilic', &
193           &                         this%description_philic_str)
194    end if
195
196    ! Read the raw scattering data at selected wavelengths if
197    ! available in the input file
198    if (file%exists('mass_ext_mono_hydrophobic')) then
199      this%use_monochromatic = .true.
200      call file%get('mass_ext_mono_hydrophobic', this%mass_ext_mono_phobic)
201      call file%get('ssa_mono_hydrophobic',      this%ssa_mono_phobic)
202      call file%get('asymmetry_mono_hydrophobic',this%g_mono_phobic)
203      call file%get('lidar_ratio_mono_hydrophobic',this%lidar_ratio_mono_phobic)
204      if (this%use_hydrophilic) then
205        call file%get('mass_ext_mono_hydrophilic', this%mass_ext_mono_philic)
206        call file%get('ssa_mono_hydrophilic',      this%ssa_mono_philic)
207        call file%get('asymmetry_mono_hydrophilic',this%g_mono_philic)
208        call file%get('lidar_ratio_mono_hydrophilic',this%lidar_ratio_mono_philic)
209      end if
210    else
211      this%use_monochromatic = .false.
212    end if
213
214    ! Close aerosol scattering file
215    call file%close()
216
217    ! Get array sizes
218    this%n_bands_lw    = size(this%mass_ext_lw_phobic, 1)
219    this%n_bands_sw    = size(this%mass_ext_sw_phobic, 1)
220    if (this%use_monochromatic) then
221      this%n_mono_wl   = size(this%mass_ext_mono_phobic, 1)
222    else
223      this%n_mono_wl   = 0
224    end if
225    this%n_type_phobic = size(this%mass_ext_lw_phobic, 2)
226
227    if (this%use_hydrophilic) then
228      this%n_type_philic = size(this%mass_ext_lw_philic, 3)
229      this%nrh           = size(this%mass_ext_lw_philic, 2)
230
231      ! Check agreement of dimensions
232      if (size(this%mass_ext_lw_philic,1) /= this%n_bands_lw) then
233        write(nulerr,'(a,a)') '*** Error: mass extinction for hydrophilic and hydrophobic ', &
234             &                'aerosol have different numbers of longwave bands'
235        call radiation_abort()
236      end if
237      if (size(this%mass_ext_sw_philic,1) /= this%n_bands_sw) then
238        write(nulerr,'(a,a)') '*** Error: mass extinction for hydrophilic and hydrophobic ', &
239             &                'aerosol have different numbers of shortwave bands'
240        call radiation_abort()
241      end if
242      if (size(this%rh_lower) /= this%nrh) then
243        write(nulerr,'(a)') '*** Error: size(relative_humidity1) /= size(mass_ext_sw_hydrophilic,2)'
244        call radiation_abort()
245      end if
246
247    else
248      this%n_type_philic = 0
249      this%nrh           = 0
250    end if
251
252    ! Allocate memory for mapping arrays
253    this%ntype = ntype
254    allocate(this%iclass(ntype))
255    allocate(this%itype(ntype))
256
257    this%iclass = IAerosolClassUndefined
258    this%itype  = 0
259
260    if (lhook) call dr_hook('radiation_aerosol_optics_data:setup',1,hook_handle)
261
262  end subroutine setup_aerosol_optics
263
264
265  !---------------------------------------------------------------------
266  ! Map user type "itype" onto stored hydrophobic type "i_type_phobic"
267  subroutine set_hydrophobic_type(this, itype, i_type_phobic)
268
269    use yomhook,     only : lhook, dr_hook
270
271    class(aerosol_optics_type), intent(inout) :: this
272    integer, intent(in)                       :: itype, i_type_phobic
273    real(jprb)                                :: hook_handle
274
275    if (lhook) call dr_hook('radiation_aerosol_optics_data:set_hydrophobic_type',0,hook_handle)
276
277    if (itype < 1 .or. itype > this%ntype) then
278      write(nulerr,'(a,i0)') '*** Error: aerosol type must be in the range 1 to ', &
279           &                  this%ntype
280       call radiation_abort('Error setting up aerosols')
281    end if
282    if (i_type_phobic < 1 .or. i_type_phobic > this%n_type_phobic) then
283      write(nulerr,'(a,i0)') '*** Error: hydrophobic type must be in the range 1 to ', &
284           &                  this%n_type_phobic
285      call radiation_abort('Error setting up aerosols')
286    end if
287
288    this%iclass(itype) = IAerosolClassHydrophobic
289    this%itype (itype) = i_type_phobic
290
291    if (lhook) call dr_hook('radiation_aerosol_optics_data:set_hydrophobic_type',1,hook_handle)
292
293  end subroutine set_hydrophobic_type
294
295
296  !---------------------------------------------------------------------
297  ! Map user type "itype" onto stored hydrophilic type "i_type_philic"
298  subroutine set_hydrophilic_type(this, itype, i_type_philic)
299
300    use yomhook,     only : lhook, dr_hook
301
302    class(aerosol_optics_type), intent(inout) :: this
303    integer, intent(in)                       :: itype, i_type_philic
304    real(jprb)                                :: hook_handle
305
306    if (lhook) call dr_hook('radiation_aerosol_optics_data:set_hydrophilic_type',0,hook_handle)
307
308    if (.not. this%use_hydrophilic) then
309      write(nulerr,'(a)') '*** Error: attempt to set hydrophilic aerosol type when no such types present'
310      call radiation_abort('Error setting up aerosols')     
311    end if
312
313    if (itype < 1 .or. itype > this%ntype) then
314      write(nulerr,'(a,i0)') '*** Error: aerosol type must be in the range 1 to ', &
315            &          this%ntype
316      call radiation_abort('Error setting up aerosols')
317    end if
318    if (i_type_philic < 1 .or. i_type_philic > this%n_type_philic) then
319      write(nulerr,'(a,i0)') '*** Error: hydrophilic type must be in the range 1 to ', &
320           &                 this%n_type_philic
321      call radiation_abort('Error setting up aerosols')
322    end if
323
324    this%iclass(itype) = IAerosolClassHydrophilic
325    this%itype (itype) = i_type_philic
326
327    if (lhook) call dr_hook('radiation_aerosol_optics_data:set_hydrophilic_type',1,hook_handle)
328
329  end subroutine set_hydrophilic_type
330
331
332  !---------------------------------------------------------------------
333  ! Set a user type "itype" to be ignored in the radiation scheme
334  subroutine set_empty_type(this, itype)
335
336    use yomhook,     only : lhook, dr_hook
337
338    class(aerosol_optics_type), intent(inout) :: this
339    integer, intent(in)                       :: itype
340    real(jprb)                                :: hook_handle
341
342    if (lhook) call dr_hook('radiation_aerosol_optics_data:set_empty_type',0,hook_handle)
343
344    if (itype < 1 .or. itype > this%ntype) then
345      write(nulerr,'(a,i0)') '*** Error: aerosol type must be in the range 1 to ', &
346           &                 this%ntype
347      call radiation_abort('Error setting up aerosols')
348    end if
349
350    this%iclass(itype) = IAerosolClassIgnored
351
352    if (lhook) call dr_hook('radiation_aerosol_optics_data:set_empty_type',1,hook_handle)
353
354  end subroutine set_empty_type
355
356
357  !---------------------------------------------------------------------
358  ! Set user types "itypes" to map onto the stored hydrophobic and
359  ! hydrophilic types according to its sign and value, with a value of
360  ! 0 indicating that this type is to be ignored.  Thus if itypes=(/
361  ! 3, 4, -6, 0 /) then user types 1 and 2 map on to hydrophobic types
362  ! 3 and 4, user type 3 maps on to hydrophilic type 6 and user type 4
363  ! is ignored.
364  subroutine set_types(this, itypes)
365
366    use yomhook,     only : lhook, dr_hook
367
368    class(aerosol_optics_type), intent(inout) :: this
369    integer, dimension(:), intent(in)         :: itypes
370
371    integer :: jtype
372    integer :: istart, iend
373    real(jprb)                                :: hook_handle
374
375    if (lhook) call dr_hook('radiation_aerosol_optics_data:set_types',0,hook_handle)
376
377    istart = lbound(itypes,1)
378    iend   = ubound(itypes,1)
379
380    do jtype = istart, iend
381      if (itypes(jtype) == 0) then
382        call this%set_empty_type(jtype)
383      else if (itypes(jtype) > 0) then
384        call this%set_hydrophobic_type(jtype, itypes(jtype))
385      else
386        call this%set_hydrophilic_type(jtype, -itypes(jtype))
387      end if
388    end do
389
390    if (lhook) call dr_hook('radiation_aerosol_optics_data:set_types',1,hook_handle)
391
392  end subroutine set_types
393
394
395  !---------------------------------------------------------------------
396  ! Return an index to the relative-humdity array, or zero if no
397  ! hydrophilic types are present. This function does so little that
398  ! it is best to remove the Dr Hook call.
399  function calc_rh_index(this, rh)
400
401    !use yomhook,     only : lhook, dr_hook
402
403    class(aerosol_optics_type), intent(inout) :: this
404    real(jprb),                 intent(in)    :: rh
405    integer                                   :: calc_rh_index
406    !real(jprb)                                :: hook_handle
407
408    !if (lhook) call dr_hook('radiation_aerosol_optics_data:calc_rh_index',0,hook_handle)
409
410    if (.not. this%use_hydrophilic) then
411      calc_rh_index = 0
412    else if (rh > this%rh_lower(this%nrh)) then
413      calc_rh_index = this%nrh
414    else
415      calc_rh_index = 1
416      do while (rh > this%rh_lower(calc_rh_index + 1))
417        calc_rh_index = calc_rh_index + 1
418      end do
419    end if
420
421    !if (lhook) call dr_hook('radiation_aerosol_optics_data:calc_rh_index',1,hook_handle)
422
423  end function calc_rh_index
424
425
426  !---------------------------------------------------------------------
427  ! Print a description of the aerosol types to nulout
428  subroutine print_description(this, i_type_map)
429
430    use radiation_io, only : nulout
431
432    class(aerosol_optics_type), intent(in) :: this
433    integer,                    intent(in) :: i_type_map(:)
434
435    integer :: jtype
436
437    if (size(i_type_map) > 0) then
438      write(nulout,'(a)') 'Aerosol mapping:'
439    else
440      write(nulout,'(a)') 'No aerosol types in radiation scheme'
441    end if
442
443    do jtype = 1,size(i_type_map)
444      if (i_type_map(jtype) > 0) then
445        write(nulout,'(i4,a,a)') jtype, ' -> hydrophobic type ', &
446             &  trim(get_line(this%description_phobic_str, i_type_map(jtype)))
447      else if (i_type_map(jtype) < 0) then
448        write(nulout,'(i4,a,a)') jtype, ' -> hydrophilic type ', &
449             &  trim(get_line(this%description_philic_str, -i_type_map(jtype)))
450      else
451        write(nulout,'(i4,a)') jtype, ' is unused'
452      end if
453    end do
454   
455  end subroutine print_description
456
457
458  !---------------------------------------------------------------------
459  ! Private helper function for print_description
460  pure function get_line(str,iline) result(line_str)
461    character(len=*), intent(in)  :: str
462    integer,          intent(in)  :: iline
463    character(len=NMaxLineLength) :: line_str
464   
465    integer :: istart, iend, i_start_new, ioffset, ilength, i_line_current
466    logical :: is_fail
467   
468    i_line_current = 1
469    istart = 1
470    iend = len(str)
471    is_fail = .false.
472    line_str = ' '
473
474    ! Find index of first character
475    do while (i_line_current < iline)
476      i_start_new = scan(str(istart:iend), new_line(' '))
477      if (i_start_new == 0) then
478        is_fail = .true.
479        cycle
480      else
481        istart = istart + i_start_new
482      end if
483      i_line_current = i_line_current + 1
484    end do
485   
486    if (.not. is_fail) then
487      ! Find index of last character
488      ioffset = scan(str(istart:iend), new_line(' '))
489      if (ioffset == 0) then
490        ilength = len(trim(str(istart:iend)))
491      else
492        ilength = ioffset - 1
493      end if
494     
495      if (ilength > NMaxLineLength) then
496        ilength = NMaxLineLength
497      end if
498      iend = istart + ilength - 1
499     
500      line_str = str(istart:iend)
501    else
502      write(line_str,'(i0,a)') iline, ': <unknown>'
503    end if
504   
505  end function get_line
506 
507end module radiation_aerosol_optics_data
Note: See TracBrowser for help on using the repository browser.