Line data Source code
1 : ! This code is part of Radiative Transfer for Energetics (RTE)
2 : !
3 : ! Contacts: Robert Pincus and Eli Mlawer
4 : ! email: rrtmgp@aer.com
5 : !
6 : ! Copyright 2015- Atmospheric and Environmental Research,
7 : ! Regents of the University of Colorado,
8 : ! Trustees of Columbia University in the City of New York
9 : ! All right reserved.
10 : !
11 : ! Use and duplication is permitted under the terms of the
12 : ! BSD 3-clause license, see http://opensource.org/licenses/BSD-3-Clause
13 : ! -------------------------------------------------------------------------------------------------
14 : !>
15 : !> ## Encapsulate optical properties defined on a spectral grid of N bands.
16 : !>
17 : !> The bands are described by their limiting wavenumbers. They need not be contiguous or complete.
18 : !> A band may contain more than one spectral sub-point (g-point) in which case a mapping must be supplied.
19 : !> A name may be provided and will be prepended to error messages.
20 : !> The base class (ty_optical_props) encapsulates only this spectral discretization and must be initialized
21 : !> with the spectral information before use.
22 : !>
23 : !> Optical properties may be represented as arrays with dimensions ncol, nlay, ngpt
24 : !> (abstract class ty_optical_props_arry).
25 : !> The type holds arrays depending on how much information is needed
26 : !> There are three possibilites
27 : !>
28 : !> - ty_optical_props_1scl holds absorption optical depth tau, used in calculations accounting for extinction and emission
29 : !> - ty_optical_props_2str holds extincion optical depth tau, single-scattering albedo ssa, and
30 : !> asymmetry parameter g. These fields are what's needed for two-stream calculations.
31 : !> - ty_optical_props_nstr holds extincion optical depth tau, single-scattering albedo ssa, and
32 : !> phase function moments p with leading dimension nmom. These fields are what's needed for multi-stream calculations.
33 : !>
34 : !> These classes must be allocated before use. Initialization and allocation can be combined.
35 : !> The classes have a validate() function that checks all arrays for valid values (e.g. tau > 0.)
36 : !>
37 : !> Optical properties can be delta-scaled (though this is currently implemented only for two-stream arrays)
38 : !>
39 : !> Optical properties can increment or "add themselves to" a set of properties represented with arrays
40 : !> as long as both sets have the same underlying band structure. Properties defined by band
41 : !> may be added to properties defined by g-point; the same value is assumed for all g-points with each band.
42 : !>
43 : !> Subsets of optical properties held as arrays may be extracted along the column dimension.
44 : !>
45 : ! Making the documentation below visible in ford, by replace ! with !>, messes up Markdown parsing above
46 : !@note
47 : !example of a note with links to other modules and variables
48 : !
49 : ! 1. [[mo_rte_config(module):check_extents(variable)]] in module [[mo_rte_config]]
50 : !>
51 : !@endnote
52 : !> -------------------------------------------------------------------------------------------------
53 : module mo_optical_props
54 : use mo_rte_kind, only: wp
55 : use mo_rte_config, only: check_extents, check_values
56 : use mo_rte_util_array_validation, &
57 : only: any_vals_less_than, any_vals_outside, extents_are
58 : use mo_optical_props_kernels, only: &
59 : increment_1scalar_by_1scalar, increment_1scalar_by_2stream, increment_1scalar_by_nstream, &
60 : increment_2stream_by_1scalar, increment_2stream_by_2stream, increment_2stream_by_nstream, &
61 : increment_nstream_by_1scalar, increment_nstream_by_2stream, increment_nstream_by_nstream, &
62 : inc_1scalar_by_1scalar_bybnd, inc_1scalar_by_2stream_bybnd, inc_1scalar_by_nstream_bybnd, &
63 : inc_2stream_by_1scalar_bybnd, inc_2stream_by_2stream_bybnd, inc_2stream_by_nstream_bybnd, &
64 : inc_nstream_by_1scalar_bybnd, inc_nstream_by_2stream_bybnd, inc_nstream_by_nstream_bybnd, &
65 : delta_scale_2str_kernel, &
66 : extract_subset
67 : implicit none
68 : private
69 : integer, parameter, public :: name_len = 32
70 : ! -------------------------------------------------------------------------------------------------
71 : !
72 : !> ### Base class for optical properties
73 : !> Describes the spectral discretization including the wavenumber limits
74 : !> of each band (spectral region) and the mapping between g-points and bands
75 : !
76 : ! -------------------------------------------------------------------------------------------------
77 : type, public :: ty_optical_props
78 : integer, dimension(:,:), allocatable, private :: band2gpt ! (begin g-point, end g-point) = band2gpt(2,band)
79 : integer, dimension(:), allocatable, private :: gpt2band ! band = gpt2band(g-point)
80 : real(wp), dimension(:,:), allocatable, private :: band_lims_wvn ! (upper and lower wavenumber by band) = band_lims_wvn(2,band)
81 : character(len=name_len), private :: name = ""
82 : contains
83 : generic, public :: init => init_base, init_base_from_copy
84 : procedure, private :: init_base
85 : procedure, private :: init_base_from_copy
86 : procedure, public :: is_initialized => is_initialized_base
87 : procedure, private :: is_initialized_base
88 : procedure, public :: finalize => finalize_base
89 : procedure, private :: finalize_base
90 : procedure, public :: get_nband
91 : procedure, public :: get_ngpt
92 : procedure, public :: get_gpoint_bands
93 : procedure, public :: convert_band2gpt
94 : procedure, public :: convert_gpt2band
95 : procedure, public :: get_band_lims_gpoint
96 : procedure, public :: get_band_lims_wavenumber
97 : procedure, public :: get_band_lims_wavelength
98 : procedure, public :: bands_are_equal
99 : procedure, public :: gpoints_are_equal
100 : procedure, public :: expand
101 : procedure, public :: set_name
102 : procedure, public :: get_name
103 : end type ty_optical_props
104 : !----------------------------------------------------------------------------------------
105 : !
106 : !>### Optical properties as arrays, normally dimensioned ncol, nlay, ngpt/nbnd
107 : !> The abstract base class for arrays defines what procedures will be available
108 : !
109 : ! -------------------------------------------------------------------------------------------------
110 : type, extends(ty_optical_props), abstract, public :: ty_optical_props_arry
111 : real(wp), dimension(:,:,:), allocatable :: tau !! optical depth (ncol, nlay, ngpt)
112 : contains
113 : procedure, public :: get_ncol
114 : procedure, public :: get_nlay
115 : !>
116 : !> Increment another set of values
117 : !>
118 : procedure, public :: increment
119 :
120 : !>
121 : !> Deferred procedures -- each must be implemented in each child class with
122 : !> arguments following the abstract interface (defined below)
123 : !>
124 : procedure(validate_abstract), deferred, public :: validate
125 : procedure(delta_scale_abstract), deferred, public :: delta_scale
126 : procedure(subset_range_abstract), deferred, public :: get_subset
127 : end type ty_optical_props_arry
128 : !
129 : ! Interfaces for the methods to be implemented
130 : !
131 : abstract interface
132 : !>
133 : !> Validation function looks only at internal data
134 : !>
135 : function validate_abstract(this) result(err_message)
136 : import ty_optical_props_arry
137 : class(ty_optical_props_arry), intent(in) :: this
138 : character(len=128) :: err_message
139 : end function validate_abstract
140 :
141 : !>
142 : !> Delta-scaling
143 : !>
144 : function delta_scale_abstract(this, for) result(err_message)
145 : import ty_optical_props_arry
146 : import wp
147 : class(ty_optical_props_arry), intent(inout) :: this
148 : real(wp), dimension(:,:,:), optional, &
149 : intent(in ) :: for
150 : !> Forward scattering fraction; g**2 if not provided
151 : character(len=128) :: err_message
152 : end function delta_scale_abstract
153 :
154 : !>
155 : !> Subsetting -- currently there are only routines with start col and count
156 : !>
157 : function subset_range_abstract(full, start, n, subset) result(err_message)
158 : import ty_optical_props_arry
159 : class(ty_optical_props_arry), intent(inout) :: full
160 : integer, intent(in ) :: start, n
161 : class(ty_optical_props_arry), intent(inout) :: subset
162 : character(128) :: err_message
163 : end function subset_range_abstract
164 : end interface
165 : !>----------------------------------------------------------------------------------------
166 : !>
167 : !> ty_optical_props_arry represents optical properties as arrays with dimensions
168 : !> column, height, spectral point
169 : !> - Class 1-scalar only (extinction) optical depth
170 : !> - Class two-stream adds arrays for single scattering albedo ssa and
171 : !> asymmetry parameter needed in two-stream methods
172 : !> - Class n-stream adds arrays for single scattering albedo ssa and
173 : !> phase function moments (index 1 = g) for use with discrete ordinate methods
174 : !>
175 : !> -------------------------------------------------------------------------------------------------
176 : type, public, extends(ty_optical_props_arry) :: ty_optical_props_1scl
177 : contains
178 : procedure, public :: validate => validate_1scalar
179 : procedure, public :: get_subset => subset_1scl_range
180 : procedure, public :: delta_scale => delta_scale_1scl
181 : procedure, public :: finalize_1scl
182 :
183 : procedure, private :: alloc_only_1scl
184 : procedure, private :: init_and_alloc_1scl
185 : procedure, private :: copy_and_alloc_1scl
186 : generic, public :: alloc_1scl => alloc_only_1scl, init_and_alloc_1scl, copy_and_alloc_1scl
187 : end type ty_optical_props_1scl
188 :
189 : ! --- 2 stream ------------------------------------------------------------------------
190 : type, public, extends(ty_optical_props_arry) :: ty_optical_props_2str
191 : real(wp), dimension(:,:,:), allocatable :: ssa !! single-scattering albedo (ncol, nlay, ngpt)
192 : real(wp), dimension(:,:,:), allocatable :: g !! asymmetry parameter (ncol, nlay, ngpt)
193 : contains
194 : procedure, public :: validate => validate_2stream
195 : procedure, public :: get_subset => subset_2str_range
196 : procedure, public :: delta_scale => delta_scale_2str
197 : procedure, public :: finalize_2str
198 :
199 : procedure, private :: alloc_only_2str
200 : procedure, private :: init_and_alloc_2str
201 : procedure, private :: copy_and_alloc_2str
202 : generic, public :: alloc_2str => alloc_only_2str, init_and_alloc_2str, copy_and_alloc_2str
203 : end type ty_optical_props_2str
204 :
205 : ! --- n stream ------------------------------------------------------------------------
206 : type, public, extends(ty_optical_props_arry) :: ty_optical_props_nstr
207 : real(wp), dimension(:,:,:), allocatable :: ssa !! single-scattering albedo (ncol, nlay, ngpt)
208 : real(wp), dimension(:,:,:,:), allocatable :: p !! phase-function moments (nmom, ncol, nlay, ngpt)
209 : contains
210 : procedure, public :: validate => validate_nstream
211 : procedure, public :: get_subset => subset_nstr_range
212 : procedure, public :: delta_scale => delta_scale_nstr
213 : procedure, public :: get_nmom
214 : procedure, public :: finalize_nstr
215 :
216 : procedure, private :: alloc_only_nstr
217 : procedure, private :: init_and_alloc_nstr
218 : procedure, private :: copy_and_alloc_nstr
219 : generic, public :: alloc_nstr => alloc_only_nstr, init_and_alloc_nstr, copy_and_alloc_nstr
220 : end type ty_optical_props_nstr
221 : ! -------------------------------------------------------------------------------------------------
222 : contains
223 : ! -------------------------------------------------------------------------------------------------
224 : !
225 : ! Routines for the base class: initialization, validity checking, finalization
226 : !
227 : !> -------------------------------------------------------------------------------------------------
228 : !>
229 : !> Base class: Initialization
230 : !> Values are assumed to be defined in bands a mapping between bands and g-points is provided
231 : !>
232 : !> -------------------------------------------------------------------------------------------------
233 4155405 : function init_base(this, band_lims_wvn, band_lims_gpt, name) result(err_message)
234 : class(ty_optical_props), intent(inout) :: this
235 : real(wp), dimension(:,:), intent(in ) :: band_lims_wvn
236 : integer, dimension(:,:), &
237 : optional, intent(in ) :: band_lims_gpt
238 : character(len=*), optional, intent(in ) :: name
239 : character(len = 128) :: err_message
240 :
241 : integer :: iband
242 4155405 : integer, dimension(2, size(band_lims_wvn, 2)) :: band_lims_gpt_lcl
243 : ! -------------------------
244 : !
245 : ! Error checking -- are the arrays the size we expect, contain positive values?
246 : !
247 4155405 : err_message = ""
248 4155405 : if(size(band_lims_wvn,1) /= 2) &
249 0 : err_message = "optical_props%init(): band_lims_wvn 1st dim should be 2"
250 4155405 : if (check_values) then
251 4155405 : if(any_vals_less_than(band_lims_wvn, 0._wp) ) &
252 0 : err_message = "optical_props%init(): band_lims_wvn has values < 0., respectively"
253 : end if
254 4155405 : if(err_message /="") return
255 4155405 : if(present(band_lims_gpt)) then
256 3020006 : if (check_extents) then
257 3020006 : if(.not. extents_are(band_lims_gpt, 2, size(band_lims_wvn,2))) &
258 0 : err_message = "optical_props%init(): band_lims_gpt size inconsistent with band_lims_wvn"
259 : end if
260 3020006 : if (check_values) then
261 143299922 : if(any(band_lims_gpt < 1) ) &
262 0 : err_message = "optical_props%init(): band_lims_gpt has values < 1"
263 : end if
264 3020006 : if(err_message /= "") return
265 :
266 143299922 : band_lims_gpt_lcl(:,:) = band_lims_gpt(:,:)
267 : else
268 : !
269 : ! Assume that values are defined by band, one g-point per band
270 : !
271 18523257 : do iband = 1, size(band_lims_wvn, 2)
272 53298973 : band_lims_gpt_lcl(1:2,iband) = iband
273 : end do
274 : end if
275 : !
276 : ! Assignment
277 : !
278 4155405 : if(allocated(this%band2gpt )) deallocate(this%band2gpt)
279 4155405 : if(allocated(this%band_lims_wvn)) deallocate(this%band_lims_wvn)
280 0 : allocate(this%band2gpt (2,size(band_lims_wvn,2)), &
281 20777025 : this%band_lims_wvn(2,size(band_lims_wvn,2)))
282 200754300 : this%band2gpt = band_lims_gpt_lcl
283 200754300 : this%band_lims_wvn = band_lims_wvn
284 4155405 : if(present(name)) this%name = trim(name)
285 :
286 : !
287 : ! Make a map between g-points and bands
288 : ! Efficient only when g-point indexes start at 1 and are contiguous.
289 : !
290 4155405 : if(allocated(this%gpt2band)) deallocate(this%gpt2band)
291 204909705 : allocate(this%gpt2band(maxval(band_lims_gpt_lcl)))
292 68303235 : do iband=1,size(band_lims_gpt_lcl,dim=2)
293 459770869 : this%gpt2band(band_lims_gpt_lcl(1,iband):band_lims_gpt_lcl(2,iband)) = iband
294 : end do
295 : end function init_base
296 : !-------------------------------------------------------------------------------------------------
297 746136 : function init_base_from_copy(this, spectral_desc) result(err_message)
298 : class(ty_optical_props), intent(inout) :: this
299 : class(ty_optical_props), intent(in ) :: spectral_desc
300 : character(len = 128) :: err_message
301 :
302 746136 : if(.not. spectral_desc%is_initialized()) then
303 0 : err_message = "optical_props%init(): can't initialize based on un-initialized input"
304 0 : return
305 : else
306 : err_message = this%init(spectral_desc%get_band_lims_wavenumber(), &
307 746136 : spectral_desc%get_band_lims_gpoint())
308 : end if
309 : end function init_base_from_copy
310 : !>-------------------------------------------------------------------------------------------------
311 : !>
312 : !> Base class: return true if initialized, false otherwise
313 : !>
314 : !> -------------------------------------------------------------------------------------------------
315 17134072852 : pure function is_initialized_base(this)
316 : class(ty_optical_props), intent(in) :: this
317 : logical :: is_initialized_base
318 :
319 17134072852 : is_initialized_base = allocated(this%band2gpt)
320 17134072852 : end function is_initialized_base
321 : !>-------------------------------------------------------------------------------------------------
322 : !>
323 : !> Base class: finalize (deallocate memory)
324 : !>
325 : !> -------------------------------------------------------------------------------------------------
326 9699768 : subroutine finalize_base(this)
327 : class(ty_optical_props), intent(inout) :: this
328 :
329 9699768 : if(allocated(this%band2gpt)) deallocate(this%band2gpt)
330 9699768 : if(allocated(this%gpt2band)) deallocate(this%gpt2band)
331 9699768 : if(allocated(this%band_lims_wvn)) &
332 3406197 : deallocate(this%band_lims_wvn)
333 9699768 : this%name = ""
334 9699768 : end subroutine finalize_base
335 : ! ------------------------------------------------------------------------------------------
336 : !
337 : ! Routines for array classes: initialization, allocation, and finalization
338 : ! Initialization and allocation can be combined by supplying either
339 : !
340 : !> ------------------------------------------------------------------------------------------
341 : !>
342 : !> Straight allocation routines
343 : !>
344 : !> --- 1 scalar ------------------------------------------------------------------------
345 2238408 : function alloc_only_1scl(this, ncol, nlay) result(err_message)
346 : class(ty_optical_props_1scl) :: this
347 : integer, intent(in) :: ncol, nlay
348 : character(len=128) :: err_message
349 :
350 2238408 : err_message = ""
351 2238408 : if(.not. this%is_initialized()) then
352 0 : err_message = "optical_props%alloc: spectral discretization hasn't been provided"
353 0 : return
354 : end if
355 6715224 : if(any([ncol, nlay] <= 0)) then
356 0 : err_message = "optical_props%alloc: must provide positive extents for ncol, nlay"
357 : else
358 2238408 : if(allocated(this%tau)) deallocate(this%tau)
359 11192040 : allocate(this%tau(ncol,nlay,this%get_ngpt()))
360 : end if
361 : end function alloc_only_1scl
362 :
363 : !> --- 2 stream ------------------------------------------------------------------------
364 1167789 : function alloc_only_2str(this, ncol, nlay) result(err_message)
365 : class(ty_optical_props_2str) :: this
366 : integer, intent(in) :: ncol, nlay
367 : character(len=128) :: err_message
368 :
369 1167789 : err_message = ""
370 1167789 : if(.not. this%is_initialized()) &
371 0 : err_message = "optical_props%alloc: spectral discretization hasn't been provided"
372 3503367 : if(any([ncol, nlay] <= 0)) &
373 0 : err_message = "optical_props%alloc: must provide positive extents for ncol, nlay"
374 1167789 : if(err_message /= "") return
375 :
376 1167789 : if(allocated(this%tau)) deallocate(this%tau)
377 5838945 : allocate(this%tau(ncol,nlay,this%get_ngpt()))
378 1167789 : if(allocated(this%ssa)) deallocate(this%ssa)
379 1167789 : if(allocated(this%g )) deallocate(this%g )
380 10510101 : allocate(this%ssa(ncol,nlay,this%get_ngpt()), this%g(ncol,nlay,this%get_ngpt()))
381 : end function alloc_only_2str
382 :
383 : !> --- n stream ------------------------------------------------------------------------
384 0 : function alloc_only_nstr(this, nmom, ncol, nlay) result(err_message)
385 : class(ty_optical_props_nstr) :: this
386 : integer, intent(in) :: nmom ! number of moments
387 : integer, intent(in) :: ncol, nlay
388 : character(len=128) :: err_message
389 :
390 0 : err_message = ""
391 0 : if(.not. this%is_initialized()) &
392 0 : err_message = "optical_props%alloc: spectral discretization hasn't been provided"
393 0 : if(any([ncol, nlay] <= 0)) &
394 0 : err_message = "optical_props%alloc: must provide positive extents for ncol, nlay"
395 0 : if(err_message /= "") return
396 :
397 0 : if(allocated(this%tau)) deallocate(this%tau)
398 0 : allocate(this%tau(ncol,nlay,this%get_ngpt()))
399 0 : if(allocated(this%ssa)) deallocate(this%ssa)
400 0 : if(allocated(this%p )) deallocate(this%p )
401 0 : allocate(this%ssa(ncol,nlay,this%get_ngpt()), this%p(nmom,ncol,nlay,this%get_ngpt()))
402 : end function alloc_only_nstr
403 : ! ------------------------------------------------------------------------------------------
404 : !
405 : ! Combined allocation/initialization routines
406 : !
407 : !> ------------------------------------------------------------------------------------------
408 : !>
409 : !> Initialization by specifying band limits and possibly g-point/band mapping
410 : !>
411 : !> ---------------------------------------------------------------------------
412 746136 : function init_and_alloc_1scl(this, ncol, nlay, band_lims_wvn, band_lims_gpt, name) result(err_message)
413 : class(ty_optical_props_1scl) :: this
414 : integer, intent(in) :: ncol, nlay
415 : real(wp), dimension(:,:), intent(in) :: band_lims_wvn
416 : integer, dimension(:,:), &
417 : optional, intent(in) :: band_lims_gpt
418 : character(len=*), optional, intent(in) :: name
419 : character(len=128) :: err_message
420 :
421 : err_message = this%ty_optical_props%init(band_lims_wvn, &
422 2238408 : band_lims_gpt, name)
423 746136 : if(err_message /= "") return
424 746136 : err_message = this%alloc_1scl(ncol, nlay)
425 : end function init_and_alloc_1scl
426 : ! ---------------------------------------------------------------------------
427 389263 : function init_and_alloc_2str(this, ncol, nlay, band_lims_wvn, band_lims_gpt, name) result(err_message)
428 : class(ty_optical_props_2str) :: this
429 : integer, intent(in) :: ncol, nlay
430 : real(wp), dimension(:,:), intent(in) :: band_lims_wvn
431 : integer, dimension(:,:), &
432 : optional, intent(in) :: band_lims_gpt
433 : character(len=*), optional, intent(in) :: name
434 : character(len=128) :: err_message
435 :
436 : err_message = this%ty_optical_props%init(band_lims_wvn, &
437 1167789 : band_lims_gpt, name)
438 389263 : if(err_message /= "") return
439 389263 : err_message = this%alloc_2str(ncol, nlay)
440 : end function init_and_alloc_2str
441 : ! ---------------------------------------------------------------------------
442 0 : function init_and_alloc_nstr(this, nmom, ncol, nlay, band_lims_wvn, band_lims_gpt, name) result(err_message)
443 : class(ty_optical_props_nstr) :: this
444 : integer, intent(in) :: nmom, ncol, nlay
445 : real(wp), dimension(:,:), intent(in) :: band_lims_wvn
446 : integer, dimension(:,:), &
447 : optional, intent(in) :: band_lims_gpt
448 : character(len=*), optional, intent(in) :: name
449 : character(len=128) :: err_message
450 :
451 : err_message = this%ty_optical_props%init(band_lims_wvn, &
452 0 : band_lims_gpt, name)
453 0 : if(err_message /= "") return
454 0 : err_message = this%alloc_nstr(nmom, ncol, nlay)
455 : end function init_and_alloc_nstr
456 : !>-------------------------------------------------------------------------------------------------
457 : !>
458 : !> Initialization from an existing spectral discretization/ty_optical_props
459 : !>
460 : !>-------------------------------------------------------------------------------------------------
461 1492272 : function copy_and_alloc_1scl(this, ncol, nlay, spectral_desc, name) result(err_message)
462 : class(ty_optical_props_1scl) :: this
463 : integer, intent(in) :: ncol, nlay
464 : class(ty_optical_props ), intent(in) :: spectral_desc
465 : character(len=*), optional, intent(in) :: name
466 : character(len=128) :: err_message
467 :
468 1492272 : err_message = ""
469 1492272 : if(this%ty_optical_props%is_initialized()) call this%ty_optical_props%finalize()
470 : err_message = this%ty_optical_props%init(spectral_desc%get_band_lims_wavenumber(), &
471 2984544 : spectral_desc%get_band_lims_gpoint(), name)
472 1492272 : if(err_message /= "") return
473 1492272 : err_message = this%alloc_1scl(ncol, nlay)
474 : end function copy_and_alloc_1scl
475 : ! ---------------------------------------------------------------------------
476 778526 : function copy_and_alloc_2str(this, ncol, nlay, spectral_desc, name) result(err_message)
477 : class(ty_optical_props_2str) :: this
478 : integer, intent(in) :: ncol, nlay
479 : class(ty_optical_props ), intent(in) :: spectral_desc
480 : character(len=*), optional, intent(in) :: name
481 : character(len=128) :: err_message
482 :
483 778526 : err_message = ""
484 778526 : if(this%ty_optical_props%is_initialized()) call this%ty_optical_props%finalize()
485 : err_message = this%ty_optical_props%init(spectral_desc%get_band_lims_wavenumber(), &
486 1557052 : spectral_desc%get_band_lims_gpoint(), name)
487 778526 : if(err_message /= "") return
488 778526 : err_message = this%alloc_2str(ncol, nlay)
489 : end function copy_and_alloc_2str
490 : ! ---------------------------------------------------------------------------
491 0 : function copy_and_alloc_nstr(this, nmom, ncol, nlay, spectral_desc, name) result(err_message)
492 : class(ty_optical_props_nstr) :: this
493 : integer, intent(in) :: nmom, ncol, nlay
494 : class(ty_optical_props ), intent(in) :: spectral_desc
495 : character(len=*), optional, intent(in) :: name
496 : character(len=128) :: err_message
497 :
498 0 : err_message = ""
499 0 : if(this%ty_optical_props%is_initialized()) call this%ty_optical_props%finalize()
500 : err_message = this%ty_optical_props%init(spectral_desc%get_band_lims_wavenumber(), &
501 0 : spectral_desc%get_band_lims_gpoint(), name)
502 0 : if(err_message /= "") return
503 0 : err_message = this%alloc_nstr(nmom, ncol, nlay)
504 : end function copy_and_alloc_nstr
505 : !> ------------------------------------------------------------------------------------------
506 : !>
507 : !> Finalize routines
508 : !>
509 : !> ------------------------------------------------------------------------------------------
510 0 : function finalize_1scl(this) result(err_message)
511 : class(ty_optical_props_1scl) :: this
512 : character(len=128) :: err_message
513 :
514 0 : if(allocated(this%tau)) deallocate(this%tau)
515 0 : err_message = ""
516 0 : end function finalize_1scl
517 : ! ---------------------------------------------------------------------------
518 0 : function finalize_2str(this) result(err_message)
519 : class(ty_optical_props_2str) :: this
520 : character(len=128) :: err_message
521 :
522 0 : if(allocated(this%tau)) deallocate(this%tau)
523 0 : if(allocated(this%ssa)) deallocate(this%ssa)
524 0 : if(allocated(this%g )) deallocate(this%g )
525 0 : err_message = ""
526 0 : end function finalize_2str
527 : ! ---------------------------------------------------------------------------
528 0 : function finalize_nstr(this) result(err_message)
529 : class(ty_optical_props_nstr) :: this
530 : character(len=128) :: err_message
531 :
532 0 : if(allocated(this%tau)) deallocate(this%tau)
533 0 : if(allocated(this%ssa)) deallocate(this%ssa)
534 0 : if(allocated(this%p )) deallocate(this%p )
535 0 : err_message = ""
536 0 : end function finalize_nstr
537 : ! ------------------------------------------------------------------------------------------
538 : !
539 : ! Routines for array classes: delta-scaling, validation (ensuring all values can be used )
540 : !
541 : !> ------------------------------------------------------------------------------------------
542 : !> --- delta scaling
543 : !> ------------------------------------------------------------------------------------------
544 0 : function delta_scale_1scl(this, for) result(err_message)
545 : class(ty_optical_props_1scl), intent(inout) :: this
546 : real(wp), dimension(:,:,:), optional, &
547 : intent(in ) :: for
548 : character(128) :: err_message
549 : !
550 : ! Nothing to do
551 : !
552 0 : err_message = ""
553 0 : end function delta_scale_1scl
554 : ! ------------------------------------------------------------------------------------------
555 389263 : function delta_scale_2str(this, for) result(err_message)
556 : class(ty_optical_props_2str), intent(inout) :: this
557 : real(wp), dimension(:,:,:), optional, &
558 : intent(in ) :: for
559 : !! Forward scattering fraction; g**2 if not provided
560 : character(128) :: err_message
561 :
562 : integer :: ncol, nlay, ngpt
563 : ! --------------------------------
564 389263 : ncol = this%get_ncol()
565 389263 : nlay = this%get_nlay()
566 389263 : ngpt = this%get_ngpt()
567 389263 : err_message = ""
568 :
569 389263 : if(present(for)) then
570 0 : if (check_extents) then
571 0 : if(.not. extents_are(for, ncol, nlay, ngpt)) then
572 0 : err_message = "delta_scale: dimension of 'for' don't match optical properties arrays"
573 0 : return
574 : end if
575 : end if
576 0 : if (check_values) then
577 0 : if(any_vals_outside(for, 0._wp, 1._wp)) then
578 0 : err_message = "delta_scale: values of 'for' out of bounds [0,1]"
579 0 : return
580 : end if
581 : end if
582 0 : call delta_scale_2str_kernel(ncol, nlay, ngpt, this%tau, this%ssa, this%g, for)
583 : else
584 389263 : call delta_scale_2str_kernel(ncol, nlay, ngpt, this%tau, this%ssa, this%g)
585 : end if
586 :
587 : end function delta_scale_2str
588 : ! ------------------------------------------------------------------------------------------
589 0 : function delta_scale_nstr(this, for) result(err_message)
590 : class(ty_optical_props_nstr), intent(inout) :: this
591 : real(wp), dimension(:,:,:), optional, &
592 : intent(in ) :: for
593 : character(128) :: err_message
594 :
595 0 : err_message = 'delta_scale_nstr: Not yet implemented'
596 0 : end function delta_scale_nstr
597 : !> ------------------------------------------------------------------------------------------
598 : !>
599 : !> --- Validation
600 : !>
601 : !> ------------------------------------------------------------------------------------------
602 2984544 : function validate_1scalar(this) result(err_message)
603 : class(ty_optical_props_1scl), intent(in) :: this
604 : character(len=128) :: err_message
605 :
606 2984544 : err_message = ''
607 2984544 : if(.not. allocated(this%tau)) then
608 0 : err_message = "validate: tau not allocated/initialized"
609 0 : return
610 : end if
611 2984544 : if (check_values) then
612 2984544 : if(any_vals_less_than(this%tau, 0._wp)) &
613 0 : err_message = "validate: tau values out of range"
614 : end if
615 2984544 : if(len_trim(err_message) > 0 .and. len_trim(this%get_name()) > 0) &
616 2984544 : err_message = trim(this%get_name()) // ': ' // trim(err_message)
617 :
618 : end function validate_1scalar
619 : ! ------------------------------------------------------------------------------------------
620 1167789 : function validate_2stream(this) result(err_message)
621 : class(ty_optical_props_2str), intent(in) :: this
622 : character(len=128) :: err_message
623 :
624 : integer :: d1, d2, d3
625 :
626 1167789 : err_message = ''
627 : !
628 : ! Array allocation status, sizing
629 : !
630 1167789 : if(check_extents) then
631 4671156 : if(.not. all([allocated(this%tau), allocated(this%ssa), allocated(this%g)])) then
632 0 : err_message = "validate: arrays not allocated/initialized"
633 0 : return
634 : end if
635 1167789 : d1 = size(this%tau, 1)
636 1167789 : d2 = size(this%tau, 2)
637 1167789 : d3 = size(this%tau, 3)
638 1167789 : if(.not. extents_are(this%ssa, d1, d2, d3) .or. &
639 : .not. extents_are(this%g , d1, d2, d3)) &
640 0 : err_message = "validate: arrays not sized consistently"
641 : end if
642 : !
643 : ! Valid values
644 : !
645 1167789 : if (check_values) then
646 1167789 : if(any_vals_less_than(this%tau, 0._wp)) &
647 0 : err_message = "validate: tau values out of range"
648 1167789 : if(any_vals_outside (this%ssa, 0._wp, 1._wp)) &
649 0 : err_message = "validate: ssa values out of range"
650 1167789 : if(any_vals_outside (this%g , -1._wp, 1._wp)) &
651 0 : err_message = "validate: g values out of range"
652 : end if
653 :
654 1167789 : if(len_trim(err_message) > 0 .and. len_trim(this%get_name()) > 0) &
655 1167789 : err_message = trim(this%get_name()) // ': ' // trim(err_message)
656 :
657 : end function validate_2stream
658 :
659 : ! ------------------------------------------------------------------------------------------
660 0 : function validate_nstream(this) result(err_message)
661 : class(ty_optical_props_nstr), intent(in) :: this
662 : character(len=128) :: err_message
663 :
664 : integer :: d1, d2, d3, d4
665 :
666 0 : err_message = ''
667 : !
668 : ! Array allocation status, sizing
669 : !
670 0 : if(.not. all([allocated(this%tau), allocated(this%ssa), allocated(this%p)])) then
671 0 : err_message = "validate: arrays not allocated/initialized"
672 0 : return
673 : end if
674 0 : d1 = size(this%tau, 1)
675 0 : d2 = size(this%tau, 2)
676 0 : d3 = size(this%tau, 3)
677 0 : d4 = size(this%p, 1)
678 0 : if (check_extents) then
679 0 : if(.not. extents_are(this%ssa, d1, d2, d3) .or. &
680 : .not. extents_are(this%p , d4, d1, d2, d3)) &
681 0 : err_message = "validate: arrays not sized consistently"
682 : end if
683 : !
684 : ! Valid values
685 : !
686 0 : if (check_values) then
687 0 : if(any_vals_less_than(this%tau, 0._wp)) &
688 0 : err_message = "validate: tau values out of range"
689 0 : if(any_vals_outside (this%ssa, 0._wp, 1._wp)) &
690 0 : err_message = "validate: ssa values out of range"
691 0 : if(any_vals_outside (this%p(1,:,:,:), -1._wp, 1._wp)) &
692 0 : err_message = "validate: p(1,:,:,:) = g values out of range"
693 : end if
694 :
695 0 : if(len_trim(err_message) > 0 .and. len_trim(this%get_name()) > 0) &
696 0 : err_message = trim(this%get_name()) // ': ' // trim(err_message)
697 : end function validate_nstream
698 :
699 : ! ------------------------------------------------------------------------------------------
700 : !
701 : ! Routines for array classes: subsetting of optical properties arrays along x (col) direction
702 : !
703 : ! Allocate class, then arrays; copy. Could probably be more efficient if
704 : ! classes used pointers internally.
705 : !
706 : ! This set takes start position and number as scalars
707 : !
708 : ! ------------------------------------------------------------------------------------------
709 :
710 0 : function subset_1scl_range(full, start, n, subset) result(err_message)
711 : class(ty_optical_props_1scl), intent(inout) :: full
712 : integer, intent(in ) :: start, n
713 : class(ty_optical_props_arry), intent(inout) :: subset
714 : character(128) :: err_message
715 :
716 : integer :: ncol, nlay, ngpt, nmom
717 :
718 0 : err_message = ""
719 0 : if(.not. full%is_initialized()) then
720 0 : err_message = "optical_props%subset: Asking for a subset of uninitialized data"
721 0 : return
722 : end if
723 0 : ncol = full%get_ncol()
724 0 : nlay = full%get_nlay()
725 0 : ngpt = full%get_ngpt()
726 0 : if(start < 1 .or. start + n-1 > full%get_ncol()) &
727 0 : err_message = "optical_props%subset: Asking for columns outside range"
728 0 : if(err_message /= "") return
729 :
730 0 : if(subset%is_initialized()) call subset%finalize()
731 0 : err_message = subset%init(full)
732 : ! Seems like the deallocation statements should be needed under Fortran 2003
733 : ! but Intel compiler doesn't run without them
734 0 : if(allocated(subset%tau)) deallocate(subset%tau)
735 : select type (subset)
736 : class is (ty_optical_props_1scl)
737 0 : err_message = subset%alloc_1scl(n, nlay)
738 0 : if(err_message /= "") return
739 : class is (ty_optical_props_2str)
740 0 : if(allocated(subset%ssa)) deallocate(subset%ssa)
741 0 : if(allocated(subset%g )) deallocate(subset%g )
742 0 : err_message = subset%alloc_2str(n, nlay)
743 0 : if(err_message /= "") return
744 0 : subset%ssa(1:n,:,:) = 0._wp
745 0 : subset%g (1:n,:,:) = 0._wp
746 : class is (ty_optical_props_nstr)
747 0 : if(allocated(subset%ssa)) deallocate(subset%ssa)
748 0 : if(allocated(subset%p )) then
749 0 : nmom = subset%get_nmom()
750 0 : deallocate(subset%p )
751 : else
752 0 : nmom = 1
753 : end if
754 0 : err_message = subset%alloc_nstr(nmom, n, nlay)
755 0 : if(err_message /= "") return
756 0 : subset%ssa(1:n,:,:) = 0._wp
757 0 : subset%p(:,1:n,:,:) = 0._wp
758 : end select
759 0 : call extract_subset(ncol, nlay, ngpt, full%tau, start, start+n-1, subset%tau)
760 :
761 : end function subset_1scl_range
762 : ! ------------------------------------------------------------------------------------------
763 0 : function subset_2str_range(full, start, n, subset) result(err_message)
764 : class(ty_optical_props_2str), intent(inout) :: full
765 : integer, intent(in ) :: start, n
766 : class(ty_optical_props_arry), intent(inout) :: subset
767 : character(128) :: err_message
768 :
769 : integer :: ncol, nlay, ngpt, nmom
770 :
771 0 : err_message = ""
772 0 : if(.not. full%is_initialized()) then
773 0 : err_message = "optical_props%subset: Asking for a subset of uninitialized data"
774 0 : return
775 : end if
776 0 : ncol = full%get_ncol()
777 0 : nlay = full%get_nlay()
778 0 : ngpt = full%get_ngpt()
779 0 : if(start < 1 .or. start + n-1 > full%get_ncol()) &
780 0 : err_message = "optical_props%subset: Asking for columns outside range"
781 0 : if(err_message /= "") return
782 :
783 0 : if(subset%is_initialized()) call subset%finalize()
784 0 : err_message = subset%init(full)
785 : select type (subset)
786 : class is (ty_optical_props_1scl)
787 0 : err_message = subset%alloc_1scl(n, nlay)
788 0 : if(err_message /= "") return
789 0 : call extract_subset(ncol, nlay, ngpt, full%tau, full%ssa, start, start+n-1, subset%tau)
790 : class is (ty_optical_props_2str)
791 0 : if(allocated(subset%ssa)) deallocate(subset%ssa)
792 0 : if(allocated(subset%g )) deallocate(subset%g )
793 0 : err_message = subset%alloc_2str(n, nlay)
794 0 : if(err_message /= "") return
795 0 : call extract_subset(ncol, nlay, ngpt, full%tau, start, start+n-1, subset%tau)
796 0 : call extract_subset(ncol, nlay, ngpt, full%ssa, start, start+n-1, subset%ssa)
797 0 : call extract_subset(ncol, nlay, ngpt, full%g , start, start+n-1, subset%g )
798 : class is (ty_optical_props_nstr)
799 0 : if(allocated(subset%ssa)) deallocate(subset%ssa)
800 0 : if(allocated(subset%p )) then
801 0 : nmom = subset%get_nmom()
802 0 : deallocate(subset%p )
803 : else
804 0 : nmom = 1
805 : end if
806 0 : err_message = subset%alloc_nstr(nmom, n, nlay)
807 0 : if(err_message /= "") return
808 0 : call extract_subset(ncol, nlay, ngpt, full%tau, start, start+n-1, subset%tau)
809 0 : call extract_subset(ncol, nlay, ngpt, full%ssa, start, start+n-1, subset%ssa)
810 0 : subset%p(1,1:n,:,:) = full%g (start:start+n-1,:,:)
811 0 : subset%p(2:,:, :,:) = 0._wp
812 : end select
813 :
814 : end function subset_2str_range
815 : ! ------------------------------------------------------------------------------------------
816 0 : function subset_nstr_range(full, start, n, subset) result(err_message)
817 : class(ty_optical_props_nstr), intent(inout) :: full
818 : integer, intent(in ) :: start, n
819 : class(ty_optical_props_arry), intent(inout) :: subset
820 : character(128) :: err_message
821 :
822 : integer :: ncol, nlay, ngpt, nmom
823 :
824 0 : err_message = ""
825 0 : if(.not. full%is_initialized()) then
826 0 : err_message = "optical_props%subset: Asking for a subset of uninitialized data"
827 0 : return
828 : end if
829 0 : ncol = full%get_ncol()
830 0 : nlay = full%get_nlay()
831 0 : ngpt = full%get_ngpt()
832 0 : if(start < 1 .or. start + n-1 > full%get_ncol()) &
833 0 : err_message = "optical_props%subset: Asking for columns outside range"
834 0 : if(err_message /= "") return
835 :
836 0 : if(subset%is_initialized()) call subset%finalize()
837 0 : err_message = subset%init(full)
838 0 : if(allocated(subset%tau)) deallocate(subset%tau)
839 : select type (subset)
840 : class is (ty_optical_props_1scl)
841 0 : err_message = subset%alloc_1scl(n, nlay)
842 0 : if(err_message /= "") return
843 0 : call extract_subset(ncol, nlay, ngpt, full%tau, full%ssa, start, start+n-1, subset%tau)
844 : class is (ty_optical_props_2str)
845 0 : if(allocated(subset%ssa)) deallocate(subset%ssa)
846 0 : if(allocated(subset%g )) deallocate(subset%g )
847 0 : err_message = subset%alloc_2str(n, nlay)
848 0 : if(err_message /= "") return
849 0 : call extract_subset(ncol, nlay, ngpt, full%tau, start, start+n-1, subset%tau)
850 0 : call extract_subset(ncol, nlay, ngpt, full%ssa, start, start+n-1, subset%ssa)
851 0 : subset%g (1:n,:,:) = full%p(1,start:start+n-1,:,:)
852 : class is (ty_optical_props_nstr)
853 0 : if(allocated(subset%ssa)) deallocate(subset%ssa)
854 0 : if(allocated(subset%p )) deallocate(subset%p )
855 0 : err_message = subset%alloc_nstr(nmom, n, nlay)
856 0 : if(err_message /= "") return
857 0 : call extract_subset( ncol, nlay, ngpt, full%tau, start, start+n-1, subset%tau)
858 0 : call extract_subset( ncol, nlay, ngpt, full%ssa, start, start+n-1, subset%ssa)
859 0 : call extract_subset(nmom, ncol, nlay, ngpt, full%p , start, start+n-1, subset%p )
860 : end select
861 : end function subset_nstr_range
862 :
863 : !> ------------------------------------------------------------------------------------------
864 : !>
865 : !> Routines for array classes: incrementing
866 : !> a%increment(b) adds the values of a to b, changing b and leaving a untouched
867 : !>
868 : !> -----------------------------------------------------------------------------------------
869 2270798 : function increment(op_in, op_io) result(err_message)
870 : class(ty_optical_props_arry), intent(in ) :: op_in
871 : class(ty_optical_props_arry), intent(inout) :: op_io
872 : character(128) :: err_message
873 : ! -----
874 : integer :: ncol, nlay, ngpt
875 : ! -----
876 2270798 : err_message = ""
877 2270798 : if(.not. op_in%is_initialized()) &
878 0 : err_message = "ty_optical_props%increment: Incrementing optical properties aren't initialized"
879 2270798 : if(.not. op_in%is_initialized()) &
880 0 : err_message = "ty_optical_props%increment: optical properties to be incremented aren't initialized"
881 2270798 : if(err_message /= "") return
882 :
883 2270798 : ncol = op_io%get_ncol()
884 2270798 : nlay = op_io%get_nlay()
885 2270798 : ngpt = op_io%get_ngpt()
886 2270798 : if(.not. op_in%bands_are_equal(op_io)) &
887 0 : err_message = "ty_optical_props%increment: optical properties objects have different band structures"
888 6812394 : if(.not. all([op_in%get_ncol(), op_in%get_nlay()] == [ncol, nlay])) &
889 0 : err_message = "ty_optical_props%increment: optical properties objects have different ncol and/or nlay"
890 2270798 : if(err_message /= "") return
891 :
892 2270798 : if(op_in%gpoints_are_equal(op_io)) then
893 : !
894 : ! Increment by gpoint
895 : ! (or by band if both op_in and op_io are defined that way)
896 : !
897 : select type (op_io)
898 : class is (ty_optical_props_1scl)
899 : select type (op_in)
900 : class is (ty_optical_props_1scl)
901 : call increment_1scalar_by_1scalar(ncol, nlay, ngpt, &
902 : op_io%tau, &
903 746136 : op_in%tau)
904 : class is (ty_optical_props_2str)
905 : call increment_1scalar_by_2stream(ncol, nlay, ngpt, &
906 : op_io%tau, &
907 0 : op_in%tau, op_in%ssa)
908 :
909 : class is (ty_optical_props_nstr)
910 : call increment_1scalar_by_nstream(ncol, nlay, ngpt, &
911 : op_io%tau, &
912 0 : op_in%tau, op_in%ssa)
913 : end select
914 : class is (ty_optical_props_2str)
915 : select type (op_in)
916 : class is (ty_optical_props_1scl)
917 : call increment_2stream_by_1scalar(ncol, nlay, ngpt, &
918 : op_io%tau, op_io%ssa,&
919 0 : op_in%tau)
920 : class is (ty_optical_props_2str)
921 : call increment_2stream_by_2stream(ncol, nlay, ngpt, &
922 : op_io%tau, op_io%ssa, op_io%g, &
923 389263 : op_in%tau, op_in%ssa, op_in%g)
924 : class is (ty_optical_props_nstr)
925 : call increment_2stream_by_nstream(ncol, nlay, ngpt, op_in%get_nmom(), &
926 : op_io%tau, op_io%ssa, op_io%g, &
927 0 : op_in%tau, op_in%ssa, op_in%p)
928 : end select
929 :
930 : class is (ty_optical_props_nstr)
931 : select type (op_in)
932 : class is (ty_optical_props_1scl)
933 : call increment_nstream_by_1scalar(ncol, nlay, ngpt, &
934 : op_io%tau, op_io%ssa, &
935 0 : op_in%tau)
936 : class is (ty_optical_props_2str)
937 : call increment_nstream_by_2stream(ncol, nlay, ngpt, op_io%get_nmom(), &
938 : op_io%tau, op_io%ssa, op_io%p, &
939 0 : op_in%tau, op_in%ssa, op_in%g)
940 : class is (ty_optical_props_nstr)
941 : call increment_nstream_by_nstream(ncol, nlay, ngpt, op_io%get_nmom(), op_in%get_nmom(), &
942 : op_io%tau, op_io%ssa, op_io%p, &
943 0 : op_in%tau, op_in%ssa, op_in%p)
944 : end select
945 : end select
946 : else
947 : !
948 : ! Values defined by-band will have ngpt() = nband()
949 : ! We can use values by band in op_in to increment op_io
950 : ! Anything else is an error
951 : !
952 1135399 : if(op_in%get_ngpt() /= op_io%get_nband()) then
953 0 : err_message = "ty_optical_props%increment: optical properties objects have incompatible g-point structures"
954 0 : return
955 : end if
956 : !
957 : ! Increment by band
958 : !
959 : select type (op_io)
960 : class is (ty_optical_props_1scl)
961 : select type (op_in)
962 : class is (ty_optical_props_1scl)
963 : call inc_1scalar_by_1scalar_bybnd(ncol, nlay, ngpt, &
964 : op_io%tau, &
965 : op_in%tau, &
966 746136 : op_io%get_nband(), op_io%get_band_lims_gpoint())
967 : class is (ty_optical_props_2str)
968 : call inc_1scalar_by_2stream_bybnd(ncol, nlay, ngpt, &
969 : op_io%tau, &
970 : op_in%tau, op_in%ssa, &
971 0 : op_io%get_nband(), op_io%get_band_lims_gpoint())
972 : class is (ty_optical_props_nstr)
973 : call inc_1scalar_by_nstream_bybnd(ncol, nlay, ngpt, &
974 : op_io%tau, &
975 : op_in%tau, op_in%ssa, &
976 0 : op_io%get_nband(), op_io%get_band_lims_gpoint())
977 : end select
978 :
979 : class is (ty_optical_props_2str)
980 : select type (op_in)
981 : class is (ty_optical_props_1scl)
982 : call inc_2stream_by_1scalar_bybnd(ncol, nlay, ngpt, &
983 : op_io%tau, op_io%ssa, &
984 : op_in%tau, &
985 0 : op_io%get_nband(), op_io%get_band_lims_gpoint())
986 : class is (ty_optical_props_2str)
987 : call inc_2stream_by_2stream_bybnd(ncol, nlay, ngpt, &
988 : op_io%tau, op_io%ssa, op_io%g, &
989 : op_in%tau, op_in%ssa, op_in%g, &
990 389263 : op_io%get_nband(), op_io%get_band_lims_gpoint())
991 : class is (ty_optical_props_nstr)
992 : call inc_2stream_by_nstream_bybnd(ncol, nlay, ngpt, op_in%get_nmom(), &
993 : op_io%tau, op_io%ssa, op_io%g, &
994 : op_in%tau, op_in%ssa, op_in%p, &
995 0 : op_io%get_nband(), op_io%get_band_lims_gpoint())
996 : end select
997 :
998 : class is (ty_optical_props_nstr)
999 : select type (op_in)
1000 : class is (ty_optical_props_1scl)
1001 : call inc_nstream_by_1scalar_bybnd(ncol, nlay, ngpt, &
1002 : op_io%tau, op_io%ssa, &
1003 : op_in%tau, &
1004 0 : op_io%get_nband(), op_io%get_band_lims_gpoint())
1005 : class is (ty_optical_props_2str)
1006 : call inc_nstream_by_2stream_bybnd(ncol, nlay, ngpt, op_io%get_nmom(), &
1007 : op_io%tau, op_io%ssa, op_io%p, &
1008 : op_in%tau, op_in%ssa, op_in%g, &
1009 0 : op_io%get_nband(), op_io%get_band_lims_gpoint())
1010 : class is (ty_optical_props_nstr)
1011 : call inc_nstream_by_nstream_bybnd(ncol, nlay, ngpt, op_io%get_nmom(), op_in%get_nmom(), &
1012 : op_io%tau, op_io%ssa, op_io%p, &
1013 : op_in%tau, op_in%ssa, op_in%p, &
1014 0 : op_io%get_nband(), op_io%get_band_lims_gpoint())
1015 : end select
1016 : end select
1017 : end if
1018 : end function increment
1019 : !> -----------------------------------------------------------------------------------------------
1020 : !>
1021 : !> Routines for array classes: problem sizes
1022 : !>
1023 : !> -----------------------------------------------------------------------------------------------
1024 27184796 : pure function get_arry_extent(this, dim)
1025 : class(ty_optical_props_arry), intent(in ) :: this
1026 : integer, intent(in ) :: dim
1027 : integer :: get_arry_extent
1028 :
1029 27184796 : if(allocated(this%tau)) then
1030 27184796 : get_arry_extent = size(this%tau, dim)
1031 : else
1032 : get_arry_extent = 0
1033 : end if
1034 27184796 : end function get_arry_extent
1035 : ! ------------------------------------------------------------------------------------------
1036 12813872 : pure function get_ncol(this)
1037 : class(ty_optical_props_arry), intent(in ) :: this
1038 : integer :: get_ncol
1039 :
1040 12813872 : get_ncol = get_arry_extent(this, 1)
1041 12813872 : end function get_ncol
1042 : ! ------------------------------------------------------------------------------------------
1043 14370924 : pure function get_nlay(this)
1044 : class(ty_optical_props_arry), intent(in ) :: this
1045 : integer :: get_nlay
1046 :
1047 14370924 : get_nlay = get_arry_extent(this, 2)
1048 14370924 : end function get_nlay
1049 : ! ------------------------------------------------------------------------------------------
1050 0 : pure function get_nmom(this)
1051 : class(ty_optical_props_nstr), intent(in ) :: this
1052 : integer :: get_nmom
1053 :
1054 0 : if(allocated(this%p)) then
1055 0 : get_nmom = size(this%p, 1)
1056 : else
1057 : get_nmom = 0
1058 : end if
1059 0 : end function get_nmom
1060 : ! -----------------------------------------------------------------------------------------------
1061 : !
1062 : ! Routines for base class: spectral discretization
1063 : !
1064 : !> -----------------------------------------------------------------------------------------------
1065 : !>
1066 : !> Number of bands
1067 : !>
1068 27674301 : pure function get_nband(this)
1069 : class(ty_optical_props), intent(in) :: this
1070 : integer :: get_nband
1071 :
1072 27674301 : if(this%is_initialized()) then
1073 27674301 : get_nband = size(this%band2gpt,dim=2)
1074 : else
1075 : get_nband = 0
1076 : end if
1077 27674301 : end function get_nband
1078 : !> -----------------------------------------------------------------------------------------------
1079 : !>
1080 : !> Number of g-points
1081 : !>
1082 28777310 : pure function get_ngpt(this)
1083 : class(ty_optical_props), intent(in) :: this
1084 : integer :: get_ngpt
1085 :
1086 28777310 : if(this%is_initialized()) then
1087 1349353946 : get_ngpt = maxval(this%band2gpt)
1088 : else
1089 : get_ngpt = 0
1090 : end if
1091 28777310 : end function get_ngpt
1092 : !>--------------------------------------------------------------------------------------------------------------------
1093 : !>
1094 : !> The first and last g-point of all bands at once
1095 : !> dimension (2, nbands)
1096 : !>
1097 10607854 : pure function get_band_lims_gpoint(this)
1098 : class(ty_optical_props), intent(in) :: this
1099 : integer, dimension(size(this%band2gpt,dim=1), size(this%band2gpt,dim=2)) &
1100 : :: get_band_lims_gpoint
1101 :
1102 496429066 : get_band_lims_gpoint = this%band2gpt
1103 : end function get_band_lims_gpoint
1104 : !>--------------------------------------------------------------------------------------------------------------------
1105 : !>
1106 : !> First and last g-point of a specific band
1107 : !>
1108 0 : pure function convert_band2gpt(this, band)
1109 : class(ty_optical_props), intent(in) :: this
1110 : integer, intent(in) :: band
1111 : integer, dimension(2) :: convert_band2gpt
1112 :
1113 0 : if(this%is_initialized()) then
1114 0 : convert_band2gpt(:) = this%band2gpt(:,band)
1115 : else
1116 0 : convert_band2gpt(:) = 0
1117 : end if
1118 : end function convert_band2gpt
1119 : !>--------------------------------------------------------------------------------------------------------------------
1120 : !>
1121 : !> Lower and upper wavenumber of all bands
1122 : !> (upper and lower wavenumber by band) = band_lims_wvn(2,band)
1123 : !>
1124 17780193 : pure function get_band_lims_wavenumber(this)
1125 : class(ty_optical_props), intent(in) :: this
1126 : real(wp), dimension(size(this%band_lims_wvn,1), size(this%band_lims_wvn,2)) &
1127 : :: get_band_lims_wavenumber
1128 :
1129 17780193 : if(this%is_initialized()) then
1130 836186571 : get_band_lims_wavenumber(:,:) = this%band_lims_wvn(:,:)
1131 : else
1132 0 : get_band_lims_wavenumber(:,:) = 0._wp
1133 : end if
1134 : end function get_band_lims_wavenumber
1135 : !>--------------------------------------------------------------------------------------------------------------------
1136 : !>
1137 : !> Lower and upper wavelength of all bands
1138 : !>
1139 0 : pure function get_band_lims_wavelength(this)
1140 : class(ty_optical_props), intent(in) :: this
1141 : real(wp), dimension(size(this%band_lims_wvn,1), size(this%band_lims_wvn,2)) &
1142 : :: get_band_lims_wavelength
1143 :
1144 0 : if(this%is_initialized()) then
1145 0 : get_band_lims_wavelength(:,:) = 1._wp/this%band_lims_wvn(:,:)
1146 : else
1147 0 : get_band_lims_wavelength(:,:) = 0._wp
1148 : end if
1149 : end function get_band_lims_wavelength
1150 : !>--------------------------------------------------------------------------------------------------------------------
1151 : !> Bands for all the g-points at once
1152 : !> dimension (ngpt)
1153 : !>
1154 3020006 : pure function get_gpoint_bands(this)
1155 : class(ty_optical_props), intent(in) :: this
1156 : integer, dimension(size(this%gpt2band,dim=1)) &
1157 : :: get_gpoint_bands
1158 :
1159 3020006 : if(this%is_initialized()) then
1160 377099782 : get_gpoint_bands(:) = this%gpt2band(:)
1161 : else
1162 0 : get_gpoint_bands(:) = 0
1163 : end if
1164 : end function get_gpoint_bands
1165 : !>--------------------------------------------------------------------------------------------------------------------
1166 : !>
1167 : !> Band associated with a specific g-point
1168 : !>
1169 17039887227 : pure function convert_gpt2band(this, gpt)
1170 : class(ty_optical_props), intent(in) :: this
1171 : integer, intent(in) :: gpt
1172 : integer :: convert_gpt2band
1173 :
1174 17039887227 : if(this%is_initialized()) then
1175 17039887227 : convert_gpt2band = this%gpt2band(gpt)
1176 : else
1177 : convert_gpt2band = 0
1178 : end if
1179 17039887227 : end function convert_gpt2band
1180 : !>--------------------------------------------------------------------------------------------------------------------
1181 : !>
1182 : !> Expand an array of dimension arr_in(nband) to dimension arr_out(ngpt)
1183 : !>
1184 0 : pure function expand(this, arr_in) result(arr_out)
1185 : class(ty_optical_props), intent(in) :: this
1186 : real(wp), dimension(:), intent(in) :: arr_in ! (nband)
1187 : real(wp), dimension(size(this%gpt2band)) :: arr_out
1188 :
1189 : integer :: iband
1190 :
1191 0 : do iband=1,this%get_nband()
1192 0 : arr_out(this%band2gpt(1,iband):this%band2gpt(2,iband)) = arr_in(iband)
1193 : end do
1194 0 : end function expand
1195 : !>--------------------------------------------------------------------------------------------------------------------
1196 : !>
1197 : !> Are the bands of two objects the same? (same number, same wavelength limits)
1198 : !>
1199 4541596 : pure function bands_are_equal(this, that)
1200 : class(ty_optical_props), intent(in) :: this, that
1201 : logical :: bands_are_equal
1202 :
1203 : bands_are_equal = this%get_nband() == that%get_nband() .and. &
1204 4541596 : this%get_nband() > 0
1205 4541596 : if(.not. bands_are_equal) return
1206 : bands_are_equal = &
1207 : all(abs(this%get_band_lims_wavenumber() - that%get_band_lims_wavenumber()) < &
1208 213195892 : 5._wp * spacing(this%get_band_lims_wavenumber()))
1209 4541596 : end function bands_are_equal
1210 : !>--------------------------------------------------------------------------------------------------------------------
1211 : !>
1212 : !> Is the g-point structure of two objects the same?
1213 : !> (same bands, same number of g-points, same mapping between bands and g-points)
1214 : !>
1215 2270798 : pure function gpoints_are_equal(this, that)
1216 : class(ty_optical_props), intent(in) :: this, that
1217 : logical :: gpoints_are_equal
1218 :
1219 : gpoints_are_equal = this%bands_are_equal(that) .and. &
1220 2270798 : this%get_ngpt() == that%get_ngpt()
1221 2270798 : if(.not. gpoints_are_equal) return
1222 : gpoints_are_equal = &
1223 140238263 : all(this%get_gpoint_bands() == that%get_gpoint_bands())
1224 1135399 : end function gpoints_are_equal
1225 : !> -----------------------------------------------------------------------------------------------
1226 : !>
1227 : !> --- Setting/getting the name
1228 : !>
1229 : !> -----------------------------------------------------------------------------------------------
1230 0 : subroutine set_name(this, name)
1231 : class(ty_optical_props), intent(inout) :: this
1232 : character(len=*), intent(in ) :: name
1233 :
1234 0 : this%name = trim(name)
1235 0 : end subroutine set_name
1236 : ! --------------------------------------------------------
1237 4152333 : function get_name(this)
1238 : class(ty_optical_props), intent(in ) :: this
1239 : character(len=name_len) :: get_name
1240 :
1241 4152333 : get_name = trim(this%name)
1242 4152333 : end function get_name
1243 : ! ------------------------------------------------------------------------------------------
1244 :
1245 2270798 : end module mo_optical_props
|
