source: LMDZ6/branches/contrails/libf/phylmd/ecrad.v1.5.1/radiation_aerosol_optics_data.F90 @ 5431

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

Merge LMDZ_ECRad branch back into trunk!

File size: 29.4 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     ! Wavenumber (cm-1) upper and lower bounds of each spectral
60     ! interval, which if used in the RRTMG gas optics scheme should
61     ! match its band bounds
62     real(jprb), allocatable, dimension(:) :: wavenumber1_sw, wavenumber2_sw
63     real(jprb), allocatable, dimension(:) :: wavenumber1_lw, wavenumber2_lw
64
65     ! Scattering properties are provided separately in the shortwave
66     ! and longwave for hydrophobic and hydrophilic aerosols.
67     ! Hydrophobic aerosols are dimensioned (nband,n_type_phobic):
68     real(jprb), allocatable, dimension(:,:) :: &
69          &  mass_ext_sw_phobic, & ! Mass-extinction coefficient (m2 kg-1)
70          &  ssa_sw_phobic,      & ! Single scattering albedo
71          &  g_sw_phobic,        & ! Asymmetry factor
72!          &  ssa_g_sw_phobic,    & ! ssa*g
73          &  mass_ext_lw_phobic, & ! Mass-extinction coefficient (m2 kg-1)
74!          &  mass_abs_lw_phobic, & ! Mass-absorption coefficient (m2 kg-1)
75          &  ssa_lw_phobic,      & ! Single scattering albedo
76          &  g_lw_phobic           ! Asymmetry factor
77
78     ! Hydrophilic aerosols are dimensioned (nband,nrh,n_type_philic):
79     real(jprb), allocatable, dimension(:,:,:) :: &
80          &  mass_ext_sw_philic, & ! Mass-extinction coefficient (m2 kg-1)
81          &  ssa_sw_philic,      & ! Single scattering albedo
82          &  g_sw_philic,        & ! Asymmetry factor
83 !         &  ssa_g_sw_philic,    & ! ssa*g
84          &  mass_ext_lw_philic, & ! Mass-extinction coefficient (m2 kg-1)
85 !         &  mass_abs_lw_philic, & ! Mass-absorption coefficient (m2 kg-1)
86          &  ssa_lw_philic,      & ! Single scattering albedo
87          &  g_lw_philic           ! Asymmetry factor
88
89     ! Wavelengths at which monochromatic properties are stored,
90     ! dimensioned (n_mono_wl), units metres
91     real(jprb), allocatable :: wavelength_mono(:)
92
93     ! Scattering properties at selected monochromatic wavelengths
94     ! (n_mono_wl,n_type_phobic)
95     real(jprb), allocatable, dimension(:,:) :: &
96          &  mass_ext_mono_phobic, & ! Mass-extinction coefficient (m2 kg-1)
97          &  ssa_mono_phobic,      & ! Single scattering albedo
98          &  g_mono_phobic,        & ! Asymmetry factor
99          &  lidar_ratio_mono_phobic ! Lidar Ratio
100     ! ...hydrophilic aerosols dimensioned (n_mono_wl,nrh,n_type_philic):
101     real(jprb), allocatable, dimension(:,:,:) :: &
102          &  mass_ext_mono_philic, & ! Mass-extinction coefficient (m2 kg-1)
103          &  ssa_mono_philic,      & ! Single scattering albedo
104          &  g_mono_philic,        & ! Asymmetry factor
105          &  lidar_ratio_mono_philic ! Lidar Ratio
106
107     ! For hydrophilic aerosols, the lower bounds of the relative
108     ! humidity bins is a vector of length nrh:
109     real(jprb), allocatable, dimension(:) :: &
110          &  rh_lower    ! Dimensionless (1.0 = 100% humidity)
111
112     ! Strings describing the aerosol types
113     character(len=NMaxStringLength) :: description_phobic_str = ' '
114     character(len=NMaxStringLength) :: description_philic_str = ' '
115
116     ! The number of user-defined aerosol types
117     integer :: ntype
118
119     ! The number of hydrophobic and hydrophilic types read from the
120     ! aerosol optics file
121     integer :: n_type_phobic = 0
122     integer :: n_type_philic = 0
123
124     ! Number of relative humidity bins
125     integer :: nrh = 0
126
127     ! Number of longwave and shortwave bands of the data in the file,
128     ! and monochromatic wavelengths
129     integer :: n_bands_lw = 0, n_bands_sw = 0, n_mono_wl = 0
130
131     ! Do we have any hydrophilic types?
132     logical :: use_hydrophilic = .true.
133
134     ! Do we have monochromatic optical properties
135     logical :: use_monochromatic = .false.
136
137   contains
138     procedure :: setup => setup_aerosol_optics
139     procedure :: save  => save_aerosol_optics
140     procedure :: allocate
141     procedure :: initialize_types
142     procedure :: set_hydrophobic_type
143     procedure :: set_hydrophilic_type
144     procedure :: set_empty_type
145     procedure :: set_types
146     procedure :: calc_rh_index
147     procedure :: print_description
148
149  end type aerosol_optics_type
150
151contains
152
153
154  !---------------------------------------------------------------------
155  ! Setup aerosol optics coefficients by reading them from a file
156  subroutine setup_aerosol_optics(this, file_name, iverbose)
157
158    use yomhook,              only : lhook, dr_hook
159    use easy_netcdf,          only : netcdf_file
160    use radiation_io,         only : nulerr, radiation_abort
161
162    class(aerosol_optics_type), intent(inout) :: this
163    character(len=*), intent(in)              :: file_name
164    integer, intent(in), optional             :: iverbose
165
166    ! The NetCDF file containing the aerosol optics data
167    type(netcdf_file)  :: file
168
169    real(jprb), allocatable :: wavelength_tmp(:)
170    integer            :: iverb
171
172    real(jprb)         :: hook_handle
173
174    if (lhook) call dr_hook('radiation_aerosol_optics_data:setup',0,hook_handle)
175
176    if (present(iverbose)) then
177       iverb = iverbose
178    else
179       iverb = 2
180    end if
181
182    ! Open the aerosol scattering file and configure the way it is
183    ! read
184    call file%open(trim(file_name), iverbose=iverb)
185
186    if (file%exists('mass_ext_sw_hydrophilic')) then
187      this%use_hydrophilic = .true.
188    else
189      this%use_hydrophilic = .false.
190    end if
191
192    ! Read the wavenumber bounds
193    call file%get('wavenumber1_sw', this%wavenumber1_sw)
194    call file%get('wavenumber2_sw', this%wavenumber2_sw)
195    call file%get('wavenumber1_lw', this%wavenumber1_lw)
196    call file%get('wavenumber2_lw', this%wavenumber2_lw)
197
198    ! Read the raw scattering data
199    call file%get('mass_ext_sw_hydrophobic', this%mass_ext_sw_phobic)
200    call file%get('ssa_sw_hydrophobic',      this%ssa_sw_phobic)
201    call file%get('asymmetry_sw_hydrophobic',this%g_sw_phobic)
202    call file%get('mass_ext_lw_hydrophobic', this%mass_ext_lw_phobic)
203    call file%get('ssa_lw_hydrophobic',      this%ssa_lw_phobic)
204    call file%get('asymmetry_lw_hydrophobic',this%g_lw_phobic)
205
206    call file%get_global_attribute('description_hydrophobic', &
207         &                         this%description_phobic_str)
208
209    ! Precompute ssa*g for the shortwave and mass-absorption for the
210    ! longwave - TBD FIX
211    !allocate(this%ssa_g_sw_phobic(...
212
213    if (this%use_hydrophilic) then
214      call file%get('mass_ext_sw_hydrophilic', this%mass_ext_sw_philic)
215      call file%get('ssa_sw_hydrophilic',      this%ssa_sw_philic)
216      call file%get('asymmetry_sw_hydrophilic',this%g_sw_philic)
217      call file%get('mass_ext_lw_hydrophilic', this%mass_ext_lw_philic)
218      call file%get('ssa_lw_hydrophilic',      this%ssa_lw_philic)
219      call file%get('asymmetry_lw_hydrophilic',this%g_lw_philic)
220
221      call file%get('relative_humidity1',      this%rh_lower)
222
223      call file%get_global_attribute('description_hydrophilic', &
224           &                         this%description_philic_str)
225    end if
226
227    ! Read the monochromatic scattering data at selected wavelengths
228    ! if available in the input file
229    if (file%exists('mass_ext_mono_hydrophobic')) then
230      this%use_monochromatic = .true.
231
232      if (allocated(this%wavelength_mono)) then
233        ! User has provided required monochromatic wavelengths, which
234        ! must match those in the file (in the more recent "general"
235        ! aerosol optics, interpolation provides optical properties at
236        ! the requested wavelengths)
237        call file%get('wavelength_mono', wavelength_tmp)
238        if (size(wavelength_tmp) /= size(this%wavelength_mono)) then
239          write(nulerr,'(a,i0,a,i0,a)') '*** Error: ', size(this%wavelength_mono), &
240               &  ' monochromatic wavelengths requested but ', &
241               &  size(wavelength_tmp), ' in file'
242          call radiation_abort('Radiation configuration error')
243        end if
244        if (any(abs(this%wavelength_mono-wavelength_tmp) &
245               &  / this%wavelength_mono > 0.01_jprb)) then
246          write(nulerr,'(a,a)') '*** Error: requested monochromatic wavelengths', &
247               &  'must all be within 1% of values in file'
248          call radiation_abort('Radiation configuration error')
249        end if
250      else
251        ! User has not provided required wavelengths, so we save the
252        ! monochromatic wavelengths in the file
253        call file%get('wavelength_mono', this%wavelength_mono)
254      end if
255
256      call file%get('mass_ext_mono_hydrophobic', this%mass_ext_mono_phobic)
257      call file%get('ssa_mono_hydrophobic',      this%ssa_mono_phobic)
258      call file%get('asymmetry_mono_hydrophobic',this%g_mono_phobic)
259      call file%get('lidar_ratio_mono_hydrophobic',this%lidar_ratio_mono_phobic)
260      if (this%use_hydrophilic) then
261        call file%get('mass_ext_mono_hydrophilic', this%mass_ext_mono_philic)
262        call file%get('ssa_mono_hydrophilic',      this%ssa_mono_philic)
263        call file%get('asymmetry_mono_hydrophilic',this%g_mono_philic)
264        call file%get('lidar_ratio_mono_hydrophilic',this%lidar_ratio_mono_philic)
265      end if
266    else
267      this%use_monochromatic = .false.
268    end if
269
270    ! Close aerosol scattering file
271    call file%close()
272
273    ! Get array sizes
274    this%n_bands_lw    = size(this%mass_ext_lw_phobic, 1)
275    this%n_bands_sw    = size(this%mass_ext_sw_phobic, 1)
276    if (this%use_monochromatic) then
277      this%n_mono_wl   = size(this%mass_ext_mono_phobic, 1)
278    else
279      this%n_mono_wl   = 0
280    end if
281    this%n_type_phobic = size(this%mass_ext_lw_phobic, 2)
282
283    if (this%use_hydrophilic) then
284      this%n_type_philic = size(this%mass_ext_lw_philic, 3)
285      this%nrh           = size(this%mass_ext_lw_philic, 2)
286
287      ! Check agreement of dimensions
288      if (size(this%mass_ext_lw_philic,1) /= this%n_bands_lw) then
289        write(nulerr,'(a,a)') '*** Error: mass extinction for hydrophilic and hydrophobic ', &
290             &                'aerosol have different numbers of longwave bands'
291        call radiation_abort('Radiation configuration error')
292      end if
293      if (size(this%mass_ext_sw_philic,1) /= this%n_bands_sw) then
294        write(nulerr,'(a,a)') '*** Error: mass extinction for hydrophilic and hydrophobic ', &
295             &                'aerosol have different numbers of shortwave bands'
296        call radiation_abort('Radiation configuration error')
297      end if
298      if (size(this%rh_lower) /= this%nrh) then
299        write(nulerr,'(a)') '*** Error: size(relative_humidity1) /= size(mass_ext_sw_hydrophilic,2)'
300        call radiation_abort('Radiation configuration error')
301      end if
302
303    else
304      this%n_type_philic = 0
305      this%nrh           = 0
306    end if
307
308    if (lhook) call dr_hook('radiation_aerosol_optics_data:setup',1,hook_handle)
309
310  end subroutine setup_aerosol_optics
311
312
313  !---------------------------------------------------------------------
314  ! Initialize the arrays describing the user's aerosol types
315  subroutine initialize_types(this, ntype)
316
317    class(aerosol_optics_type), intent(inout) :: this
318    integer,                    intent(in)    :: ntype
319   
320    ! Allocate memory for mapping arrays
321    this%ntype = ntype
322    allocate(this%iclass(ntype))
323    allocate(this%itype(ntype))
324
325    this%iclass = IAerosolClassUndefined
326    this%itype  = 0
327
328  end subroutine initialize_types
329
330  !---------------------------------------------------------------------
331  ! Allocate arrays for aerosol optics data type
332  subroutine allocate(this, n_type_phobic, n_type_philic, nrh, &
333       &              n_bands_lw, n_bands_sw, n_mono_wl)
334
335    use yomhook,     only : lhook, dr_hook
336
337    class(aerosol_optics_type), intent(inout) :: this
338    integer, intent(in) :: n_type_phobic, n_type_philic, nrh
339    integer, intent(in) :: n_bands_lw, n_bands_sw, n_mono_wl
340
341    real(jprb) :: hook_handle
342
343    if (lhook) call dr_hook('radiation_aerosol_optics_data:allocate',0,hook_handle)
344
345    this%n_type_phobic = n_type_phobic
346    this%n_type_philic = n_type_philic
347    this%nrh           = nrh
348    this%n_bands_lw    = n_bands_lw
349    this%n_bands_sw    = n_bands_sw
350    this%n_mono_wl     = n_mono_wl
351
352    if (n_type_philic > 0) then
353      this%use_hydrophilic = .true.
354    else
355      this%use_hydrophilic = .false.
356    end if
357
358    if (n_bands_sw > 0) then
359      allocate(this%mass_ext_sw_phobic(n_bands_sw, n_type_phobic))
360      allocate(this%ssa_sw_phobic(n_bands_sw, n_type_phobic))
361      allocate(this%g_sw_phobic(n_bands_sw, n_type_phobic))
362    end if
363    if (n_bands_lw > 0) then
364      allocate(this%mass_ext_lw_phobic(n_bands_lw, n_type_phobic))
365      allocate(this%ssa_lw_phobic(n_bands_lw, n_type_phobic))
366      allocate(this%g_lw_phobic(n_bands_lw, n_type_phobic))
367    end if
368    if (n_mono_wl > 0) then
369      allocate(this%mass_ext_mono_phobic(n_mono_wl, n_type_phobic))
370      allocate(this%ssa_mono_phobic(n_mono_wl, n_type_phobic))
371      allocate(this%g_mono_phobic(n_mono_wl, n_type_phobic))
372      allocate(this%lidar_ratio_mono_phobic(n_mono_wl, n_type_phobic))
373    end if
374
375    if (n_type_philic > 0 .and. nrh > 0) then
376      if (n_bands_sw > 0) then
377        allocate(this%mass_ext_sw_philic(n_bands_sw, nrh, n_type_philic))
378        allocate(this%ssa_sw_philic(n_bands_sw, nrh, n_type_philic))
379        allocate(this%g_sw_philic(n_bands_sw, nrh, n_type_philic))
380      end if
381      if (n_bands_lw > 0) then
382        allocate(this%mass_ext_lw_philic(n_bands_lw, nrh, n_type_philic))
383        allocate(this%ssa_lw_philic(n_bands_lw, nrh, n_type_philic))
384        allocate(this%g_lw_philic(n_bands_lw, nrh, n_type_philic))
385      end if
386      if (n_mono_wl > 0) then
387        allocate(this%mass_ext_mono_philic(n_mono_wl, nrh, n_type_philic))
388        allocate(this%ssa_mono_philic(n_mono_wl, nrh, n_type_philic))
389        allocate(this%g_mono_philic(n_mono_wl, nrh, n_type_philic))
390        allocate(this%lidar_ratio_mono_philic(n_mono_wl, nrh, n_type_philic))
391      end if
392    end if
393
394    if (lhook) call dr_hook('radiation_aerosol_optics_data:allocate',1,hook_handle)
395
396  end subroutine allocate
397
398
399  !---------------------------------------------------------------------
400  subroutine save_aerosol_optics(this, file_name, iverbose)
401
402    use yomhook,     only : lhook, dr_hook
403    use easy_netcdf, only : netcdf_file
404
405    class(aerosol_optics_type), intent(inout) :: this
406    character(len=*),           intent(in)    :: file_name
407    integer,          optional, intent(in)    :: iverbose
408
409    ! Object for output NetCDF file
410    type(netcdf_file) :: out_file
411
412    real(jprb) :: hook_handle
413
414    if (lhook) call dr_hook('radiation_aerosol_optics_data:save',0,hook_handle)
415
416    ! Create the file
417    call out_file%create(trim(file_name), iverbose=iverbose)
418
419    ! Define dimensions
420    call out_file%define_dimension("band_lw", this%n_bands_lw)
421    call out_file%define_dimension("band_sw", this%n_bands_sw)
422    call out_file%define_dimension("hydrophilic", this%n_type_philic)
423    call out_file%define_dimension("hydrophobic", this%n_type_phobic)
424    call out_file%define_dimension("relative_humidity", this%nrh)
425    !if (this%use_monochromatic) then
426    !  call out_file%define_dimension("wavelength_mono", this%n_mono_wl)
427    !end if
428
429    ! Put global attributes
430    call out_file%put_global_attributes( &
431         &   title_str="Aerosol optical properties in the spectral intervals of the gas-optics scheme for ecRad", &
432         &   source_str="ecRad offline radiation model")
433    call out_file%put_global_attribute( &
434         &  "description_hydrophobic", this%description_phobic_str)
435    call out_file%put_global_attribute( &
436         &  "description_hydrophilic", this%description_philic_str)
437
438    ! Define variables
439    call out_file%define_variable("mass_ext_sw_hydrophobic", units_str="m2 kg-1", &
440         &  long_name="Shortwave mass-extinction coefficient of hydrophobic aerosols", &
441         &  dim2_name="hydrophobic", dim1_name="band_sw")
442    call out_file%define_variable("ssa_sw_hydrophobic", units_str="1", &
443         &  long_name="Shortwave single scattering albedo of hydrophobic aerosols", &
444         &  dim2_name="hydrophobic", dim1_name="band_sw")
445    call out_file%define_variable("asymmetry_sw_hydrophobic", units_str="1", &
446         &  long_name="Shortwave asymmetry factor of hydrophobic aerosols", &
447         &  dim2_name="hydrophobic", dim1_name="band_sw")
448
449    call out_file%define_variable("mass_ext_lw_hydrophobic", units_str="m2 kg-1", &
450         &  long_name="Longwave mass-extinction coefficient of hydrophobic aerosols", &
451         &  dim2_name="hydrophobic", dim1_name="band_lw")
452    call out_file%define_variable("ssa_lw_hydrophobic", units_str="1", &
453         &  long_name="Longwave single scattering albedo of hydrophobic aerosols", &
454         &  dim2_name="hydrophobic", dim1_name="band_lw")
455    call out_file%define_variable("asymmetry_lw_hydrophobic", units_str="1", &
456         &  long_name="Longwave asymmetry factor of hydrophobic aerosols", &
457         &  dim2_name="hydrophobic", dim1_name="band_lw")
458
459    call out_file%define_variable("mass_ext_sw_hydrophilic", units_str="m2 kg-1", &
460         &  long_name="Shortwave mass-extinction coefficient of hydrophilic aerosols", &
461         &  dim3_name="hydrophilic", dim2_name="relative_humidity", dim1_name="band_sw")
462    call out_file%define_variable("ssa_sw_hydrophilic", units_str="1", &
463         &  long_name="Shortwave single scattering albedo of hydrophilic aerosols", &
464         &  dim3_name="hydrophilic", dim2_name="relative_humidity", dim1_name="band_sw")
465    call out_file%define_variable("asymmetry_sw_hydrophilic", units_str="1", &
466         &  long_name="Shortwave asymmetry factor of hydrophilic aerosols", &
467         &  dim3_name="hydrophilic", dim2_name="relative_humidity", dim1_name="band_sw")
468
469    call out_file%define_variable("mass_ext_lw_hydrophilic", units_str="m2 kg-1", &
470         &  long_name="Longwave mass-extinction coefficient of hydrophilic aerosols", &
471         &  dim3_name="hydrophilic", dim2_name="relative_humidity", dim1_name="band_lw")
472    call out_file%define_variable("ssa_lw_hydrophilic", units_str="1", &
473         &  long_name="Longwave single scattering albedo of hydrophilic aerosols", &
474         &  dim3_name="hydrophilic", dim2_name="relative_humidity", dim1_name="band_lw")
475    call out_file%define_variable("asymmetry_lw_hydrophilic", units_str="1", &
476         &  long_name="Longwave asymmetry factor of hydrophilic aerosols", &
477         &  dim3_name="hydrophilic", dim2_name="relative_humidity", dim1_name="band_lw")
478
479    ! Write variables
480    call out_file%put("mass_ext_sw_hydrophobic", this%mass_ext_sw_phobic)
481    call out_file%put("ssa_sw_hydrophobic", this%ssa_sw_phobic)
482    call out_file%put("asymmetry_sw_hydrophobic", this%g_sw_phobic)
483    call out_file%put("mass_ext_lw_hydrophobic", this%mass_ext_lw_phobic)
484    call out_file%put("ssa_lw_hydrophobic", this%ssa_lw_phobic)
485    call out_file%put("asymmetry_lw_hydrophobic", this%g_lw_phobic)
486    call out_file%put("mass_ext_sw_hydrophilic", this%mass_ext_sw_philic)
487    call out_file%put("ssa_sw_hydrophilic", this%ssa_sw_philic)
488    call out_file%put("asymmetry_sw_hydrophilic", this%g_sw_philic)
489    call out_file%put("mass_ext_lw_hydrophilic", this%mass_ext_lw_philic)
490    call out_file%put("ssa_lw_hydrophilic", this%ssa_lw_philic)
491    call out_file%put("asymmetry_lw_hydrophilic", this%g_lw_philic)
492
493    call out_file%close()
494
495    if (lhook) call dr_hook('radiation_aerosol_optics_data:save',1,hook_handle)
496
497  end subroutine save_aerosol_optics
498
499
500  !---------------------------------------------------------------------
501  ! Map user type "itype" onto stored hydrophobic type "i_type_phobic"
502  subroutine set_hydrophobic_type(this, itype, i_type_phobic)
503
504    use yomhook,     only : lhook, dr_hook
505
506    class(aerosol_optics_type), intent(inout) :: this
507    integer, intent(in)                       :: itype, i_type_phobic
508    real(jprb)                                :: hook_handle
509
510    if (lhook) call dr_hook('radiation_aerosol_optics_data:set_hydrophobic_type',0,hook_handle)
511
512    if (itype < 1 .or. itype > this%ntype) then
513      write(nulerr,'(a,i0)') '*** Error: aerosol type must be in the range 1 to ', &
514           &                  this%ntype
515       call radiation_abort('Error setting up aerosols')
516    end if
517    if (i_type_phobic < 1 .or. i_type_phobic > this%n_type_phobic) then
518      write(nulerr,'(a,i0)') '*** Error: hydrophobic type must be in the range 1 to ', &
519           &                  this%n_type_phobic
520      call radiation_abort('Error setting up aerosols')
521    end if
522
523    this%iclass(itype) = IAerosolClassHydrophobic
524    this%itype (itype) = i_type_phobic
525
526    if (lhook) call dr_hook('radiation_aerosol_optics_data:set_hydrophobic_type',1,hook_handle)
527
528  end subroutine set_hydrophobic_type
529
530
531  !---------------------------------------------------------------------
532  ! Map user type "itype" onto stored hydrophilic type "i_type_philic"
533  subroutine set_hydrophilic_type(this, itype, i_type_philic)
534
535    use yomhook,     only : lhook, dr_hook
536
537    class(aerosol_optics_type), intent(inout) :: this
538    integer, intent(in)                       :: itype, i_type_philic
539    real(jprb)                                :: hook_handle
540
541    if (lhook) call dr_hook('radiation_aerosol_optics_data:set_hydrophilic_type',0,hook_handle)
542
543    if (.not. this%use_hydrophilic) then
544      write(nulerr,'(a)') '*** Error: attempt to set hydrophilic aerosol type when no such types present'
545      call radiation_abort('Error setting up aerosols')     
546    end if
547
548    if (itype < 1 .or. itype > this%ntype) then
549      write(nulerr,'(a,i0)') '*** Error: aerosol type must be in the range 1 to ', &
550            &          this%ntype
551      call radiation_abort('Error setting up aerosols')
552    end if
553    if (i_type_philic < 1 .or. i_type_philic > this%n_type_philic) then
554      write(nulerr,'(a,i0)') '*** Error: hydrophilic type must be in the range 1 to ', &
555           &                 this%n_type_philic
556      call radiation_abort('Error setting up aerosols')
557    end if
558
559    this%iclass(itype) = IAerosolClassHydrophilic
560    this%itype (itype) = i_type_philic
561
562    if (lhook) call dr_hook('radiation_aerosol_optics_data:set_hydrophilic_type',1,hook_handle)
563
564  end subroutine set_hydrophilic_type
565
566
567  !---------------------------------------------------------------------
568  ! Set a user type "itype" to be ignored in the radiation scheme
569  subroutine set_empty_type(this, itype)
570
571    use yomhook,     only : lhook, dr_hook
572
573    class(aerosol_optics_type), intent(inout) :: this
574    integer, intent(in)                       :: itype
575    real(jprb)                                :: hook_handle
576
577    if (lhook) call dr_hook('radiation_aerosol_optics_data:set_empty_type',0,hook_handle)
578
579    if (itype < 1 .or. itype > this%ntype) then
580      write(nulerr,'(a,i0)') '*** Error: aerosol type must be in the range 1 to ', &
581           &                 this%ntype
582      call radiation_abort('Error setting up aerosols')
583    end if
584
585    this%iclass(itype) = IAerosolClassIgnored
586
587    if (lhook) call dr_hook('radiation_aerosol_optics_data:set_empty_type',1,hook_handle)
588
589  end subroutine set_empty_type
590
591
592  !---------------------------------------------------------------------
593  ! Set user types "itypes" to map onto the stored hydrophobic and
594  ! hydrophilic types according to its sign and value, with a value of
595  ! 0 indicating that this type is to be ignored.  Thus if itypes=(/
596  ! 3, 4, -6, 0 /) then user types 1 and 2 map on to hydrophobic types
597  ! 3 and 4, user type 3 maps on to hydrophilic type 6 and user type 4
598  ! is ignored.
599  subroutine set_types(this, itypes)
600
601    use yomhook,     only : lhook, dr_hook
602
603    class(aerosol_optics_type), intent(inout) :: this
604    integer, dimension(:), intent(in)         :: itypes
605
606    integer :: jtype
607    integer :: istart, iend
608    real(jprb)                                :: hook_handle
609
610    if (lhook) call dr_hook('radiation_aerosol_optics_data:set_types',0,hook_handle)
611
612    istart = lbound(itypes,1)
613    iend   = ubound(itypes,1)
614
615    do jtype = istart, iend
616      if (itypes(jtype) == 0) then
617        call this%set_empty_type(jtype)
618      else if (itypes(jtype) > 0) then
619        call this%set_hydrophobic_type(jtype, itypes(jtype))
620      else
621        call this%set_hydrophilic_type(jtype, -itypes(jtype))
622      end if
623    end do
624
625    if (lhook) call dr_hook('radiation_aerosol_optics_data:set_types',1,hook_handle)
626
627  end subroutine set_types
628
629
630  !---------------------------------------------------------------------
631  ! Return an index to the relative-humdity array, or zero if no
632  ! hydrophilic types are present. This function does so little that
633  ! it is best to remove the Dr Hook call.
634  function calc_rh_index(this, rh)
635
636    !use yomhook,     only : lhook, dr_hook
637
638    class(aerosol_optics_type), intent(inout) :: this
639    real(jprb),                 intent(in)    :: rh
640    integer                                   :: calc_rh_index
641    !real(jprb)                                :: hook_handle
642
643    !if (lhook) call dr_hook('radiation_aerosol_optics_data:calc_rh_index',0,hook_handle)
644
645    if (.not. this%use_hydrophilic) then
646      calc_rh_index = 0
647    else if (rh > this%rh_lower(this%nrh)) then
648      calc_rh_index = this%nrh
649    else
650      calc_rh_index = 1
651      do while (rh > this%rh_lower(calc_rh_index + 1))
652        calc_rh_index = calc_rh_index + 1
653      end do
654    end if
655
656    !if (lhook) call dr_hook('radiation_aerosol_optics_data:calc_rh_index',1,hook_handle)
657
658  end function calc_rh_index
659
660
661  !---------------------------------------------------------------------
662  ! Print a description of the aerosol types to nulout
663  subroutine print_description(this, i_type_map)
664
665    use radiation_io, only : nulout
666
667    class(aerosol_optics_type), intent(in) :: this
668    integer,                    intent(in) :: i_type_map(:)
669
670    integer :: jtype
671
672    if (size(i_type_map) > 0) then
673      write(nulout,'(a)') 'Aerosol mapping:'
674    else
675      write(nulout,'(a)') 'No aerosol types in radiation scheme'
676    end if
677
678    do jtype = 1,size(i_type_map)
679      if (i_type_map(jtype) > 0) then
680        write(nulout,'(i4,a,a)') jtype, ' -> hydrophobic type ', &
681             &  trim(get_line(this%description_phobic_str, i_type_map(jtype)))
682      else if (i_type_map(jtype) < 0) then
683        write(nulout,'(i4,a,a)') jtype, ' -> hydrophilic type ', &
684             &  trim(get_line(this%description_philic_str, -i_type_map(jtype)))
685      else
686        write(nulout,'(i4,a)') jtype, ' is unused'
687      end if
688    end do
689   
690  end subroutine print_description
691
692
693  !---------------------------------------------------------------------
694  ! Private helper function for print_description
695  pure function get_line(str,iline) result(line_str)
696    character(len=*), intent(in)  :: str
697    integer,          intent(in)  :: iline
698    character(len=NMaxLineLength) :: line_str
699   
700    integer :: istart, iend, i_start_new, ioffset, ilength, i_line_current
701    logical :: is_fail
702   
703    i_line_current = 1
704    istart = 1
705    iend = len(str)
706    is_fail = .false.
707    line_str = ' '
708
709    ! Find index of first character
710    do while (i_line_current < iline)
711      i_start_new = scan(str(istart:iend), new_line(' '))
712      if (i_start_new == 0) then
713        is_fail = .true.
714        cycle
715      else
716        istart = istart + i_start_new
717      end if
718      i_line_current = i_line_current + 1
719    end do
720   
721    if (.not. is_fail) then
722      ! Find index of last character
723      ioffset = scan(str(istart:iend), new_line(' '))
724      if (ioffset == 0) then
725        ilength = len(trim(str(istart:iend)))
726      else
727        ilength = ioffset - 1
728      end if
729     
730      if (ilength > NMaxLineLength) then
731        ilength = NMaxLineLength
732      end if
733      iend = istart + ilength - 1
734     
735      line_str = str(istart:iend)
736    else
737      write(line_str,'(i0,a)') iline, ': <unknown>'
738    end if
739   
740  end function get_line
741 
742end module radiation_aerosol_optics_data
Note: See TracBrowser for help on using the repository browser.