source: LMDZ6/branches/contrails/libf/phylmd/ecrad.v1.5.1/radiation_aerosol_optics_description.F90 @ 5452

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

Merge LMDZ_ECRad branch back into trunk!

File size: 12.2 KB
Line 
1! radiation_aerosol_optics_description.F90 - Type to store aerosol optics metadata
2!
3! (C) Copyright 2022- 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
16module radiation_aerosol_optics_description
17
18  use parkind1,      only : jprb
19
20  implicit none
21  public
22
23  !---------------------------------------------------------------------
24  ! This type holds the metadata from an aerosol optical property
25  ! file, enabling the user to request the index to the aerosol type
26  ! with particular properties.  Note that string information is held
27  ! in the form of 2D arrays of single characters, so comparison to
28  ! character strings requires the to_string helper function at the
29  ! end of this file.
30  type aerosol_optics_description_type
31
32    ! Two-character code describing the aerosol family, dimensioned
33    ! (2,naer), e.g.
34    !   SS: Sea salt
35    !   OM: Organic matter
36    !   SU: Sulfate
37    !   OB: Secondary organic biogenic
38    !   OA: Secondary organic anthropogenic
39    !   AM: Fine-mode ammonium sulfate
40    !   NI: Nitrate
41    !   DD: Desert dust
42    !   BC: Black carbon
43    character(len=1), allocatable :: code_phobic(:,:)
44    character(len=1), allocatable :: code_philic(:,:)
45
46    ! Size bin, typically 1-2 or 1-3 in order fine to coarse, or zero
47    ! if no division by size is used, dimensioned (naer)
48    integer, allocatable :: bin_phobic(:)
49    integer, allocatable :: bin_philic(:)
50
51    ! Character string characterizing the optical model, e.g. OPAC,
52    ! GACP, GLOMAP, Dubovik2002 etc.
53    character(len=1), allocatable :: optical_model_phobic(:,:)
54    character(len=1), allocatable :: optical_model_philic(:,:)
55
56    ! The user can call preferred_optical_model to specify that a
57    ! certain optical model for a certain aerosol family is to be
58    ! preferred when get_index is called
59    logical, allocatable :: is_preferred_phobic(:)
60    logical, allocatable :: is_preferred_philic(:)
61
62  contains
63    procedure :: read
64    procedure :: preferred_optical_model
65    procedure :: get_index
66
67  end type aerosol_optics_description_type
68
69contains
70
71  !---------------------------------------------------------------------
72  ! Read optical property file file_name into an
73  ! aerosol_optics_description_type object
74  subroutine read(this, file_name, iverbose)
75
76    use yomhook,              only : lhook, dr_hook
77    use easy_netcdf,          only : netcdf_file
78
79    class(aerosol_optics_description_type), intent(inout) :: this
80    character(len=*), intent(in)              :: file_name
81    integer, intent(in), optional             :: iverbose
82   
83    ! The NetCDF file containing the aerosol optics data
84    type(netcdf_file)  :: file
85
86    real(jprb)         :: hook_handle
87
88    if (lhook) call dr_hook('radiation_aerosol_optics_description:load',0,hook_handle)
89
90    ! Open the aerosol scattering file and configure the way it is
91    ! read
92    call file%open(trim(file_name), iverbose=iverbose)
93
94    ! Read metadata variables
95    call file%get('code_hydrophilic', this%code_philic)
96    call file%get('code_hydrophobic', this%code_phobic)
97    call file%get('bin_hydrophilic',  this%bin_philic)
98    call file%get('bin_hydrophobic',  this%bin_phobic)
99    call file%get('optical_model_hydrophilic', this%optical_model_philic)
100    call file%get('optical_model_hydrophobic', this%optical_model_phobic)
101
102    ! Allocate logical arrays of the appropriate size and set to FALSE
103    allocate(this%is_preferred_philic(size(this%bin_philic)))
104    allocate(this%is_preferred_phobic(size(this%bin_phobic)))
105    this%is_preferred_philic = .false.
106    this%is_preferred_phobic = .false.
107
108    call file%close()
109
110    if (lhook) call dr_hook('radiation_aerosol_optics_description:load',1,hook_handle)
111
112  end subroutine read
113
114  !---------------------------------------------------------------------
115  ! Specify the preferred optical model for a particular aerosol
116  ! family, e.g. "call
117  ! aer_desc%preferred_optical_model('DD','Woodward2001')" would mean
118  ! that subsequent calls to get_index in which the optical model is
119  ! not specified would return the Woodward model rather than the
120  ! first matching model in the file.  The check is only done on the
121  ! first len(optical_model_str) characters, so "Woodward" and
122  ! "Woodward2001" would both match the Woodward2001 model.
123  subroutine preferred_optical_model(this, code_str, optical_model_str)
124
125    use yomhook,              only : lhook, dr_hook
126
127    class(aerosol_optics_description_type), intent(inout) :: this
128    character(len=2), intent(in) :: code_str
129    character(len=*), intent(in) :: optical_model_str
130
131    ! Aerosol loop counter
132    integer :: ja
133
134    real(jprb)         :: hook_handle
135
136    if (lhook) call dr_hook('radiation_aerosol_optics_description:preferred_optical_model',0,hook_handle)
137
138    ! Loop over hydrophilic types
139    do ja = 1,size(this%bin_philic)
140      ! Check if we have a match
141      if (to_string(this%code_philic(:,ja)) == code_str &
142           &  .and. to_string(this%optical_model_philic(1:len(optical_model_str),ja)) &
143           &          == optical_model_str) then
144        this%is_preferred_philic(ja) = .true.
145      end if
146    end do
147    ! Repeat for the hydrophobic types
148    do ja = 1,size(this%bin_phobic)
149      if (to_string(this%code_phobic(:,ja)) == code_str &
150           &  .and. to_string(this%optical_model_phobic(1:len(optical_model_str),ja)) &
151           &          == optical_model_str) then
152        this%is_preferred_phobic(ja) = .true.
153      end if
154    end do
155
156    if (lhook) call dr_hook('radiation_aerosol_optics_description:preferred_optical_model',1,hook_handle)
157
158  end subroutine preferred_optical_model
159
160  !---------------------------------------------------------------------
161  ! Return the index to the aerosol optical properties corresponding
162  ! to the aerosol family in code_str (e.g. SS, DD etc), whether or
163  ! not the requested aerosol is hydrophilic in the logical
164  ! lhydrophilic, and optionally the size bin ibin and optical model
165  ! in optical_model_str. The return value may be used to populate the
166  ! radiation_config%i_aerosol_map vector, where a positive number is
167  ! a hydrophobic index, a negative number is a hydrophilic index and
168  ! zero indicates that the aerosol type was not found in the file.
169  ! This is a valid entry in i_aerosol_map meaning the aerosol is
170  ! ignored, but the calling routine to get_index might wish to throw
171  ! a warning or error. This routine works by assigning a score based
172  ! on the closeness of the match.
173  function get_index(this, code_str, lhydrophilic, ibin, optical_model_str)
174   
175    use yomhook,              only : lhook, dr_hook
176    use easy_netcdf,          only : netcdf_file
177    use radiation_io,         only : nulout
178
179    class(aerosol_optics_description_type), intent(in) :: this
180    character(len=2), intent(in) :: code_str
181    logical, intent(in) :: lhydrophilic
182    integer, intent(in), optional :: ibin
183    character(len=*), intent(in), optional :: optical_model_str
184
185    ! Score of the currently selected aerosol index, and the score of
186    ! the current one under consideration
187    integer :: score, current_score
188
189    ! Loop index for aerosol type
190    integer :: ja
191
192    ! Return value
193    integer :: get_index
194
195    ! Issue a warning if there is more than one equal match
196    logical :: is_ambiguous
197
198    real(jprb) :: hook_handle
199
200    if (lhook) call dr_hook('radiation_aerosol_optics_description:get_index',0,hook_handle)
201
202    ! Initial values
203    get_index = 0
204    score = 0
205    is_ambiguous = .false.
206
207    if (lhydrophilic) then
208      ! Loop over hydrophilic aerosol types
209      do ja = 1,size(this%bin_philic)
210        current_score = 0
211        if (to_string(this%code_philic(:,ja)) == code_str) then
212          ! Aerosol code matches
213          if (present(ibin) .and. this%bin_philic(ja) > 0) then
214            if (ibin > 0) then
215              if (ibin == this%bin_philic(ja)) then
216                ! Requested bin number matches
217                current_score = 4
218              else
219                ! Requested bin number does not match
220                current_score = -1
221              end if
222            else
223              ! Bin number is zero: no request
224              current_score = 2
225            end if
226          else
227            ! No bin number present
228            current_score = 2
229          end if
230          if (present(optical_model_str)) then
231            if (to_string(this%optical_model_philic(1:len(optical_model_str),ja)) &
232                 &  == optical_model_str) then
233              ! Requested optical model matches
234              if (current_score >= 0) then
235                current_score = current_score + 4
236              end if
237            else
238              ! Requested optical model does not match
239              current_score = -1
240            end if
241          else if (current_score >= 0) then
242            ! No requested optical model
243            current_score = current_score + 2
244          end if
245          if (current_score > 0 .and. this%is_preferred_philic(ja)) then
246            current_score = current_score + 1
247          end if
248          if (current_score > score) then
249            ! Better score than any existing aerosol type
250            get_index = -ja
251            score = current_score
252            is_ambiguous = .false.
253          else if (current_score > 0 .and. current_score == score) then
254            is_ambiguous = .true.
255          end if
256        end if
257      end do
258    else
259      ! Loop over hydrophobic aerosol types
260      do ja = 1,size(this%bin_phobic)
261        current_score = 0
262        if (to_string(this%code_phobic(:,ja)) == code_str) then
263          ! Aerosol code matches
264          if (present(ibin) .and. this%bin_phobic(ja) > 0) then
265            if (ibin > 0) then
266              if (ibin == this%bin_phobic(ja)) then
267                ! Requested bin number matches
268                current_score = 4
269              else
270                ! Requested bin number does not match
271                current_score = -1
272              end if
273            else
274              ! Bin number is zero: no request
275              current_score = 2
276            end if
277          else
278            ! No bin number requested or present
279            current_score = 2
280          end if
281          if (present(optical_model_str)) then
282            if (to_string(this%optical_model_phobic(1:len(optical_model_str),ja)) &
283                 &  == optical_model_str) then
284              ! Requested optical model matches
285              if (current_score >= 0) then
286                current_score = current_score + 4
287              end if
288            else
289              ! Requested optical model does not match
290              current_score = -1
291            end if
292          else if (current_score >= 0) then
293            ! No requested optical model
294            current_score = current_score + 2
295          end if
296          if (current_score > 0 .and. this%is_preferred_phobic(ja)) then
297            current_score = current_score + 1
298          end if
299          if (current_score > score) then
300            ! Better score than any existing aerosol type
301            get_index = ja
302            score = current_score
303            is_ambiguous = .false.
304          else if (current_score > 0 .and. current_score == score) then
305            is_ambiguous = .true.
306          end if         
307        end if
308      end do
309    end if
310
311    if (is_ambiguous) then
312      write(nulout,'(a,a2,a,l1,a)') 'Warning: get_index("', code_str, '",', lhydrophilic, &
313           &  ',...) does not unambiguously identify an aerosol optical property index'
314    end if
315
316    if (lhook) call dr_hook('radiation_aerosol_optics_description:get_index',1,hook_handle)
317
318  end function get_index
319
320  !---------------------------------------------------------------------
321  ! Utility function to convert an array of single characters to a
322  ! character string (yes Fortran's string handling is a bit rubbish)
323  pure function to_string(arr) result(str)
324    character, intent(in)  :: arr(:)
325    character(len=size(arr)) :: str
326    integer :: jc
327    do jc = 1,size(arr)
328      str(jc:jc) = arr(jc)
329    end do
330  end function to_string
331
332end module radiation_aerosol_optics_description
Note: See TracBrowser for help on using the repository browser.