LCOV - code coverage report
Current view: | top level - physics/rrtmgp/ext/rte-kernels - mo_optical_props_kernels.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 52 | 250 | 20.8 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 5 | 23 | 21.7 %

              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, Trustees of Columbia University.  All right reserved.
           8             : !
           9             : ! Use and duplication is permitted under the terms of the
          10             : !    BSD 3-clause license, see http://opensource.org/licenses/BSD-3-Clause
          11             : ! -------------------------------------------------------------------------------------------------
          12             : !
          13             : !> ## Kernels for arrays of optical properties:
          14             : !>     - delta-scaling
          15             : !>     - adding two sets of properties
          16             : !>     - extracting subsets along the column dimension
          17             : !
          18             : ! -------------------------------------------------------------------------------------------------
          19             :
          20             : module mo_optical_props_kernels
          21             :   use, intrinsic :: iso_c_binding
          22             :   use mo_rte_kind, only: wp, wl
          23             :   implicit none
          24             :
          25             :   public
          26             :
          27             :   !> Delta-scale two-stream optical properties
          28             :   interface delta_scale_2str_kernel
          29             :     module procedure delta_scale_2str_f_k, delta_scale_2str_k
          30             :   end interface
          31             :
          32             :   !> Subsetting, meaning extracting some portion of the 3D domain
          33             :   interface extract_subset
          34             :     module procedure extract_subset_dim1_3d, extract_subset_dim2_4d
          35             :     module procedure extract_subset_absorption_tau
          36             :   end interface extract_subset
          37             :
          38             :   real(wp), parameter, private :: eps = 3.0_wp*tiny(1.0_wp)
          39             : contains
          40             :   ! -------------------------------------------------------------------------------------------------
          41             :   !
          42             :   ! Delta-scaling is provided only for two-stream properties at present
          43             :   !
          44             :   ! -------------------------------------------------------------------------------------------------
          45             :   !> Delta-scale two-stream optical properties given user-provided value of \(f\) (forward scattering)
          46             :   !
          47           0 :   pure subroutine delta_scale_2str_f_k(ncol, nlay, ngpt, tau, ssa, g, f) &
          48             :       bind(C, name="rte_delta_scale_2str_f_k")
          49             :     integer,                               intent(in   ) :: ncol, nlay, ngpt
          50             :       !! Array sizes
          51             :     real(wp), dimension(ncol, nlay, ngpt), intent(inout) ::  tau, ssa, g
          52             :       !! Optical depth, single-scattering albedo, asymmetry parameter
          53             :     real(wp), dimension(ncol, nlay, ngpt), intent(in   ) ::  f
          54             :       !! User-provided forward-scattering fraction
          55             :
          56             :     real(wp) :: wf
          57             :     integer  :: icol, ilay, igpt
          58             :
          59           0 :     do igpt = 1, ngpt
          60           0 :       do ilay = 1, nlay
          61           0 :         do icol = 1, ncol
          62           0 :           wf = ssa(icol,ilay,igpt) * f(icol,ilay,igpt)
          63           0 :           tau(icol,ilay,igpt) = (1._wp - wf) * tau(icol,ilay,igpt)
          64           0 :           ssa(icol,ilay,igpt) = (ssa(icol,ilay,igpt) - wf) /  max(eps,(1.0_wp - wf))
          65             :           g  (icol,ilay,igpt) = (g  (icol,ilay,igpt) - f(icol,ilay,igpt)) / &
          66           0 :                                         max(eps,(1._wp - f(icol,ilay,igpt)))
          67             :         end do
          68             :       end do
          69             :     end do
          70             :
          71           0 :   end subroutine delta_scale_2str_f_k
          72             :   ! ---------------------------------
          73             :   !> Delta-scale assuming forward-scatternig fraction is the square of the asymmetry parameter
          74             :   !>    i.e. \(f = g^2\)
          75             :   !
          76      389263 :   pure subroutine delta_scale_2str_k(ncol, nlay, ngpt, tau, ssa, g) &
          77             :       bind(C, name="rte_delta_scale_2str_k")
          78             :     integer,                               intent(in   ) :: ncol, nlay, ngpt
          79             :       !! Array sizes
          80             :     real(wp), dimension(ncol, nlay, ngpt), intent(inout) ::  tau, ssa, g
          81             :       !! Optical depth, single-scattering albedo, asymmetry parameter
          82             :
          83             :     real(wp) :: f, wf
          84             :     integer  :: icol, ilay, igpt
          85             :
          86    43986719 :     do igpt = 1, ngpt
          87  4142147583 :       do ilay = 1, nlay
          88 65796884720 :         do icol = 1, ncol
          89 61655126400 :           f  = g  (icol,ilay,igpt) * g  (icol,ilay,igpt)
          90 61655126400 :           wf = ssa(icol,ilay,igpt) * f
          91 61655126400 :           tau(icol,ilay,igpt) = (1._wp - wf) * tau(icol,ilay,igpt)
          92 61655126400 :           ssa(icol,ilay,igpt) = (ssa(icol,ilay,igpt) - wf) /  max(eps,(1.0_wp - wf))
          93 65753287264 :           g  (icol,ilay,igpt) = (g  (icol,ilay,igpt) -  f) /  max(eps,(1.0_wp -  f))
          94             :         end do
          95             :       end do
          96             :     end do
          97             :
          98      389263 :   end subroutine delta_scale_2str_k
          99             :   ! -------------------------------------------------------------------------------------------------
         100             :   !
         101             :   ! Addition of optical properties: the first set are incremented by the second set.
         102             :   !
         103             :   !   There are three possible representations of optical properties (scalar = optical depth only;
         104             :   !   two-stream = tau, single-scattering albedo, and asymmetry factor g, and
         105             :   !   n-stream = tau, ssa, and phase function moments p.) Thus we need nine routines, three for
         106             :   !   each choice of representation on the left hand side times three representations of the
         107             :   !   optical properties to be added.
         108             :   !
         109             :   !   There are two sets of these nine routines. In the first the two sets of optical
         110             :   !   properties are defined at the same spectral resolution. There is also a set of routines
         111             :   !   to add properties defined at lower spectral resolution to a set defined at higher spectral
         112             :   !   resolution (adding properties defined by band to those defined by g-point)
         113             :   !
         114             :   ! -------------------------------------------------------------------------------------------------
         115             :   !> increase one absorption optical depth by a second value
         116      746136 :   pure subroutine increment_1scalar_by_1scalar(ncol, nlay, ngpt, &
         117      746136 :                                                tau1,             &
         118      746136 :                                                tau2) bind(C, name="rte_increment_1scalar_by_1scalar")
         119             :     integer,                             intent(in   ) :: ncol, nlay, ngpt !! array sizes
         120             :     real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1             !! optical properties to be modified
         121             :     real(wp), dimension(ncol,nlay,ngpt), intent(in   ) :: tau2             !! optical properties to be added to original
         122             :
         123             :     integer  :: icol, ilay, igpt
         124             :
         125    96251544 :     do igpt = 1, ngpt
         126  9073759896 :       do ilay = 1, nlay
         127 >14999*10^7 :         do icol = 1, ncol
         128 >14990*10^7 :           tau1(icol,ilay,igpt) = tau1(icol,ilay,igpt) + tau2(icol,ilay,igpt)
         129             :         end do
         130             :       end do
         131             :     end do
         132      746136 :   end subroutine increment_1scalar_by_1scalar
         133             :   ! ---------------------------------
         134             :   !> increase absorption optical depth with extinction optical depth (2-stream form)
         135           0 :   pure subroutine increment_1scalar_by_2stream(ncol, nlay, ngpt, &
         136           0 :                                                tau1,             &
         137           0 :                                                tau2, ssa2) bind(C, name="rte_increment_1scalar_by_2stream")
         138             :     integer,                             intent(in   ) :: ncol, nlay, ngpt !! array sizes
         139             :     real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1             !! optical properties to be modified
         140             :     real(wp), dimension(ncol,nlay,ngpt), intent(in   ) :: tau2, ssa2       !! optical properties to be added to original
         141             :
         142             :     integer  :: icol, ilay, igpt
         143             :
         144           0 :     do igpt = 1, ngpt
         145           0 :       do ilay = 1, nlay
         146           0 :         do icol = 1, ncol
         147           0 :           tau1(icol,ilay,igpt) = tau1(icol,ilay,igpt) + &
         148           0 :                                  tau2(icol,ilay,igpt) * (1._wp - ssa2(icol,ilay,igpt))
         149             :         end do
         150             :       end do
         151             :     end do
         152           0 :   end subroutine increment_1scalar_by_2stream
         153             :   ! ---------------------------------
         154             :   !> increase absorption optical depth with extinction optical depth (n-stream form)
         155           0 :   pure subroutine increment_1scalar_by_nstream(ncol, nlay, ngpt, &
         156           0 :                                                tau1,             &
         157           0 :                                                tau2, ssa2) bind(C, name="rte_increment_1scalar_by_nstream")
         158             :     integer,                             intent(in   ) :: ncol, nlay, ngpt !! array sizes
         159             :     real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1             !! optical properties to be modified
         160             :     real(wp), dimension(ncol,nlay,ngpt), intent(in   ) :: tau2, ssa2       !! optical properties to be added to original
         161             :
         162             :     integer  :: icol, ilay, igpt
         163             :
         164           0 :     do igpt = 1, ngpt
         165           0 :       do ilay = 1, nlay
         166           0 :         do icol = 1, ncol
         167           0 :           tau1(icol,ilay,igpt) = tau1(icol,ilay,igpt) + &
         168           0 :                                  tau2(icol,ilay,igpt) * (1._wp - ssa2(icol,ilay,igpt))
         169             :         end do
         170             :       end do
         171             :     end do
         172           0 :   end subroutine increment_1scalar_by_nstream
         173             :   ! ---------------------------------
         174             :   ! ---------------------------------
         175             :   !> increment two-stream optical properties \(\tau, \omega_0, g\) with absorption optical depth
         176           0 :   pure subroutine increment_2stream_by_1scalar(ncol, nlay, ngpt, &
         177           0 :                                                tau1, ssa1,       &
         178           0 :                                                tau2) bind(C, name="rte_increment_2stream_by_1scalar")
         179             :     integer,                             intent(in   ) :: ncol, nlay, ngpt !! array sizes
         180             :     real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1       !! optical properties to be modified
         181             :     real(wp), dimension(ncol,nlay,ngpt), intent(in   ) :: tau2             !! optical properties to be added to original
         182             :
         183             :     integer  :: icol, ilay, igpt
         184             :     real(wp) :: tau12
         185             :
         186           0 :     do igpt = 1, ngpt
         187           0 :       do ilay = 1, nlay
         188           0 :         do icol = 1, ncol
         189           0 :           tau12 = tau1(icol,ilay,igpt) + tau2(icol,ilay,igpt)
         190           0 :           ssa1(icol,ilay,igpt) = tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) / max(eps,tau12)
         191           0 :           tau1(icol,ilay,igpt) = tau12
         192             :           ! g is unchanged
         193             :         end do
         194             :       end do
         195             :     end do
         196           0 :   end subroutine increment_2stream_by_1scalar
         197             :   ! ---------------------------------
         198             :   !> increment two-stream optical properties \(\tau, \omega_0, g\) with a second set
         199      389263 :   pure subroutine increment_2stream_by_2stream(ncol, nlay, ngpt, &
         200      389263 :                                                tau1, ssa1, g1,   &
         201      389263 :                                                tau2, ssa2, g2) bind(C, name="rte_increment_2stream_by_2stream")
         202             :     integer,                             intent(in   ) :: ncol, nlay, ngpt !! array sizes
         203             :     real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1, g1   !! optical properties to be modified
         204             :     real(wp), dimension(ncol,nlay,ngpt), intent(in   ) :: tau2, ssa2, g2   !! optical properties to be added to original
         205             :
         206             :     integer :: icol, ilay, igpt
         207             :     real(wp) :: tau12, tauscat12
         208             :
         209    43986719 :     do igpt = 1, ngpt
         210  4142147583 :       do ilay = 1, nlay
         211 65796884720 :         do icol = 1, ncol
         212             :           ! t=tau1 + tau2
         213 61655126400 :           tau12 = tau1(icol,ilay,igpt) + tau2(icol,ilay,igpt)
         214             :           ! w=(tau1*ssa1 + tau2*ssa2) / t
         215             :           tauscat12 = tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) + &
         216 61655126400 :                       tau2(icol,ilay,igpt) * ssa2(icol,ilay,igpt)
         217             :           g1(icol,ilay,igpt) = &
         218             :             (tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) * g1(icol,ilay,igpt) + &
         219             :              tau2(icol,ilay,igpt) * ssa2(icol,ilay,igpt) * g2(icol,ilay,igpt)) &
         220 61655126400 :               / max(eps,tauscat12)
         221 61655126400 :           ssa1(icol,ilay,igpt) = tauscat12 / max(eps,tau12)
         222 65753287264 :           tau1(icol,ilay,igpt) = tau12
         223             :         end do
         224             :       end do
         225             :     end do
         226      389263 :   end subroutine increment_2stream_by_2stream
         227             :   ! ---------------------------------
         228             :   !> increment two-stream optical properties \(\tau, \omega_0, g\) with _n_-stream
         229           0 :   pure subroutine increment_2stream_by_nstream(ncol, nlay, ngpt, nmom2, &
         230           0 :                                                tau1, ssa1, g1,          &
         231           0 :                                                tau2, ssa2, p2) bind(C, name="rte_increment_2stream_by_nstream")
         232             :     integer,                             intent(in   ) :: ncol, nlay, ngpt, nmom2  !! array sizes
         233             :     real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1, g1           !! optical properties to be modified
         234             :     real(wp), dimension(ncol,nlay,ngpt), intent(in   ) :: tau2, ssa2               !! optical properties to be added to original
         235             :     real(wp), dimension(nmom2, &
         236             :                         ncol,nlay,ngpt), intent(in   ) :: p2                       !! moments of the phase function to be added
         237             :
         238             :     integer  :: icol, ilay, igpt
         239             :     real(wp) :: tau12, tauscat12
         240             :
         241           0 :     do igpt = 1, ngpt
         242           0 :       do ilay = 1, nlay
         243           0 :         do icol = 1, ncol
         244             :           ! t=tau1 + tau2
         245           0 :           tau12 = tau1(icol,ilay,igpt) + tau2(icol,ilay,igpt)
         246             :           ! w=(tau1*ssa1 + tau2*ssa2) / t
         247             :           tauscat12 = &
         248             :              tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) + &
         249           0 :              tau2(icol,ilay,igpt) * ssa2(icol,ilay,igpt)
         250             :           g1(icol,ilay,igpt) = &
         251             :             (tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) * g1(   icol,ilay,igpt)+ &
         252           0 :              tau2(icol,ilay,igpt) * ssa2(icol,ilay,igpt) * p2(1, icol,ilay,igpt)) / max(eps,tauscat12)
         253           0 :           ssa1(icol,ilay,igpt) = tauscat12 / max(eps,tau12)
         254           0 :           tau1(icol,ilay,igpt) = tau12
         255             :         end do
         256             :       end do
         257             :     end do
         258           0 :   end subroutine increment_2stream_by_nstream
         259             :   ! ---------------------------------
         260             :   ! ---------------------------------
         261             :   !> increment _n_-stream optical properties \(\tau, \omega_0, p\) with absorption optical depth
         262           0 :   pure subroutine increment_nstream_by_1scalar(ncol, nlay, ngpt, &
         263           0 :                                                tau1, ssa1,       &
         264           0 :                                                tau2) bind(C, name="rte_increment_nstream_by_1scalar")
         265             :     integer,                             intent(in   ) :: ncol, nlay, ngpt  !! array sizes
         266             :     real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1        !! optical properties to be modified
         267             :     real(wp), dimension(ncol,nlay,ngpt), intent(in   ) :: tau2              !! optical properties to be added to original
         268             :
         269             :     integer  :: icol, ilay, igpt
         270             :     real(wp) :: tau12
         271             :
         272           0 :     do igpt = 1, ngpt
         273           0 :       do ilay = 1, nlay
         274           0 :         do icol = 1, ncol
         275           0 :           tau12 = tau1(icol,ilay,igpt) + tau2(icol,ilay,igpt)
         276           0 :           ssa1(icol,ilay,igpt) = tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) / max(eps,tau12)
         277           0 :           tau1(icol,ilay,igpt) = tau12
         278             :           ! p is unchanged
         279             :         end do
         280             :       end do
         281             :     end do
         282           0 :   end subroutine increment_nstream_by_1scalar
         283             :   ! ---------------------------------
         284             :   !> increment _n_-stream optical properties \(\tau, \omega_0, p\) with two-stream values
         285           0 :   pure subroutine increment_nstream_by_2stream(ncol, nlay, ngpt, nmom1, &
         286           0 :                                                tau1, ssa1, p1,          &
         287           0 :                                                tau2, ssa2, g2) bind(C, name="rte_increment_nstream_by_2stream")
         288             :     integer,                              intent(in   ) :: ncol, nlay, ngpt, nmom1  !! array sizes
         289             :     real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1                !! optical properties to be modified
         290             :     real(wp), dimension(nmom1, &
         291             :                         ncol,nlay,ngpt), intent(inout) :: p1                        !! moments of the phase function be modified
         292             :     real(wp), dimension(ncol,nlay,ngpt), intent(in   ) :: tau2, ssa2, g2            !! optical properties to be added to original
         293             :
         294             :     integer  :: icol, ilay, igpt
         295             :     real(wp) :: tau12, tauscat12
         296           0 :     real(wp), dimension(nmom1) :: temp_moms ! TK
         297             :     integer  :: imom  !TK
         298             :
         299           0 :     do igpt = 1, ngpt
         300           0 :       do ilay = 1, nlay
         301           0 :         do icol = 1, ncol
         302           0 :           tau12 = tau1(icol,ilay,igpt) + tau2(icol,ilay,igpt)
         303             :           tauscat12 = &
         304             :              tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) + &
         305           0 :              tau2(icol,ilay,igpt) * ssa2(icol,ilay,igpt)
         306             :           !
         307             :           ! Here assume Henyey-Greenstein
         308             :           !
         309           0 :           temp_moms(1) = g2(icol,ilay,igpt)
         310           0 :           do imom = 2, nmom1
         311           0 :             temp_moms(imom) = temp_moms(imom-1) * g2(icol,ilay,igpt)
         312             :           end do
         313             :           p1(1:nmom1, icol,ilay,igpt) = &
         314             :               (tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) * p1(1:nmom1, icol,ilay,igpt) + &
         315           0 :                tau2(icol,ilay,igpt) * ssa2(icol,ilay,igpt) * temp_moms(1:nmom1)  ) / max(eps,tauscat12)
         316           0 :           ssa1(icol,ilay,igpt) = tauscat12 / max(eps,tau12)
         317           0 :           tau1(icol,ilay,igpt) = tau12
         318             :         end do
         319             :       end do
         320             :     end do
         321           0 :   end subroutine increment_nstream_by_2stream
         322             :   ! ---------------------------------
         323             :   !> increment one set of _n_-stream optical properties with another set
         324           0 :   pure subroutine increment_nstream_by_nstream(ncol, nlay, ngpt, nmom1, nmom2, &
         325           0 :                                                tau1, ssa1, p1,                 &
         326           0 :                                                tau2, ssa2, p2) bind(C, name="rte_increment_nstream_by_nstream")
         327             :     integer,                              intent(in   ) :: ncol, nlay, ngpt, nmom1, nmom2  !! array sizes
         328             :     real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1   !! optical properties to be modified
         329             :     real(wp), dimension(nmom1, &
         330             :                         ncol,nlay,ngpt), intent(inout) :: p1           !! moments of the phase function be modified
         331             :     real(wp), dimension(ncol,nlay,ngpt), intent(in   ) :: tau2, ssa2   !! optical properties to be added to original
         332             :     real(wp), dimension(nmom2, &
         333             :                         ncol,nlay,ngpt), intent(in   ) :: p2           !! moments of the phase function to be added
         334             :
         335             :     integer  :: icol, ilay, igpt, mom_lim
         336             :     real(wp) :: tau12, tauscat12
         337             :
         338           0 :     mom_lim = min(nmom1, nmom2)
         339           0 :     do igpt = 1, ngpt
         340           0 :       do ilay = 1, nlay
         341           0 :         do icol = 1, ncol
         342           0 :           tau12 = tau1(icol,ilay,igpt) + tau2(icol,ilay,igpt)
         343             :           tauscat12 = &
         344             :              tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) + &
         345           0 :              tau2(icol,ilay,igpt) * ssa2(icol,ilay,igpt)
         346             :           !
         347             :           ! If op2 has more moments than op1 these are ignored;
         348             :           !   if it has fewer moments the higher orders are assumed to be 0
         349             :           !
         350             :           p1(1:mom_lim, icol,ilay,igpt) = &
         351           0 :               (tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) * p1(1:mom_lim, icol,ilay,igpt) + &
         352           0 :                tau2(icol,ilay,igpt) * ssa2(icol,ilay,igpt) * p2(1:mom_lim, icol,ilay,igpt)) / max(eps,tauscat12)
         353           0 :           ssa1(icol,ilay,igpt) = tauscat12 / max(eps,tau12)
         354           0 :           tau1(icol,ilay,igpt) = tau12
         355             :         end do
         356             :       end do
         357             :     end do
         358           0 :   end subroutine increment_nstream_by_nstream
         359             :   ! -------------------------------------------------------------------------------------------------
         360             :   !
         361             :   ! Incrementing when the second set of optical properties is defined at lower spectral resolution
         362             :   !   (e.g. by band instead of by gpoint)
         363             :   !
         364             :   ! -------------------------------------------------------------------------------------------------
         365             :   !> increase one absorption optical depth defined on g-points by a second value defined on bands
         366      746136 :   pure subroutine inc_1scalar_by_1scalar_bybnd(ncol, nlay, ngpt, &
         367      746136 :                                                tau1,             &
         368      746136 :                                                tau2,             &
         369      746136 :                                                nbnd, gpt_lims) bind(C, name="rte_inc_1scalar_by_1scalar_bybnd")
         370             :     integer,                             intent(in   ) :: ncol, nlay, ngpt, nbnd  !! array sizes
         371             :     real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1     !! optical properties to be modified (defined on g-points)
         372             :     real(wp), dimension(ncol,nlay,nbnd), intent(in   ) :: tau2     !! optical properties to be added to original (defined on bands)
         373             :     integer,  dimension(2,nbnd),         intent(in   ) :: gpt_lims !! Starting and ending gpoint for each band
         374             :
         375             :     integer :: ibnd, igpt
         376             :
         377    12684312 :     do ibnd = 1, nbnd
         378   108189720 :       do igpt = gpt_lims(1, ibnd), gpt_lims(2, ibnd)
         379 >15001*10^7 :         tau1(:,:,igpt) = tau1(:,:,igpt) + tau2(:,:,ibnd)
         380             :       end do
         381             :     end do
         382      746136 :   end subroutine inc_1scalar_by_1scalar_bybnd
         383             :   ! ---------------------------------
         384             :   !> increase absorption optical depth defined on g-points  with extinction optical depth (2-stream form) defined on bands
         385           0 :   pure subroutine inc_1scalar_by_2stream_bybnd(ncol, nlay, ngpt, &
         386           0 :                                                tau1,             &
         387           0 :                                                tau2, ssa2,       &
         388           0 :                                                nbnd, gpt_lims) bind(C, name="rte_inc_1scalar_by_2stream_bybnd")
         389             :     integer,                             intent(in   ) :: ncol, nlay, ngpt, nbnd  !! array sizes
         390             :     real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1        !! optical properties to be modified (defined on g-points)
         391             :     real(wp), dimension(ncol,nlay,nbnd), intent(in   ) :: tau2, ssa2  !! optical properties to be added to original (defined on bands)
         392             :     integer,  dimension(2,nbnd),         intent(in   ) :: gpt_lims    !! Starting and ending gpoint for each band
         393             :
         394             :     integer :: ibnd, igpt
         395             :
         396           0 :     do ibnd = 1, nbnd
         397           0 :       do igpt = gpt_lims(1, ibnd), gpt_lims(2, ibnd)
         398           0 :         tau1(:,:,igpt) = tau1(:,:,igpt) + tau2(:,:,ibnd) * (1._wp - ssa2(:,:,ibnd))
         399             :       end do
         400             :     end do
         401           0 :   end subroutine inc_1scalar_by_2stream_bybnd
         402             :   ! ---------------------------------
         403             :   !> increase absorption optical depth defined on g-points  with extinction optical depth (n-stream form) defined on bands
         404           0 :   pure subroutine inc_1scalar_by_nstream_bybnd(ncol, nlay, ngpt, &
         405           0 :                                                tau1,             &
         406           0 :                                                tau2, ssa2,       &
         407           0 :                                                nbnd, gpt_lims) bind(C, name="rte_inc_1scalar_by_nstream_bybnd")
         408             :     integer,                             intent(in   ) :: ncol, nlay, ngpt, nbnd  !! array sizes
         409             :     real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1       !! optical properties to be modified (defined on g-points)
         410             :     real(wp), dimension(ncol,nlay,nbnd), intent(in   ) :: tau2, ssa2 !! optical properties to be added to original (defined on bands)
         411             :     integer,  dimension(2,nbnd),         intent(in   ) :: gpt_lims   !! Starting and ending gpoint for each band
         412             :
         413             :     integer :: ibnd, igpt
         414             :
         415           0 :     do ibnd = 1, nbnd
         416           0 :       do igpt = gpt_lims(1, ibnd), gpt_lims(2, ibnd)
         417           0 :         tau1(:,:,igpt) = tau1(:,:,igpt) + tau2(:,:,ibnd) * (1._wp - ssa2(:,:,ibnd))
         418             :       end do
         419             :     end do
         420           0 :   end subroutine inc_1scalar_by_nstream_bybnd
         421             :
         422             :     ! ---------------------------------
         423             :   !> increment two-stream optical properties \(\tau, \omega_0, g\) defined on g-points with absorption optical depth defined on bands
         424           0 :   pure subroutine inc_2stream_by_1scalar_bybnd(ncol, nlay, ngpt, &
         425           0 :                                                tau1, ssa1,       &
         426           0 :                                                tau2,             &
         427           0 :                                                nbnd, gpt_lims) bind(C, name="rte_inc_2stream_by_1scalar_bybnd")
         428             :     integer,                             intent(in   ) :: ncol, nlay, ngpt, nbnd  !! array sizes
         429             :     real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1 !! optical properties to be modified (defined on g-points)
         430             :     real(wp), dimension(ncol,nlay,nbnd), intent(in   ) :: tau2       !! optical properties to be added to original (defined on bands)
         431             :     integer,  dimension(2,nbnd),         intent(in   ) :: gpt_lims   !! Starting and ending gpoint for each band
         432             :
         433             :     integer  :: icol, ilay, igpt, ibnd
         434             :     real(wp) :: tau12
         435             :
         436           0 :     do ibnd = 1, nbnd
         437           0 :       do igpt = gpt_lims(1, ibnd), gpt_lims(2, ibnd)
         438           0 :         do ilay = 1, nlay
         439           0 :           do icol = 1, ncol
         440           0 :             tau12 = tau1(icol,ilay,igpt) + tau2(icol,ilay,ibnd)
         441           0 :             ssa1(icol,ilay,igpt) = tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) / max(eps,tau12)
         442           0 :             tau1(icol,ilay,igpt) = tau12
         443             :             ! g is unchanged
         444             :           end do
         445             :         end do
         446             :       end do
         447             :     end do
         448           0 :   end subroutine inc_2stream_by_1scalar_bybnd
         449             :   ! ---------------------------------
         450             :   !> increment 2-stream optical properties defined on g-points with another set defined on bands
         451      389263 :   pure subroutine inc_2stream_by_2stream_bybnd(ncol, nlay, ngpt, &
         452      389263 :                                                tau1, ssa1, g1,   &
         453      389263 :                                                tau2, ssa2, g2,   &
         454      389263 :                                                nbnd, gpt_lims) bind(C, name="rte_inc_2stream_by_2stream_bybnd")
         455             :     integer,                             intent(in   ) :: ncol, nlay, ngpt, nbnd  !! array sizes
         456             :     real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1, g1 !! optical properties to be modified (defined on g-points)
         457             :     real(wp), dimension(ncol,nlay,nbnd), intent(in   ) :: tau2, ssa2, g2 !! optical properties to be added to original (defined on bands)
         458             :     integer,  dimension(2,nbnd),         intent(in   ) :: gpt_lims       !! Starting and ending gpoint for each band
         459             :
         460             :     integer  :: icol, ilay, igpt, ibnd
         461             :     real(wp) :: tau12, tauscat12
         462             :
         463     5838945 :     do ibnd = 1, nbnd
         464    49436401 :       do igpt = gpt_lims(1, ibnd), gpt_lims(2, ibnd)
         465  4147208002 :         do ilay = 1, nlay
         466 65796884720 :           do icol = 1, ncol
         467             :             ! t=tau1 + tau2
         468 61655126400 :             tau12 = tau1(icol,ilay,igpt) + tau2(icol,ilay,ibnd)
         469             :             ! w=(tau1*ssa1 + tau2*ssa2) / t
         470             :             tauscat12 = &
         471             :                tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) + &
         472 61655126400 :                tau2(icol,ilay,ibnd) * ssa2(icol,ilay,ibnd)
         473             :             g1(icol,ilay,igpt) = &
         474             :               (tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) * g1(icol,ilay,igpt) + &
         475 61655126400 :                tau2(icol,ilay,ibnd) * ssa2(icol,ilay,ibnd) * g2(icol,ilay,ibnd)) / max(eps,tauscat12)
         476 61655126400 :             ssa1(icol,ilay,igpt) = tauscat12 / max(eps,tau12)
         477 65753287264 :             tau1(icol,ilay,igpt) = tau12
         478             :           end do
         479             :         end do
         480             :       end do
         481             :     end do
         482      389263 :   end subroutine inc_2stream_by_2stream_bybnd
         483             :   ! ---------------------------------
         484             :   !> increment 2-stream optical properties defined on g-points with _n_-stream properties set defined on bands
         485           0 :   pure subroutine inc_2stream_by_nstream_bybnd(ncol, nlay, ngpt, nmom2, &
         486           0 :                                                tau1, ssa1, g1,          &
         487           0 :                                                tau2, ssa2, p2,          &
         488           0 :                                                nbnd, gpt_lims) bind(C, name="rte_inc_2stream_by_nstream_bybnd")
         489             :     integer,                             intent(in   ) :: ncol, nlay, ngpt, nmom2, nbnd  !! array sizes
         490             :     real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1, g1 !! optical properties to be modified (defined on g-points)
         491             :     real(wp), dimension(ncol,nlay,nbnd), intent(in   ) :: tau2, ssa2     !! optical properties to be added to original (defined on bands)
         492             :     real(wp), dimension(nmom2, &
         493             :                         ncol,nlay,nbnd), intent(in   ) :: p2             !! moments of the phase function to be added
         494             :     integer,  dimension(2,nbnd),         intent(in   ) :: gpt_lims       !! Starting and ending gpoint for each band
         495             :
         496             :     integer  :: icol, ilay, igpt, ibnd
         497             :     real(wp) :: tau12, tauscat12
         498             :
         499           0 :     do ibnd = 1, nbnd
         500           0 :       do igpt = gpt_lims(1, ibnd), gpt_lims(2, ibnd)
         501           0 :         do ilay = 1, nlay
         502           0 :           do icol = 1, ncol
         503             :             ! t=tau1 + tau2
         504           0 :             tau12 = tau1(icol,ilay,igpt) + tau2(icol,ilay,ibnd)
         505             :             ! w=(tau1*ssa1 + tau2*ssa2) / t
         506             :             tauscat12 = &
         507             :                tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) + &
         508           0 :                tau2(icol,ilay,ibnd) * ssa2(icol,ilay,ibnd)
         509             :             g1(icol,ilay,igpt) = &
         510             :               (tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) * g1(   icol,ilay,igpt)+ &
         511           0 :                tau2(icol,ilay,ibnd) * ssa2(icol,ilay,ibnd) * p2(1, icol,ilay,ibnd)) / max(eps,tauscat12)
         512           0 :             ssa1(icol,ilay,igpt) = tauscat12 / max(eps,tau12)
         513           0 :             tau1(icol,ilay,igpt) = tau12
         514             :           end do
         515             :         end do
         516             :       end do
         517             :     end do
         518           0 :   end subroutine inc_2stream_by_nstream_bybnd
         519             :   ! ---------------------------------
         520             :   ! ---------------------------------
         521             :   !> increment _n_-stream optical properties defined on g-points with absorption optical depth defined on bands
         522           0 :   pure subroutine inc_nstream_by_1scalar_bybnd(ncol, nlay, ngpt, &
         523           0 :                                                tau1, ssa1,       &
         524           0 :                                                tau2,             &
         525           0 :                                                nbnd, gpt_lims) bind(C, name="rte_inc_nstream_by_1scalar_bybnd")
         526             :     integer,                             intent(in   ) :: ncol, nlay, ngpt, nbnd  !! array sizes
         527             :     real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1 !! optical properties to be modified (defined on g-points)
         528             :     real(wp), dimension(ncol,nlay,nbnd), intent(in   ) :: tau2       !! optical properties to be added to original (defined on bands)
         529             :     integer,  dimension(2,nbnd),         intent(in   ) :: gpt_lims   !! Starting and ending gpoint for each band
         530             :
         531             :     integer  :: icol, ilay, igpt, ibnd
         532             :     real(wp) :: tau12
         533             :
         534           0 :     do ibnd = 1, nbnd
         535           0 :       do igpt = gpt_lims(1, ibnd), gpt_lims(2, ibnd)
         536           0 :         do ilay = 1, nlay
         537           0 :           do icol = 1, ncol
         538           0 :             tau12 = tau1(icol,ilay,igpt) + tau2(icol,ilay,ibnd)
         539           0 :             ssa1(icol,ilay,igpt) = tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) / max(eps,tau12)
         540           0 :             tau1(icol,ilay,igpt) = tau12
         541             :             ! p is unchanged
         542             :           end do
         543             :         end do
         544             :       end do
         545             :     end do
         546           0 :   end subroutine inc_nstream_by_1scalar_bybnd
         547             :   ! ---------------------------------
         548             :   !> increment n-stream optical properties defined on g-points with 2-stream properties set defined on bands
         549           0 :   pure subroutine inc_nstream_by_2stream_bybnd(ncol, nlay, ngpt, nmom1, &
         550           0 :                                                tau1, ssa1, p1,          &
         551           0 :                                                tau2, ssa2, g2,          &
         552           0 :                                                nbnd, gpt_lims) bind(C, name="rte_inc_nstream_by_2stream_bybnd")
         553             :     integer,                             intent(in   ) :: ncol, nlay, ngpt, nmom1, nbnd  !! array sizes
         554             :     real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1     !! optical properties to be modified (defined on g-points)
         555             :     real(wp), dimension(nmom1, &
         556             :                         ncol,nlay,ngpt), intent(inout) :: p1             !! moments of the phase function be modified
         557             :     real(wp), dimension(ncol,nlay,nbnd), intent(in   ) :: tau2, ssa2, g2 !! optical properties to be added to original (defined on bands)
         558             :     integer,  dimension(2,nbnd),         intent(in   ) :: gpt_lims       !! Starting and ending gpoint for each band
         559             :
         560             :     integer  :: icol, ilay, igpt, ibnd
         561             :     real(wp) :: tau12, tauscat12
         562           0 :     real(wp), dimension(nmom1) :: temp_moms ! TK
         563             :     integer  :: imom  !TK
         564             :
         565           0 :     do ibnd = 1, nbnd
         566           0 :       do igpt = gpt_lims(1, ibnd), gpt_lims(2, ibnd)
         567           0 :         do ilay = 1, nlay
         568           0 :           do icol = 1, ncol
         569           0 :             tau12 = tau1(icol,ilay,igpt) + tau2(icol,ilay,ibnd)
         570             :             tauscat12 = &
         571             :                tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) + &
         572           0 :                tau2(icol,ilay,ibnd) * ssa2(icol,ilay,ibnd)
         573             :             !
         574             :             ! Here assume Henyey-Greenstein
         575             :             !
         576           0 :             temp_moms(1) = g2(icol,ilay,ibnd)
         577           0 :             do imom = 2, nmom1
         578           0 :               temp_moms(imom) = temp_moms(imom-1) * g2(icol,ilay,ibnd)
         579             :             end do
         580             :             p1(1:nmom1, icol,ilay,igpt) = &
         581             :                 (tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) * p1(1:nmom1, icol,ilay,igpt) + &
         582           0 :                  tau2(icol,ilay,ibnd) * ssa2(icol,ilay,ibnd) * temp_moms(1:nmom1)  ) / max(eps,tauscat12)
         583           0 :             ssa1(icol,ilay,igpt) = tauscat12 / max(eps,tau12)
         584           0 :             tau1(icol,ilay,igpt) = tau12
         585             :           end do
         586             :         end do
         587             :       end do
         588             :     end do
         589           0 :   end subroutine inc_nstream_by_2stream_bybnd
         590             :   ! ---------------------------------
         591             :   !> increment _n_-stream optical properties defined on g-points with a second set defined on bands
         592           0 :   pure subroutine inc_nstream_by_nstream_bybnd(ncol, nlay, ngpt, nmom1, nmom2, &
         593           0 :                                                tau1, ssa1, p1,                 &
         594           0 :                                                tau2, ssa2, p2,                 &
         595           0 :                                                nbnd, gpt_lims) bind(C, name="rte_inc_nstream_by_nstream_bybnd")
         596             :     integer,                             intent(in   ) :: ncol, nlay, ngpt, nmom1, nmom2, nbnd  !! array sizes
         597             :     real(wp), dimension(ncol,nlay,ngpt), intent(inout) :: tau1, ssa1 !! optical properties to be modified (defined on g-points)
         598             :     real(wp), dimension(nmom1, &
         599             :                         ncol,nlay,ngpt), intent(inout) :: p1         !! moments of the phase function be modified
         600             :     real(wp), dimension(ncol,nlay,nbnd), intent(in   ) :: tau2, ssa2 !! optical properties to be added to original (defined on bands)
         601             :     real(wp), dimension(nmom2, &
         602             :                         ncol,nlay,nbnd), intent(in   ) :: p2         !! moments of the phase function to be added
         603             :     integer,  dimension(2,nbnd),         intent(in   ) :: gpt_lims   !! Starting and ending gpoint for each band
         604             :
         605             :     integer  :: icol, ilay, igpt, ibnd, mom_lim
         606             :     real(wp) :: tau12, tauscat12
         607             :
         608           0 :     mom_lim = min(nmom1, nmom2)
         609           0 :     do ibnd = 1, nbnd
         610           0 :       do igpt = gpt_lims(1, ibnd), gpt_lims(2, ibnd)
         611           0 :         do ilay = 1, nlay
         612           0 :           do icol = 1, ncol
         613           0 :             tau12 = tau1(icol,ilay,igpt) + tau2(icol,ilay,ibnd)
         614             :             tauscat12 = &
         615             :                tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) + &
         616           0 :                tau2(icol,ilay,ibnd) * ssa2(icol,ilay,ibnd)
         617             :             !
         618             :             ! If op2 has more moments than op1 these are ignored;
         619             :             !   if it has fewer moments the higher orders are assumed to be 0
         620             :             !
         621             :             p1(1:mom_lim, icol,ilay,igpt) = &
         622           0 :                 (tau1(icol,ilay,igpt) * ssa1(icol,ilay,igpt) * p1(1:mom_lim, icol,ilay,igpt) + &
         623           0 :                  tau2(icol,ilay,ibnd) * ssa2(icol,ilay,ibnd) * p2(1:mom_lim, icol,ilay,ibnd)) / max(eps,tauscat12)
         624           0 :             ssa1(icol,ilay,igpt) = tauscat12 / max(eps,tau12)
         625           0 :             tau1(icol,ilay,igpt) = tau12
         626             :           end do
         627             :         end do
         628             :       end do
         629             :     end do
         630           0 :   end subroutine inc_nstream_by_nstream_bybnd
         631             :   ! -------------------------------------------------------------------------------------------------
         632             :   !
         633             :   ! Subsetting, meaning extracting some portion of the 3D domain
         634             :   !
         635             :   ! -------------------------------------------------------------------------------------------------
         636             :   !>
         637             :   !> Extract a subset from the first dimension (normally columns) of a 3D field.
         638             :   !>   Applicable to most variables e.g. tau, ssa, g
         639             :   !>
         640           0 :   pure subroutine extract_subset_dim1_3d(ncol, nlay, ngpt, array_in, colS, colE, array_out) &
         641             :     bind (C, name="rte_extract_subset_dim1_3d")
         642             :     integer,                             intent(in ) :: ncol, nlay, ngpt !! Array sizes
         643             :     real(wp), dimension(ncol,nlay,ngpt), intent(in ) :: array_in         !! Array to subset
         644             :     integer,                             intent(in ) :: colS, colE       !! Starting and ending index
         645             :     real(wp), dimension(colE-colS+1,&
         646             :                              nlay,ngpt), intent(out) :: array_out        !! subset of the input array
         647             :
         648             :     integer :: icol, ilay, igpt
         649           0 :     do igpt = 1, ngpt
         650           0 :       do ilay = 1, nlay
         651           0 :         do icol = colS, colE
         652           0 :           array_out(icol-colS+1, ilay, igpt) = array_in(icol, ilay, igpt)
         653             :         end do
         654             :       end do
         655             :     end do
         656             :
         657           0 :   end subroutine extract_subset_dim1_3d
         658             :   ! ---------------------------------
         659             :   !> Extract a subset from the second dimension (normally columns) of a 4D field.
         660             :   !>   Applicable to phase function moments, where the first dimension is the moment
         661           0 :   pure subroutine extract_subset_dim2_4d(nmom, ncol, nlay, ngpt, array_in, colS, colE, array_out) &
         662             :     bind (C, name="rte_extract_subset_dim2_4d")
         663             :     integer,                                  intent(in ) :: nmom, ncol, nlay, ngpt !! Array sizes
         664             :     real(wp), dimension(nmom,ncol,nlay,ngpt), intent(in ) :: array_in               !! Array to subset
         665             :     integer,                                  intent(in ) :: colS, colE             !! Starting and ending index
         666             :     real(wp), dimension(nmom,colE-colS+1,&
         667             :                                   nlay,ngpt), intent(out) :: array_out              !! subset of the input array
         668             :
         669             :     integer :: icol, ilay, igpt, imom
         670             :
         671           0 :     do igpt = 1, ngpt
         672           0 :       do ilay = 1, nlay
         673           0 :         do icol = colS, colE
         674           0 :           do imom = 1, nmom
         675           0 :             array_out(imom, icol-colS+1, ilay, igpt) = array_in(imom, icol, ilay, igpt)
         676             :           end do
         677             :         end do
         678             :       end do
         679             :     end do
         680             :
         681           0 :   end subroutine extract_subset_dim2_4d
         682             :   ! ---------------------------------
         683             :   !
         684             :   !> Extract the absorption optical thickness \(\tau_{abs} = 1 - \omega_0 \tau_{ext}\)
         685             :   !
         686           0 :   pure subroutine extract_subset_absorption_tau(ncol, nlay, ngpt, tau_in, ssa_in, &
         687           0 :                                                 colS, colE, tau_out)              &
         688             :     bind (C, name="rte_extract_subset_absorption_tau")
         689             :     integer,                             intent(in ) :: ncol, nlay, ngpt !! Array sizes
         690             :     real(wp), dimension(ncol,nlay,ngpt), intent(in ) :: tau_in, ssa_in   !! Optical thickness, single scattering albedo
         691             :     integer,                             intent(in ) :: colS, colE       !! Starting and ending index
         692             :     real(wp), dimension(colE-colS+1,&
         693             :                              nlay,ngpt), intent(out) :: tau_out          !! absorption optical thickness subset
         694             :
         695             :     integer :: icol, ilay, igpt
         696             :
         697           0 :     do igpt = 1, ngpt
         698           0 :       do ilay = 1, nlay
         699           0 :         do icol = colS, colE
         700           0 :           tau_out(icol-colS+1, ilay, igpt) = &
         701           0 :             tau_in(icol, ilay, igpt) * (1._wp - ssa_in(icol, ilay, igpt))
         702             :         end do
         703             :       end do
         704             :     end do
         705             :
         706           0 :   end subroutine extract_subset_absorption_tau
         707             : end module mo_optical_props_kernels

Generated by: LCOV version 1.14