LCOV - code coverage report
Current view: | top level - physics/rrtmgp/ext/rte-frontend - mo_optical_props.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 164 | 397 | 41.3 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 26 | 52 | 50.0 %

              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 16161565074 :   pure function is_initialized_base(this)
         316             :     class(ty_optical_props), intent(in) :: this
         317             :     logical                             :: is_initialized_base
         318             :
         319 16161565074 :     is_initialized_base = allocated(this%band2gpt)
         320 16161565074 :   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 16067379449 :   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 16067379449 :     if(this%is_initialized()) then
        1175 16067379449 :       convert_gpt2band = this%gpt2band(gpt)
        1176             :     else
        1177             :       convert_gpt2band = 0
        1178             :     end if
        1179 16067379449 :   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

Generated by: LCOV version 1.14