LCOV - code coverage report
Current view: | top level - physics/rrtmgp/ext/rte-kernels - mo_rte_solver_kernels.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 200 | 391 | 51.2 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 8 | 14 | 57.1 %

              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             : !>## Numeric calculations for radiative transfer solvers
          14             : !>  - Emission/absorption (no-scattering) calculations
          15             : !>  - solver for multi-angle Gaussian quadrature
          16             : !>  - solver for a single angle, calling
          17             : !>      - source function computation (linear-in-tau)
          18             : !>      -  transport
          19             : !>  - Extinction-only calculation (direct solar beam)
          20             : !>  - Two-stream calculations:
          21             : !>    solvers for LW and SW with different boundary conditions and source functions
          22             : !>      - source function calculation for LW, SW
          23             : !>      - two-stream calculations for LW, SW (using different assumtions about phase function)
          24             : !>      - transport (adding)
          25             : !>  - Application of boundary conditions
          26             : !
          27             : ! -------------------------------------------------------------------------------------------------
          28             : module mo_rte_solver_kernels
          29             :   use,  intrinsic :: iso_c_binding
          30             :   use mo_rte_kind,      only: wp, wl
          31             :   use mo_rte_util_array,only: zero_array
          32             :   implicit none
          33             :   private
          34             :
          35             :   public :: lw_solver_noscat, lw_solver_2stream, &
          36             :             sw_solver_noscat, sw_solver_2stream
          37             :
          38             :   real(wp), parameter :: pi = acos(-1._wp)
          39             : contains
          40             :   ! -------------------------------------------------------------------------------------------------
          41             :   !
          42             :   ! Top-level longwave kernels
          43             :   !
          44             :   ! -------------------------------------------------------------------------------------------------
          45             :   !
          46             :   !>  LW fluxes, no scattering, mu (cosine of integration angle) specified by column
          47             :   !>    Does radiation calculation at user-supplied angles; converts radiances to flux
          48             :   !>    using user-supplied weights
          49             :   !
          50             :   ! ---------------------------------------------------------------
          51     1492272 :   subroutine lw_solver_noscat_oneangle(ncol, nlay, ngpt, top_at_1, D, weight,                              &
          52     1492272 :                               tau, lay_source, lev_source_inc, lev_source_dec, sfc_emis, sfc_src, &
          53     1492272 :                               incident_flux,    &
          54     1492272 :                               flux_up, flux_dn, &
          55     1492272 :                               do_broadband, broadband_up, broadband_dn, &
          56     1492272 :                               do_Jacobians, sfc_srcJac, flux_upJac,               &
          57     1492272 :                               do_rescaling, ssa, g)
          58             :     integer,                               intent(in   ) :: ncol, nlay, ngpt ! Number of columns, layers, g-points
          59             :     logical(wl),                           intent(in   ) :: top_at_1
          60             :     real(wp), dimension(ncol,       ngpt), intent(in   ) :: D            ! secant of propagation angle  []
          61             :     real(wp),                              intent(in   ) :: weight       ! quadrature weight
          62             :     real(wp), dimension(ncol,nlay,  ngpt), intent(in   ) :: tau          ! Absorption optical thickness []
          63             :     real(wp), dimension(ncol,nlay,  ngpt), intent(in   ) :: lay_source   ! Planck source at layer average temperature [W/m2]
          64             :     ! Planck source at layer edge for radiation in increasing/decreasing ilay direction
          65             :     ! lev_source_dec applies the mapping in layer i to the Planck function at layer i
          66             :     ! lev_source_inc applies the mapping in layer i to the Planck function at layer i+1
          67             :     real(wp), dimension(ncol,nlay,  ngpt), target, &
          68             :                                            intent(in   ) :: lev_source_inc, lev_source_dec
          69             :     real(wp), dimension(ncol,       ngpt), intent(in   ) :: sfc_emis     ! Surface emissivity      []
          70             :     real(wp), dimension(ncol,       ngpt), intent(in   ) :: sfc_src      ! Surface source function [W/m2]
          71             :     real(wp), dimension(ncol,       ngpt), intent(in   ) :: incident_flux! Boundary condition for flux [W/m2]
          72             :     real(wp), dimension(ncol,nlay+1,ngpt), target, &                     ! Fluxes [W/m2]
          73             :                                            intent(  out) :: flux_up, flux_dn
          74             :     !
          75             :     ! Optional variables - arrays aren't referenced if corresponding logical  == False
          76             :     !
          77             :     logical(wl),                           intent(in   ) :: do_broadband
          78             :     real(wp), dimension(ncol,nlay+1     ), intent(  out) :: broadband_up, broadband_dn ! Spectrally-integrated fluxes [W/m2]
          79             :     logical(wl),                           intent(in   ) :: do_Jacobians
          80             :     real(wp), dimension(ncol       ,ngpt), intent(in   ) :: sfc_srcJac    ! surface temperature Jacobian of surface source function [W/m2/K]
          81             :     real(wp), dimension(ncol,nlay+1     ), intent(  out) :: flux_upJac    ! surface temperature Jacobian of Radiances [W/m2-str / K]
          82             :     logical(wl),                           intent(in   ) :: do_rescaling
          83             :     real(wp), dimension(ncol,nlay  ,ngpt), intent(in   ) :: ssa, g    ! single-scattering albedo, asymmetry parameter
          84             :     ! ------------------------------------
          85             :     ! Local variables, no g-point dependency
          86             :     !
          87             :     integer                        :: icol, ilay, igpt
          88             :     integer                        :: top_level, sfc_level
          89     2984544 :     real(wp), dimension(ncol,nlay) :: tau_loc, &  ! path length (tau/mu)
          90     2984544 :                                       trans       ! transmissivity  = exp(-tau)
          91     2984544 :     real(wp), dimension(ncol,nlay) :: source_dn, source_up
          92     2984544 :     real(wp), dimension(ncol     ) :: sfc_albedo
          93             :
          94     1492272 :     real(wp), dimension(:,:,:), pointer :: lev_source_up, lev_source_dn ! Mapping increasing/decreasing indicies to up/down
          95             :
          96             :     real(wp), parameter :: pi = acos(-1._wp)
          97             :     ! loc_fluxes hold a single g-point flux if fluxes are being integrated instead of returned
          98             :     !   with spectral detail
          99             :     real(wp), dimension(ncol,nlay+1), &
         100     2984544 :                               target  :: loc_flux_up, loc_flux_dn
         101             :     ! gpt_fluxes point to calculations for the current g-point
         102     1492272 :     real(wp), dimension(:,:), pointer :: gpt_flux_up, gpt_flux_dn
         103             :     ! -------------------------------------------------------------------------------------------------
         104             :     ! Optionally, use an approximate treatment of scattering using rescaling
         105             :     !   Implemented based on the paper
         106             :     !   Tang G, et al, 2018: https://doi.org/10.1175/JAS-D-18-0014.1
         107             :     !   a) relies on rescaling of the optical parameters based on asymetry factor and single scattering albedo
         108             :     !       scaling can be computed  by scaling_1rescl
         109             :     !   b) adds adustment term based on cloud properties (lw_transport_1rescl)
         110             :     !      adustment terms is computed based on solution of the Tang equations
         111             :     !      for "linear-in-tau" internal source (not in the paper)
         112             :     !
         113             :     ! Used when approximating scattering
         114             :     !
         115             :     real(wp)                         :: ssal, wb, scaleTau
         116     2984544 :     real(wp), dimension(ncol,nlay  ) :: An, Cn
         117     2984544 :     real(wp), dimension(ncol,nlay+1) :: gpt_flux_Jac
         118             :     ! ------------------------------------
         119             :     ! Which way is up?
         120             :     ! Level Planck sources for upward and downward radiation
         121             :     ! When top_at_1, lev_source_up => lev_source_dec
         122             :     !                lev_source_dn => lev_source_inc, and vice-versa
         123     1492272 :     if(top_at_1) then
         124     1492272 :       top_level = 1
         125     1492272 :       sfc_level = nlay+1
         126     1492272 :       lev_source_up => lev_source_dec
         127     1492272 :       lev_source_dn => lev_source_inc
         128             :     else
         129           0 :       top_level = nlay+1
         130           0 :       sfc_level = 1
         131           0 :       lev_source_up => lev_source_inc
         132           0 :       lev_source_dn => lev_source_dec
         133             :     end if
         134             :
         135             :     !
         136             :     ! Integrated fluxes need zeroing
         137             :     !
         138     1492272 :     if(do_broadband) then
         139      746136 :       call zero_array(ncol, nlay+1, broadband_up )
         140      746136 :       call zero_array(ncol, nlay+1, broadband_dn )
         141             :     end if
         142     1492272 :     if(do_Jacobians) &
         143           0 :       call zero_array(ncol, nlay+1, flux_upJac )
         144             :
         145   192503088 :     do igpt = 1, ngpt
         146   191010816 :       if(do_broadband) then
         147    95505408 :         gpt_flux_up  => loc_flux_up
         148    95505408 :         gpt_flux_dn  => loc_flux_dn
         149             :       else
         150    95505408 :         gpt_flux_up  => flux_up (:,:,igpt)
         151    95505408 :         gpt_flux_dn  => flux_dn (:,:,igpt)
         152             :       end if
         153             :       !
         154             :       ! Transport is for intensity
         155             :       !   convert flux at top of domain to intensity assuming azimuthal isotropy
         156             :       !
         157  3189436416 :       gpt_flux_dn(:,top_level) = incident_flux(:,igpt)/(2._wp * pi * weight)
         158             :       !
         159             :       ! Optical path and transmission, used in source function and transport calculations
         160             :       !
         161   191010816 :       if (do_rescaling) then
         162             :         !
         163             :         ! The scaling and scaleTau terms are independent of propagation
         164             :         !   angle D and could be pre-computed if several values of D are used
         165             :         ! We re-compute them here to keep not have to localize memory use
         166             :         !
         167           0 :         do ilay = 1, nlay
         168           0 :           do icol = 1, ncol
         169           0 :             ssal = ssa(icol, ilay, igpt)
         170             :
         171             :             ! w is the layer single scattering albedo
         172             :             ! b is phase function parameter (Eq.13 of the paper)
         173             :             ! for the similarity principle scaling scheme
         174             :             ! b = (1-g)/2 (where g is phase function avergae cosine)
         175           0 :             wb = ssal*(1._wp - g(icol, ilay, igpt)) * 0.5_wp
         176             :
         177             :             ! scaleTau=1-w(1-b) is a scaling factor of the optical thickness representing
         178             :             ! the radiative transfer equation in a nonscattering form Eq(14) of the paper
         179           0 :             scaleTau = (1._wp - ssal + wb)
         180             :
         181             :             ! Cn = 0.5*wb/(1-w(1-b)) is parameter of Eq.21-22 of the Tang paper
         182             :             ! Tang paper, p.2222 advises to replace 0.5 with 0.4 based on simulations
         183           0 :             Cn(icol,ilay) = 0.4_wp*wb/scaleTau
         184             :
         185             :             ! Eqs.15, 18ab and 19 of the paper,
         186             :             ! rescaling of the optical depth multiplied by path length
         187           0 :             tau_loc(icol,ilay) = tau(icol,ilay,igpt)*D(icol,igpt)*scaleTau
         188             :           end do
         189           0 :           trans  (:,ilay) = exp(-tau_loc(:,ilay))
         190           0 :           An(:,ilay) = (1._wp-trans(:,ilay)**2)
         191             :         end do
         192             :       else
         193 18146027520 :         do ilay = 1, nlay
         194 >29980*10^7 :           tau_loc(:,ilay) = tau(:,ilay,igpt)*D(:,igpt)
         195 >29999*10^7 :           trans  (:,ilay) = exp(-tau_loc(:,ilay))
         196             :         end do
         197             :       end if
         198             :       !
         199             :       ! Source function for diffuse radiation
         200             :       !
         201             :       call lw_source_noscat(ncol, nlay, &
         202             :                             lay_source(:,:,igpt), lev_source_up(:,:,igpt), lev_source_dn(:,:,igpt), &
         203   191010816 :                             tau_loc, trans, source_dn, source_up)
         204             :       !
         205             :       ! Transport down
         206             :       !
         207   191010816 :       call lw_transport_noscat_dn(ncol, nlay, top_at_1, trans, source_dn, gpt_flux_dn)
         208             :       !
         209             :       ! Surface albedo, surface source function, reflection and emission
         210             :       !
         211  3189436416 :       sfc_albedo(:)    = 1._wp - sfc_emis(:,igpt)
         212           0 :       gpt_flux_up   (:,sfc_level) = gpt_flux_dn(:,sfc_level)*sfc_albedo(:) + &
         213  6187862016 :                                     sfc_emis(:,igpt) * sfc_src   (:,igpt)
         214   191010816 :       if(do_Jacobians) &
         215           0 :         gpt_flux_Jac(:,sfc_level)  = sfc_emis(:,igpt) * sfc_srcJac(:,igpt)
         216             :       !
         217             :       ! Transport up, or up and down again if using rescaling
         218             :       !
         219   191010816 :       if(do_rescaling) then
         220             :         call lw_transport_1rescl(ncol, nlay, top_at_1, trans,                  &
         221             :                                  source_dn, source_up,                         &
         222             :                                  gpt_flux_up, gpt_flux_dn, An, Cn, &
         223           0 :                                  do_Jacobians, gpt_flux_Jac) ! Standing in for Jacobian, i.e. rad_up_Jac(:,:,igpt), rad_dn_Jac(:,:,igpt))
         224             :       else
         225             :         call lw_transport_noscat_up(ncol, nlay, top_at_1, trans, source_up, gpt_flux_up, &
         226   191010816 :                                     do_Jacobians, gpt_flux_Jac)
         227             :       end if
         228             :
         229   191010816 :       if(do_broadband) then
         230 >15159*10^7 :         broadband_up(:,:) = broadband_up(:,:) + gpt_flux_up(:,:)
         231 >15159*10^7 :         broadband_dn(:,:) = broadband_dn(:,:) + gpt_flux_dn(:,:)
         232             :       else
         233             :         !
         234             :         ! Convert intensity to flux assuming azimuthal isotropy and quadrature weight
         235             :         !
         236 >15159*10^7 :         gpt_flux_dn(:,:)    = 2._wp * pi * weight * gpt_flux_dn(:,:)
         237 >15159*10^7 :         gpt_flux_up(:,:)    = 2._wp * pi * weight * gpt_flux_up(:,:)
         238             :       end if
         239             :       !
         240             :       ! Only broadband-integrated Jacobians are provided
         241             :       !
         242   191010816 :       if(do_Jacobians) &
         243     1492272 :           flux_upJac(:,:) =  flux_upJac(:,:) + gpt_flux_Jac(:,:)
         244             :     end do  ! g point loop
         245             :
         246     1492272 :     if(do_broadband) then
         247  1184326056 :       broadband_up(:,:) = 2._wp * pi * weight* broadband_up(:,:)
         248  1184326056 :       broadband_dn(:,:) = 2._wp * pi * weight* broadband_dn(:,:)
         249             :     end if
         250     1492272 :     if(do_Jacobians) &
         251           0 :       flux_upJac(:,:)   = 2._wp * pi * weight * flux_upJac(:,:)
         252             :
         253     1492272 :   end subroutine lw_solver_noscat_oneangle
         254             :   ! -------------------------------------------------------------------------------------------------
         255             :   !
         256             :   !> LW transport, no scattering, multi-angle quadrature
         257             :   !>   Users provide a set of weights and quadrature angles
         258             :   !>   Routine sums over single-angle solutions for each sets of angles/weights
         259             :   !
         260             :   ! ---------------------------------------------------------------
         261     1492272 :   subroutine lw_solver_noscat(ncol, nlay, ngpt, top_at_1, &
         262     1492272 :                               nmus, Ds, weights,          &
         263     1492272 :                               tau,                        &
         264     1492272 :                               lay_source, lev_source_inc, lev_source_dec, &
         265     1492272 :                               sfc_emis, sfc_src,          &
         266     1492272 :                               inc_flux,                   &
         267     1492272 :                               flux_up, flux_dn,           &
         268     1492272 :                               do_broadband, broadband_up, broadband_dn,   &
         269     1492272 :                               do_Jacobians, sfc_srcJac, flux_upJac,       &
         270     1492272 :                               do_rescaling, ssa, g) bind(C, name="rte_lw_solver_noscat")
         271             :     integer,                               intent(in   ) :: ncol, nlay, ngpt
         272             :                                                             !! Number of columns, layers, g-points
         273             :     logical(wl),                           intent(in   ) :: top_at_1
         274             :                                                             !! ilay = 1 is the top of the atmosphere?
         275             :     integer,                               intent(in   ) :: nmus
         276             :                                                             !! number of quadrature angles
         277             :     real(wp), dimension (ncol,      ngpt, &
         278             :                                     nmus), intent(in   ) :: Ds
         279             :                                                             !! quadrature secants
         280             :     real(wp), dimension(nmus),             intent(in   ) :: weights
         281             :                                                             !! quadrature weights
         282             :     real(wp), dimension(ncol,nlay,  ngpt), intent(in   ) :: tau
         283             :                                                             !! Absorption optical thickness []
         284             :     real(wp), dimension(ncol,nlay,  ngpt), intent(in   ) :: lay_source
         285             :                                                             !! Planck source at layer average temperature [W/m2]
         286             :     real(wp), dimension(ncol,nlay,  ngpt), intent(in   ) :: lev_source_inc
         287             :                                         !! Planck source at layer edge for radiation in increasing ilay direction [W/m2]
         288             :     real(wp), dimension(ncol,nlay,  ngpt), intent(in   ) :: lev_source_dec
         289             :                                         !! Planck source at layer edge for radiation in decreasing ilay direction [W/m2]
         290             :     real(wp), dimension(ncol,       ngpt), intent(in   ) :: sfc_emis
         291             :                                                             !! Surface emissivity      []
         292             :     real(wp), dimension(ncol,       ngpt), intent(in   ) :: sfc_src
         293             :                                                             !! Surface source function [W/m2]
         294             :     real(wp), dimension(ncol,       ngpt), intent(in   ) :: inc_flux
         295             :                                                             !! Incident diffuse flux, probably 0 [W/m2]
         296             :     real(wp), dimension(ncol,nlay+1,ngpt), target, &
         297             :                                            intent(  out) :: flux_up, flux_dn
         298             :                                                             !! Fluxes [W/m2]
         299             :     !
         300             :     ! Optional variables - arrays aren't referenced if corresponding logical  == False
         301             :     !
         302             :     logical(wl),                           intent(in   ) :: do_broadband
         303             :     real(wp), dimension(ncol,nlay+1     ), target, &
         304             :                                            intent(  out) :: broadband_up, broadband_dn
         305             :                                                             !! Spectrally-integrated fluxes [W/m2]
         306             :     logical(wl),                           intent(in   ) :: do_Jacobians
         307             :                                                             !! compute Jacobian with respect to surface temeprature?
         308             :     real(wp), dimension(ncol       ,ngpt), intent(in   ) :: sfc_srcJac
         309             :                                                             !! surface temperature Jacobian of surface source function [W/m2/K]
         310             :     real(wp), dimension(ncol,nlay+1     ), target, &
         311             :                                            intent(  out) :: flux_upJac
         312             :                                                             !! surface temperature Jacobian of Radiances [W/m2-str / K]
         313             :     logical(wl),                           intent(in   ) :: do_rescaling
         314             :                                                             !! Approximate treatment of scattering (10.1175/JAS-D-18-0014.1)
         315             :     real(wp), dimension(ncol,nlay  ,ngpt), intent(in   ) :: ssa, g
         316             :                                                             !! single-scattering albedo, asymmetry parameter
         317             :     ! ------------------------------------
         318             :     !
         319             :     ! Local variables - used for a single quadrature angle
         320             :     !
         321     1492272 :     real(wp), dimension(:,:,:), pointer :: this_flux_up,      this_flux_dn
         322     1492272 :     real(wp), dimension(:,:),   pointer :: this_broadband_up, this_broadband_dn, this_flux_upJac
         323             :
         324             :     integer :: imu
         325             :     ! ------------------------------------
         326             :     !
         327             :     ! For the first angle output arrays store total flux
         328             :     !
         329             :     call lw_solver_noscat_oneangle(ncol, nlay, ngpt, &
         330           0 :                           top_at_1, Ds(:,:,1), weights(1), tau, &
         331             :                           lay_source, lev_source_inc, lev_source_dec, sfc_emis, sfc_src, &
         332             :                           inc_flux,         &
         333             :                           flux_up, flux_dn, &
         334             :                           do_broadband, broadband_up, broadband_dn, &
         335             :                           do_Jacobians, sfc_srcJac, flux_upJac,     &
         336     1492272 :                           do_rescaling, ssa, g)
         337             :     !
         338             :     ! For more than one angle use local arrays
         339             :     !
         340     1492272 :     if(nmus > 1) then
         341           0 :       if(do_broadband) then
         342           0 :         allocate(this_broadband_up(ncol,nlay+1), this_broadband_dn(ncol,nlay+1))
         343             :         ! Spectrally-resolved fluxes won't be filled in so can point to caller-supplied memory
         344           0 :         this_flux_up => flux_up
         345           0 :         this_flux_dn => flux_dn
         346             :       else
         347           0 :         allocate(this_flux_up(ncol,nlay+1,ngpt), this_flux_dn(ncol,nlay+1,ngpt))
         348             :         ! Spectrally-integrated fluxes won't be filled in so can point to caller-supplied memory
         349           0 :         this_broadband_up => broadband_up
         350           0 :         this_broadband_dn => broadband_dn
         351             :       end if
         352           0 :       if(do_Jacobians) then
         353           0 :         allocate(this_flux_upJac(ncol,nlay+1))
         354             :       else
         355           0 :         this_flux_upJac => flux_upJac
         356             :       end if
         357             :     end if
         358     1492272 :     do imu = 2, nmus
         359             :       call lw_solver_noscat_oneangle(ncol, nlay, ngpt, &
         360           0 :                             top_at_1, Ds(:,:,imu), weights(imu), tau, &
         361             :                             lay_source, lev_source_inc, lev_source_dec, sfc_emis, sfc_src, &
         362             :                             inc_flux,         &
         363             :                             this_flux_up,  this_flux_dn, &
         364             :                             do_broadband, this_broadband_up, this_broadband_dn, &
         365             :                             do_Jacobians, sfc_srcJac, this_flux_upJac,         &
         366           0 :                             do_rescaling, ssa, g)
         367           0 :       if(do_broadband) then
         368           0 :         broadband_up(:,:) = broadband_up(:,:) + this_broadband_up(:,:)
         369           0 :         broadband_dn(:,:) = broadband_dn(:,:) + this_broadband_dn(:,:)
         370             :       else
         371           0 :         flux_up   (:,:,:) = flux_up   (:,:,:) + this_flux_up   (:,:,:)
         372           0 :         flux_dn   (:,:,:) = flux_dn   (:,:,:) + this_flux_dn   (:,:,:)
         373             :       end if
         374           0 :       if (do_Jacobians) &
         375     1492272 :         flux_upJac(:,:)  = flux_upJac(:,:  ) + this_flux_upJac(:,:  )
         376             :     end do
         377     1492272 :     if(nmus > 1) then
         378           0 :       if(      do_broadband) deallocate(this_broadband_up, this_broadband_dn)
         379           0 :       if(.not. do_broadband) deallocate(this_flux_up,        this_flux_dn)
         380           0 :       if(      do_Jacobians) deallocate(this_flux_upJac)
         381             :     end if
         382     1492272 :   end subroutine lw_solver_noscat
         383             :   ! -------------------------------------------------------------------------------------------------
         384             :   !
         385             :   !> Longwave two-stream calculation:
         386             :   !>   - combine RRTMGP-specific sources at levels
         387             :   !>   - compute layer reflectance, transmittance
         388             :   !>   - compute total source function at levels using linear-in-tau
         389             :   !>   - transport
         390             :   !
         391             :   ! -------------------------------------------------------------------------------------------------
         392           0 :    subroutine lw_solver_2stream (ncol, nlay, ngpt, top_at_1, &
         393           0 :                                  tau, ssa, g,                &
         394           0 :                                  lay_source, lev_source_inc, lev_source_dec, sfc_emis, sfc_src, &
         395           0 :                                  inc_flux,                   &
         396           0 :                                  flux_up, flux_dn) bind(C, name="rte_lw_solver_2stream")
         397             :     integer,                               intent(in   ) :: ncol, nlay, ngpt
         398             :                                                             !! Number of columns, layers, g-points
         399             :     logical(wl),                           intent(in   ) :: top_at_1
         400             :                                                             !! ilay = 1 is the top of the atmosphere?
         401             :     real(wp), dimension(ncol,nlay,  ngpt), intent(in   ) :: tau, ssa, g
         402             :                                                             !! Optical thickness, single-scattering albedo, asymmetry parameter []
         403             :     real(wp), dimension(ncol,nlay,  ngpt),   intent(in   ) :: lay_source
         404             :                                                             !! Planck source at layer average temperature [W/m2]
         405             :     real(wp), dimension(ncol,nlay,  ngpt), intent(in   ) :: lev_source_inc
         406             :                                           !! Planck source at layer edge for radiation in increasing ilay direction [W/m2]
         407             :     real(wp), dimension(ncol,nlay,  ngpt), intent(in   ) :: lev_source_dec
         408             :                                           !! Planck source at layer edge for radiation in decreasing ilay direction [W/m2]
         409             :     real(wp), dimension(ncol,       ngpt), intent(in   ) :: sfc_emis
         410             :                                                             !! Surface emissivity      []
         411             :     real(wp), dimension(ncol,       ngpt), intent(in   ) :: sfc_src
         412             :                                                             !! Surface source function [W/m2]
         413             :     real(wp), dimension(ncol,       ngpt), intent(in   ) :: inc_flux
         414             :                                                             !! Incident diffuse flux, probably 0 [W/m2]
         415             :     real(wp), dimension(ncol,nlay+1,ngpt), intent(  out) :: flux_up, flux_dn
         416             :                                                             !! Fluxes [W/m2]
         417             :     ! ----------------------------------------------------------------------
         418             :     integer :: igpt, top_level
         419           0 :     real(wp), dimension(ncol,nlay  ) :: Rdif, Tdif, gamma1, gamma2
         420           0 :     real(wp), dimension(ncol       ) :: sfc_albedo
         421           0 :     real(wp), dimension(ncol,nlay+1) :: lev_source
         422           0 :     real(wp), dimension(ncol,nlay  ) :: source_dn, source_up
         423           0 :     real(wp), dimension(ncol       ) :: source_sfc
         424             :     ! ------------------------------------
         425           0 :     top_level = nlay+1
         426           0 :     if(top_at_1) top_level = 1
         427           0 :     do igpt = 1, ngpt
         428             :       !
         429             :       ! RRTMGP provides source functions at each level using the spectral mapping
         430             :       !   of each adjacent layer. Combine these for two-stream calculations
         431             :       !
         432             :       call lw_combine_sources(ncol, nlay, top_at_1, &
         433             :                               lev_source_inc(:,:,igpt), lev_source_dec(:,:,igpt), &
         434           0 :                               lev_source)
         435             :       !
         436             :       ! Cell properties: reflection, transmission for diffuse radiation
         437             :       !   Coupling coefficients needed for source function
         438             :       !
         439             :       call lw_two_stream(ncol, nlay,                                 &
         440             :                          tau (:,:,igpt), ssa(:,:,igpt), g(:,:,igpt), &
         441           0 :                          gamma1, gamma2, Rdif, Tdif)
         442             :       !
         443             :       ! Source function for diffuse radiation
         444             :       !
         445             :       call lw_source_2str(ncol, nlay, top_at_1, &
         446             :                           sfc_emis(:,igpt), sfc_src(:,igpt), &
         447             :                           lay_source(:,:,igpt), lev_source, &
         448             :                           gamma1, gamma2, Rdif, Tdif, tau(:,:,igpt), &
         449           0 :                           source_dn, source_up, source_sfc)
         450             :       !
         451             :       ! Transport
         452             :       !
         453           0 :       sfc_albedo(1:ncol) = 1._wp - sfc_emis(:,igpt)
         454             :       !
         455             :       ! Boundary condition
         456             :       !
         457           0 :       flux_dn(:,top_level,igpt) = inc_flux(:,igpt)
         458             :       call adding(ncol, nlay, top_at_1,              &
         459             :                   sfc_albedo,                        &
         460             :                   Rdif, Tdif,                        &
         461             :                   source_dn, source_up, source_sfc,  &
         462           0 :                   flux_up(:,:,igpt), flux_dn(:,:,igpt))
         463             :     end do
         464             :
         465           0 :   end subroutine lw_solver_2stream
         466             :   ! -------------------------------------------------------------------------------------------------
         467             :   !
         468             :   !   Top-level shortwave kernels
         469             :   !
         470             :   ! -------------------------------------------------------------------------------------------------
         471             :   !
         472             :   !  !> Extinction-only shortwave solver i.e. solar direct beam
         473             :   !
         474             :   ! -------------------------------------------------------------------------------------------------
         475           0 :   pure subroutine sw_solver_noscat(ncol, nlay, ngpt, top_at_1, &
         476           0 :                                    tau, mu0, inc_flux_dir, flux_dir) bind(C, name="rte_sw_solver_noscat")
         477             :     integer,                               intent(in ) :: ncol, nlay, ngpt ! Number of columns, layers, g-points
         478             :                                                           !! Number of columns, layers, g-points
         479             :     logical(wl),                           intent(in ) :: top_at_1
         480             :                                                           !! ilay = 1 is the top of the atmosphere?
         481             :     real(wp), dimension(ncol,nlay,  ngpt), intent(in ) :: tau
         482             :                                                           !! Absorption optical thickness []
         483             :     real(wp), dimension(ncol,nlay       ), intent(in ) :: mu0
         484             :                                                           !! cosine of solar zenith angle
         485             :     real(wp), dimension(ncol,       ngpt), intent(in ) :: inc_flux_dir
         486             :                                                           !! Direct beam incident flux [W/m2]
         487             :     real(wp), dimension(ncol,nlay+1,ngpt), intent(out) :: flux_dir
         488             :                                                           !! Direct-beam flux, spectral [W/m2]
         489             :
         490             :     integer :: ilev, igpt
         491             :
         492             :     ! ------------------------------------
         493             :     ! Indexing into arrays for upward and downward propagation depends on the vertical
         494             :     !   orientation of the arrays (whether the domain top is at the first or last index)
         495             :     ! We write the loops out explicitly so compilers will have no trouble optimizing them.
         496             :
         497             :     ! Downward propagation
         498           0 :     if(top_at_1) then
         499             :       ! For the flux at this level, what was the previous level, and which layer has the
         500             :       !   radiation just passed through?
         501             :       ! layer index = level index - 1
         502             :       ! previous level is up (-1)
         503           0 :       do igpt = 1, ngpt
         504           0 :         flux_dir(:,    1,igpt) = inc_flux_dir(:,igpt) * mu0(:,1)
         505           0 :         do ilev = 2, nlay+1
         506           0 :           flux_dir(:,ilev,igpt) = flux_dir(:,ilev-1,igpt) * exp(-tau(:,ilev-1,igpt)/mu0(:,ilev-1))
         507             :         end do
         508             :       end do
         509             :     else
         510             :       ! layer index = level index
         511             :       ! previous level is up (+1)
         512           0 :       do igpt = 1, ngpt
         513           0 :         flux_dir(:,nlay+1,igpt) = inc_flux_dir(:,igpt) * mu0(:,nlay)
         514           0 :         do ilev = nlay, 1, -1
         515           0 :           flux_dir(:,ilev,igpt) = flux_dir(:,ilev+1,igpt) * exp(-tau(:,ilev,igpt)/mu0(:,ilev))
         516             :         end do
         517             :       end do
         518             :     end if
         519           0 :   end subroutine sw_solver_noscat
         520             :   ! -------------------------------------------------------------------------------------------------
         521             :   !
         522             :   !> Shortwave two-stream calculation:
         523             :   !>   compute layer reflectance, transmittance
         524             :   !>   compute solar source function for diffuse radiation
         525             :   !>   transport
         526             :   !
         527             :   ! -------------------------------------------------------------------------------------------------
         528      778526 :   subroutine sw_solver_2stream (ncol, nlay, ngpt, top_at_1,  &
         529      778526 :                                 tau, ssa, g, mu0,           &
         530      778526 :                                 sfc_alb_dir, sfc_alb_dif,   &
         531      778526 :                                             inc_flux_dir,   &
         532      778526 :                                 flux_up, flux_dn, flux_dir, &
         533      778526 :                                 has_dif_bc, inc_flux_dif,   &
         534      778526 :                                 do_broadband, broadband_up, &
         535      778526 :                                 broadband_dn, broadband_dir) bind(C, name="rte_sw_solver_2stream")
         536             :     integer,                               intent(in   ) :: ncol, nlay, ngpt
         537             :                                                             !! Number of columns, layers, g-points
         538             :     logical(wl),                           intent(in   ) :: top_at_1
         539             :                                                             !! ilay = 1 is the top of the atmosphere?
         540             :     real(wp), dimension(ncol,nlay,  ngpt), intent(in   ) :: tau, ssa, g
         541             :                                                             !! Optical thickness, single-scattering albedo, asymmetry parameter []
         542             :     real(wp), dimension(ncol,nlay       ), intent(in   ) :: mu0
         543             :                                                             !! cosine of solar zenith angle
         544             :     real(wp), dimension(ncol,       ngpt), intent(in   ) :: sfc_alb_dir, sfc_alb_dif
         545             :                                                             !! Spectral surface albedo for direct and diffuse radiation
         546             :     real(wp), dimension(ncol,       ngpt), intent(in   ) :: inc_flux_dir
         547             :                                                             !! Direct beam incident flux
         548             :     real(wp), dimension(ncol,nlay+1,ngpt), target, &
         549             :                                            intent(  out) :: flux_up, flux_dn, flux_dir
         550             :                                                             !! Fluxes [W/m2]
         551             :     logical(wl),                           intent(in   ) :: has_dif_bc
         552             :                                                             !! Is a boundary condition for diffuse flux supplied?
         553             :     real(wp), dimension(ncol,       ngpt), intent(in   ) :: inc_flux_dif
         554             :                                                             !! Boundary condition for diffuse flux [W/m2]
         555             :     logical(wl),                           intent(in   ) :: do_broadband
         556             :                                                             !! Provide broadband-integrated, not spectrally-resolved, fluxes?
         557             :     real(wp), dimension(ncol,nlay+1     ), intent(  out) :: broadband_up, broadband_dn, broadband_dir
         558             :                                                             !! Broadband integrated fluxes
         559             :     ! -------------------------------------------
         560             :     integer :: igpt, top_level, top_layer
         561     1557052 :     real(wp), dimension(ncol,nlay  )  :: Rdif, Tdif
         562     1557052 :     real(wp), dimension(ncol,nlay  )  :: source_up, source_dn
         563     1557052 :     real(wp), dimension(ncol       )  :: source_srf
         564             :     ! loc_fluxes hold a single g-point flux if fluxes are being integrated instead of returned
         565             :     !   with spectral detail
         566             :     real(wp), dimension(ncol,nlay+1), &
         567     1557052 :                               target  :: loc_flux_up, loc_flux_dn, loc_flux_dir
         568             :     ! gpt_fluxes point to calculations for the current g-point
         569      778526 :     real(wp), dimension(:,:), pointer :: gpt_flux_up, gpt_flux_dn, gpt_flux_dir
         570             :     ! ------------------------------------
         571      778526 :     if(top_at_1) then
         572             :       top_level = 1
         573             :       top_layer = 1
         574             :     else
         575           0 :       top_level  = nlay+1
         576           0 :       top_layer  = nlay
         577             :     end if
         578             :     !
         579             :     ! Integrated fluxes need zeroing
         580             :     !
         581      778526 :     if(do_broadband) then
         582      389263 :       call zero_array(ncol, nlay+1, broadband_up )
         583      389263 :       call zero_array(ncol, nlay+1, broadband_dn )
         584      389263 :       call zero_array(ncol, nlay+1, broadband_dir)
         585             :     end if
         586             :
         587    87973438 :     do igpt = 1, ngpt
         588    87194912 :       if(do_broadband) then
         589    43597456 :         gpt_flux_up  => loc_flux_up
         590    43597456 :         gpt_flux_dn  => loc_flux_dn
         591    43597456 :         gpt_flux_dir => loc_flux_dir
         592             :       else
         593    43597456 :         gpt_flux_up  => flux_up (:,:,igpt)
         594    43597456 :         gpt_flux_dn  => flux_dn (:,:,igpt)
         595    43597456 :         gpt_flux_dir => flux_dir(:,:,igpt)
         596             :       end if
         597             :       !
         598             :       ! Boundary conditions direct beam...
         599             :       !
         600  1399006112 :       gpt_flux_dir(:,top_level) = inc_flux_dir(:,igpt) * mu0(:,top_layer)
         601             :       !
         602             :       ! ... and diffuse field, using 0 if no BC is provided
         603             :       !
         604    87194912 :       if(has_dif_bc) then
         605           0 :         gpt_flux_dn(:,top_level) = inc_flux_dif(:,igpt)
         606             :       else
         607  1399006112 :         gpt_flux_dn(:,top_level) = 0._wp
         608             :       end if
         609             :       !
         610             :       ! Cell properties: transmittance and reflectance for diffuse radiation
         611             :       !   Direct-beam and source for diffuse radiation
         612             :       !
         613             :       call sw_dif_and_source(ncol, nlay, top_at_1, mu0, sfc_alb_dir(:,igpt), &
         614             :                              tau(:,:,igpt), ssa(:,:,igpt), g(:,:,igpt),      &
         615             :                              Rdif, Tdif, source_dn, source_up, source_srf,  &
         616    87194912 :                              gpt_flux_dir)
         617             :       !
         618             :       ! Transport
         619             :       !
         620             :       call adding(ncol, nlay, top_at_1,            &
         621             :                   sfc_alb_dif(:,igpt), Rdif, Tdif, &
         622    87194912 :                   source_dn, source_up, source_srf, gpt_flux_up, gpt_flux_dn)
         623             :       !
         624             :       ! adding() computes only diffuse flux; flux_dn is total
         625             :       !
         626    87973438 :       if(do_broadband) then
         627 66496387776 :         broadband_up (:,:) = broadband_up (:,:) + gpt_flux_up (:,:)
         628 66496387776 :         broadband_dn (:,:) = broadband_dn (:,:) + gpt_flux_dn (:,:) + gpt_flux_dir(:,:)
         629 66496387776 :         broadband_dir(:,:) = broadband_dir(:,:)                     + gpt_flux_dir(:,:)
         630             :       else
         631 >13299*10^7 :         gpt_flux_dn(:,:) =                        gpt_flux_dn (:,:) + gpt_flux_dir(:,:)
         632             :       end if
         633             :     end do
         634      778526 :   end subroutine sw_solver_2stream
         635             :   ! -------------------------------------------------------------------------------------------------
         636             :   !
         637             :   !   Lower-level longwave kernels
         638             :   !
         639             :   ! -------------------------------------------------------------------------------------------------
         640             :   !
         641             :   ! Compute LW source function for upward and downward emission at levels using linear-in-tau assumption
         642             :   ! See Clough et al., 1992, doi: 10.1029/92JD01419, Eq 13
         643             :   !
         644             :   ! ---------------------------------------------------------------
         645   191010816 :   subroutine lw_source_noscat(ncol, nlay, lay_source, lev_source_up, lev_source_dn, tau, trans, &
         646   191010816 :                               source_dn, source_up)
         647             :     integer,                         intent(in) :: ncol, nlay
         648             :     real(wp), dimension(ncol, nlay), intent(in) :: lay_source, & ! Planck source at layer center
         649             :                                                    lev_source_up, & ! Planck source at levels (layer edges),
         650             :                                                    lev_source_dn, & !   increasing/decreasing layer index
         651             :                                                    tau,        & ! Optical path (tau/mu)
         652             :                                                    trans         ! Transmissivity (exp(-tau))
         653             :     real(wp), dimension(ncol, nlay), intent(out):: source_dn, source_up
         654             :                                                                    ! Source function at layer edges
         655             :                                                                    ! Down at the bottom of the layer, up at the top
         656             :     ! --------------------------------
         657             :     integer             :: icol, ilay
         658             :     real(wp)            :: fact
         659             :     real(wp), parameter :: tau_thresh = sqrt(sqrt(epsilon(tau)))
         660             :     ! ---------------------------------------------------------------
         661 18146027520 :     do ilay = 1, nlay
         662 >29999*10^7 :       do icol = 1, ncol
         663             :       !
         664             :       ! Weighting factor. Use 2nd order series expansion when rounding error (~tau^2)
         665             :       !   is of order epsilon (smallest difference from 1. in working precision)
         666             :       !   Thanks to Peter Blossey
         667             :       !   Updated to 3rd order series and lower threshold based on suggestion from Dmitry Alexeev (Nvidia)
         668             :       !
         669 >28185*10^7 :       if(tau(icol, ilay) > tau_thresh) then
         670 >21727*10^7 :         fact = (1._wp - trans(icol,ilay))/tau(icol,ilay) - trans(icol,ilay)
         671             :       else
         672 64576670882 :         fact = tau(icol, ilay) * (0.5_wp + tau(icol, ilay) * (- 1._wp/3._wp + tau(icol, ilay) * 1._wp/8._wp ) )
         673             :       end if
         674             :       !
         675             :       ! Equation below is developed in Clough et al., 1992, doi:10.1029/92JD01419, Eq 13
         676             :       !
         677             :       source_dn(icol,ilay) = (1._wp - trans(icol,ilay)) * lev_source_dn(icol,ilay) + &
         678 >28185*10^7 :                               2._wp * fact * (lay_source(icol,ilay) - lev_source_dn(icol,ilay))
         679             :       source_up(icol,ilay) = (1._wp - trans(icol,ilay)) * lev_source_up(icol,ilay  ) + &
         680 >29980*10^7 :                               2._wp * fact * (lay_source(icol,ilay) - lev_source_up(icol,ilay))
         681             :       end do
         682             :     end do
         683   191010816 :   end subroutine lw_source_noscat
         684             :   ! -------------------------------------------------------------------------------------------------
         685             :   !
         686             :   ! Longwave no-scattering transport - separate routines for up and down
         687             :   !
         688             :   ! -------------------------------------------------------------------------------------------------
         689   191010816 :   subroutine lw_transport_noscat_dn(ncol, nlay, top_at_1,     &
         690   191010816 :                                    trans, source_dn, radn_dn)
         691             :     integer,                          intent(in   ) :: ncol, nlay ! Number of columns, layers, g-points
         692             :     logical(wl),                      intent(in   ) :: top_at_1   !
         693             :     real(wp), dimension(ncol,nlay  ), intent(in   ) :: trans      ! transmissivity = exp(-tau)
         694             :     real(wp), dimension(ncol,nlay  ), intent(in   ) :: source_dn  ! Diffuse radiation emitted by the layer
         695             :     real(wp), dimension(ncol,nlay+1), intent(inout) :: radn_dn    ! Radiances [W/m2-str] Top level must contain incident flux boundary condition
         696             :
         697             :     ! ---------------------------------------------------
         698             :     ! Local variables
         699             :     integer :: ilev
         700             :     ! ---------------------------------------------------
         701   191010816 :     if(top_at_1) then
         702             :       !
         703             :       ! Top of domain is index 1
         704             :       !
         705 18146027520 :       do ilev = 2, nlay+1
         706 >29999*10^7 :         radn_dn(:,ilev) = trans(:,ilev-1)*radn_dn(:,ilev-1) + source_dn(:,ilev-1)
         707             :       end do
         708             :     else
         709             :       !
         710             :       ! Top of domain is index nlay+1
         711             :       !
         712           0 :       do ilev = nlay, 1, -1
         713           0 :         radn_dn(:,ilev) = trans(:,ilev  )*radn_dn(:,ilev+1) + source_dn(:,ilev)
         714             :       end do
         715             :     end if
         716   191010816 :   end subroutine lw_transport_noscat_dn
         717             :   ! -------------------------------------------------------------------------------------------------
         718   191010816 :   subroutine lw_transport_noscat_up(ncol, nlay, top_at_1,     &
         719   191010816 :                                    trans, source_up, radn_up, do_Jacobians, radn_upJac)
         720             :     integer,                          intent(in   ) :: ncol, nlay ! Number of columns, layers, g-points
         721             :     logical(wl),                      intent(in   ) :: top_at_1   !
         722             :     real(wp), dimension(ncol,nlay  ), intent(in   ) :: trans      ! transmissivity = exp(-tau)
         723             :     real(wp), dimension(ncol,nlay  ), intent(in   ) :: source_up  ! Diffuse radiation emitted by the layer
         724             :     real(wp), dimension(ncol,nlay+1), intent(inout) :: radn_up    ! Radiances [W/m2-str] Top level must contain incident flux boundary condition
         725             :     logical(wl),                      intent(in   ) :: do_Jacobians
         726             :     real(wp), dimension(ncol,nlay+1), intent(inout) :: radn_upJac       ! surface temperature Jacobian of Radiances [W/m2-str / K]
         727             :
         728             :     ! ---------------------------------------------------
         729             :     ! Local variables
         730             :     integer :: ilev
         731             :     ! ---------------------------------------------------
         732   191010816 :     if(top_at_1) then
         733             :       !
         734             :       ! Top of domain is index 1
         735             :       !
         736             :       ! Upward propagation
         737 18146027520 :       do ilev = nlay, 1, -1
         738 >31776*10^7 :         radn_up     (:,ilev) = trans(:,ilev  )*radn_up   (:,ilev+1) + source_up(:,ilev)
         739 17955016704 :         if(do_Jacobians) &
         740   191010816 :           radn_upJac(:,ilev) = trans(:,ilev  )*radn_upJac(:,ilev+1)
         741             :       end do
         742             :     else
         743             :       !
         744             :       ! Top of domain is index nlay+1
         745             :       !
         746             :       ! Upward propagation
         747           0 :       do ilev = 2, nlay+1
         748           0 :         radn_up     (:,ilev) = trans(:,ilev-1) * radn_up   (:,ilev-1) +  source_up(:,ilev-1)
         749           0 :         if(do_Jacobians) &
         750           0 :           radn_upJac(:,ilev) = trans(:,ilev-1) * radn_upJac(:,ilev-1)
         751             :       end do
         752             :     end if
         753   191010816 :   end subroutine lw_transport_noscat_up
         754             :   ! -------------------------------------------------------------------------------------------------
         755             :   ! Upward and (second) downward transport for re-scaled longwave solution
         756             :   !   adds adjustment factor based on cloud properties
         757             :   !
         758             :   !   implementation notice:
         759             :   !       the adjustmentFactor computation can be skipped where Cn <= epsilon
         760             :   ! -------------------------------------------------------------------------------------------------
         761           0 :   subroutine lw_transport_1rescl(ncol, nlay, top_at_1, &
         762           0 :                                  trans, source_dn, source_up, &
         763           0 :                                  radn_up, radn_dn, An, Cn,&
         764           0 :                                  do_Jacobians, radn_up_Jac)
         765             :     integer,                          intent(in   ) :: ncol, nlay ! Number of columns, layers, g-points
         766             :     logical(wl),                      intent(in   ) :: top_at_1   !
         767             :     real(wp), dimension(ncol,nlay  ), intent(in   ) :: trans      ! transmissivity = exp(-tau)
         768             :     real(wp), dimension(ncol,nlay  ), intent(in   ) :: source_dn, &
         769             :                                                        source_up  ! Diffuse radiation emitted by the layer
         770             :     real(wp), dimension(ncol,nlay+1), intent(inout) :: radn_up    ! Radiances [W/m2-str]
         771             :     real(wp), dimension(ncol,nlay+1), intent(inout) :: radn_dn    !Top level must contain incident flux boundary condition
         772             :     real(wp), dimension(ncol,nlay),   intent(in   ) :: An, Cn
         773             :     logical(wl),                      intent(in   ) :: do_Jacobians
         774             :     real(wp), dimension(ncol,nlay+1), intent(inout) :: radn_up_Jac ! Surface temperature Jacobians [W/m2-str/K]
         775             :     !
         776             :     ! We could in principle compute a downwelling Jacobian too, but it's small
         777             :     !   (only a small proportion of LW is scattered) and it complicates code and the API,
         778             :     !   so we will not
         779             :     !
         780             :
         781             :     ! Local variables
         782             :     integer :: ilev, icol
         783             :     ! ---------------------------------------------------
         784             :     real(wp) :: adjustmentFactor
         785           0 :     if(top_at_1) then
         786             :       !
         787             :       ! Top of domain is index 1
         788             :       !
         789             :       ! Upward propagation
         790             :       ! adjustment factor is obtained as a solution of 18b of the Tang paper
         791             :       ! eqvivalent to Eq.20 of the Tang paper but for linear-in-tau source
         792           0 :       do ilev = nlay, 1, -1
         793           0 :         do icol=1,ncol
         794           0 :           adjustmentFactor = Cn(icol,ilev)*( An(icol,ilev)*radn_dn(icol,ilev) - &
         795           0 :                    trans(icol,ilev)*source_dn(icol,ilev) - source_up(icol,ilev) )
         796           0 :           radn_up (icol,ilev) = trans(icol,ilev)*radn_up(icol,ilev+1) + source_up(icol,ilev) + &
         797           0 :                                 adjustmentFactor
         798             :         end do
         799           0 :         if(do_Jacobians) &
         800           0 :           radn_up_Jac(:,ilev) = trans(:,ilev)*radn_up_Jac(:,ilev+1)
         801             :       end do
         802             :       ! Downward propagation
         803             :       ! radn_dn_Jac(:,1) = 0._wp
         804             :       ! adjustment factor is obtained as a solution of 19 of the Tang paper
         805             :       ! eqvivalent to Eq.21 of the Tang paper but for linear-in-tau source
         806           0 :       do ilev = 1, nlay
         807             :         ! radn_dn_Jac(:,ilev+1) = trans(:,ilev)*radn_dn_Jac(:,ilev)
         808           0 :         do icol=1,ncol
         809           0 :             adjustmentFactor = Cn(icol,ilev)*( An(icol,ilev)*radn_up(icol,ilev) - &
         810           0 :                      trans(icol,ilev)*source_up(icol,ilev) - source_dn(icol,ilev) )
         811           0 :             radn_dn(icol,ilev+1) = trans(icol,ilev)*radn_dn(icol,ilev) + source_dn(icol,ilev) + &
         812           0 :                                    adjustmentFactor
         813             :             ! adjustmentFactor         = Cn(icol,ilev)*An(icol,ilev)*radn_up_Jac(icol,ilev)
         814             :             ! radn_dn_Jac(icol,ilev+1) = radn_dn_Jac(icol,ilev+1) + adjustmentFactor
         815             :         enddo
         816             :       end do
         817             :     else
         818             :       !
         819             :       ! Top of domain is index nlay+1
         820             :       !
         821             :       ! Upward propagation
         822             :       ! adjustment factor is obtained as a solution of 18b of the Tang paper
         823             :       ! eqvivalent to Eq.20 of the Tang paper but for linear-in-tau source
         824           0 :       do ilev = 1, nlay
         825           0 :         radn_up      (:,ilev+1) = trans(:,ilev) * radn_up    (:,ilev) +  source_up(:,ilev)
         826           0 :         do icol=1,ncol
         827           0 :             adjustmentFactor = Cn(icol,ilev)*( An(icol,ilev)*radn_dn(icol,ilev+1) - &
         828           0 :                      trans(icol,ilev)*source_dn(icol,ilev) - source_up(icol,ilev) )
         829             :             radn_up(icol,ilev+1) = trans(icol,ilev)*radn_up(icol,ilev) +  source_up(icol,ilev) + &
         830           0 :                                    adjustmentFactor
         831             :         enddo
         832           0 :         if(do_Jacobians) &
         833           0 :           radn_up_Jac(:,ilev+1) = trans(:,ilev) * radn_up_Jac(:,ilev)
         834             :       end do
         835             :
         836             :       ! Downward propagation
         837             :       ! adjustment factor is obtained as a solution of 19 of the Tang paper
         838             :       ! eqvivalent to Eq.21 of the Tang paper but for linear-in-tau source
         839             :       ! radn_dn_Jac(:,nlay+1) = 0._wp
         840           0 :       do ilev = nlay, 1, -1
         841             :         ! radn_dn_Jac(:,ilev) = trans(:,ilev)*radn_dn_Jac(:,ilev+1)
         842           0 :         do icol=1,ncol
         843           0 :             adjustmentFactor = Cn(icol,ilev)*( An(icol,ilev)*radn_up(icol,ilev) - &
         844           0 :                      trans(icol,ilev)*source_up(icol,ilev) - source_dn(icol,ilev) )
         845           0 :             radn_dn(icol,ilev)  = trans(icol,ilev)*radn_dn(icol,ilev+1) + source_dn(icol,ilev) + &
         846           0 :                                   adjustmentFactor
         847             :             ! adjustmentFactor    = Cn(icol,ilev)*An(icol,ilev)*radn_up_Jac(icol,ilev)
         848             :             ! radn_dn_Jac(icol,ilev) = radn_dn_Jac(icol,ilev) + adjustmentFactor
         849             :         enddo
         850             :       end do
         851             :     end if
         852           0 :   end subroutine lw_transport_1rescl
         853             : ! -------------------------------------------------------------------------------------------------
         854             :   !
         855             :   ! Longwave two-stream solutions to diffuse reflectance and transmittance for a layer
         856             :   !    with optical depth tau, single scattering albedo w0, and asymmetery parameter g.
         857             :   !
         858             :   ! Equations are developed in Meador and Weaver, 1980,
         859             :   !    doi:10.1175/1520-0469(1980)037<0630:TSATRT>2.0.CO;2
         860             :   !
         861             :   ! -------------------------------------------------------------------------------------------------
         862           0 :   pure subroutine lw_two_stream(ncol, nlay, tau, w0, g, &
         863           0 :                                 gamma1, gamma2, Rdif, Tdif)
         864             :     integer,                        intent(in)  :: ncol, nlay
         865             :     real(wp), dimension(ncol,nlay), intent(in)  :: tau, w0, g
         866             :     real(wp), dimension(ncol,nlay), intent(out) :: gamma1, gamma2, Rdif, Tdif
         867             :
         868             :     ! -----------------------
         869             :     integer  :: i, j
         870             :
         871             :     ! Variables used in Meador and Weaver
         872           0 :     real(wp) :: k(ncol)
         873             :
         874             :     ! Ancillary variables
         875           0 :     real(wp) :: RT_term(ncol)
         876           0 :     real(wp) :: exp_minusktau(ncol), exp_minus2ktau(ncol)
         877             :
         878             :     real(wp), parameter :: LW_diff_sec = 1.66  ! 1./cos(diffusivity angle)
         879             :     ! ---------------------------------
         880           0 :     do j = 1, nlay
         881           0 :       do i = 1, ncol
         882             :         !
         883             :         ! Coefficients differ from SW implementation because the phase function is more isotropic
         884             :         !   Here we follow Fu et al. 1997, doi:10.1175/1520-0469(1997)054<2799:MSPITI>2.0.CO;2
         885             :         !   and use a diffusivity sec of 1.66
         886             :         !
         887           0 :         gamma1(i,j)= LW_diff_sec * (1._wp - 0.5_wp * w0(i,j) * (1._wp + g(i,j))) ! Fu et al. Eq 2.9
         888           0 :         gamma2(i,j)= LW_diff_sec *          0.5_wp * w0(i,j) * (1._wp - g(i,j))  ! Fu et al. Eq 2.10
         889             :         ! Eq 18;  k = SQRT(gamma1**2 - gamma2**2), limited below to avoid div by 0.
         890             :         !   k = 0 for isotropic, conservative scattering; this lower limit on k
         891             :         !   gives relative error with respect to conservative solution
         892             :         !   of < 0.1% in Rdif down to tau = 10^-9
         893           0 :         k(i) = sqrt(max((gamma1(i,j) - gamma2(i,j)) * (gamma1(i,j) + gamma2(i,j)), 1.e-12_wp))
         894             :       end do
         895             :
         896             :       ! Written to encourage vectorization of exponential
         897           0 :       exp_minusktau(1:ncol) = exp(-tau(1:ncol,j)*k(1:ncol))
         898             :
         899             :       !
         900             :       ! Diffuse reflection and transmission
         901             :       !
         902           0 :       do i = 1, ncol
         903           0 :         exp_minus2ktau(i) = exp_minusktau(i) * exp_minusktau(i)
         904             :
         905             :         ! Refactored to avoid rounding errors when k, gamma1 are of very different magnitudes
         906             :         RT_term(i) = 1._wp / (k     (i  ) * (1._wp + exp_minus2ktau(i))  + &
         907           0 :                               gamma1(i,j) * (1._wp - exp_minus2ktau(i)) )
         908             :
         909             :         ! Equation 25
         910           0 :         Rdif(i,j) = RT_term(i) * gamma2(i,j) * (1._wp - exp_minus2ktau(i))
         911             :
         912             :         ! Equation 26
         913           0 :         Tdif(i,j) = RT_term(i) * 2._wp * k(i) * exp_minusktau(i)
         914             :       end do
         915             :
         916             :     end do
         917           0 :   end subroutine lw_two_stream
         918             :   ! -------------------------------------------------------------------------------------------------
         919             :   !
         920             :   ! Source function combination
         921             :   ! RRTMGP provides two source functions at each level
         922             :   !   using the spectral mapping from each of the adjascent layers.
         923             :   !   Need to combine these for use in two-stream calculation.
         924             :   !
         925             :   ! -------------------------------------------------------------------------------------------------
         926           0 :   subroutine lw_combine_sources(ncol, nlay, top_at_1, &
         927           0 :                                 lev_src_inc, lev_src_dec, lev_source)
         928             :     integer,                           intent(in ) :: ncol, nlay
         929             :     logical(wl),                       intent(in ) :: top_at_1
         930             :     real(wp), dimension(ncol, nlay  ), intent(in ) :: lev_src_inc, lev_src_dec
         931             :     real(wp), dimension(ncol, nlay+1), intent(out) :: lev_source
         932             :
         933             :     integer :: icol, ilay
         934             :     ! ---------------------------------------------------------------
         935           0 :     ilay = 1
         936           0 :     do icol = 1,ncol
         937           0 :       lev_source(icol, ilay) =        lev_src_dec(icol, ilay)
         938             :     end do
         939           0 :     do ilay = 2, nlay
         940           0 :       do icol = 1,ncol
         941           0 :         lev_source(icol, ilay) = sqrt(lev_src_dec(icol, ilay) * &
         942           0 :                                       lev_src_inc(icol, ilay-1))
         943             :       end do
         944             :     end do
         945           0 :     ilay = nlay+1
         946           0 :     do icol = 1,ncol
         947           0 :       lev_source(icol, ilay) =        lev_src_inc(icol, ilay-1)
         948             :     end do
         949             :
         950           0 :   end subroutine lw_combine_sources
         951             :   ! ---------------------------------------------------------------
         952             :   !
         953             :   ! Compute LW source function for upward and downward emission at levels using linear-in-tau assumption
         954             :   !   This version straight from ECRAD
         955             :   !   Source is provided as W/m2-str; factor of pi converts to flux units
         956             :   !
         957             :   ! ---------------------------------------------------------------
         958           0 :   subroutine lw_source_2str(ncol, nlay, top_at_1,   &
         959           0 :                             sfc_emis, sfc_src,      &
         960           0 :                             lay_source, lev_source, &
         961           0 :                             gamma1, gamma2, rdif, tdif, tau, source_dn, source_up, source_sfc) &
         962             :                             bind (C, name="rte_lw_source_2str")
         963             :     integer,                         intent(in) :: ncol, nlay
         964             :     logical(wl),                     intent(in) :: top_at_1
         965             :     real(wp), dimension(ncol      ), intent(in) :: sfc_emis, sfc_src
         966             :     real(wp), dimension(ncol, nlay), intent(in) :: lay_source,    & ! Planck source at layer center
         967             :                                                    tau,           & ! Optical depth (tau)
         968             :                                                    gamma1, gamma2,& ! Coupling coefficients
         969             :                                                    rdif, tdif       ! Layer reflectance and transmittance
         970             :     real(wp), dimension(ncol, nlay+1), target, &
         971             :                                      intent(in)  :: lev_source       ! Planck source at layer edges
         972             :     real(wp), dimension(ncol, nlay), intent(out) :: source_dn, source_up
         973             :     real(wp), dimension(ncol      ), intent(out) :: source_sfc      ! Source function for upward radation at surface
         974             :
         975             :     integer             :: icol, ilay
         976             :     real(wp)            :: Z, Zup_top, Zup_bottom, Zdn_top, Zdn_bottom
         977           0 :     real(wp), dimension(:), pointer :: lev_source_bot, lev_source_top
         978             :     ! ---------------------------------------------------------------
         979           0 :     do ilay = 1, nlay
         980           0 :       if(top_at_1) then
         981           0 :         lev_source_top => lev_source(:,ilay)
         982           0 :         lev_source_bot => lev_source(:,ilay+1)
         983             :       else
         984           0 :         lev_source_top => lev_source(:,ilay+1)
         985           0 :         lev_source_bot => lev_source(:,ilay)
         986             :       end if
         987           0 :       do icol = 1, ncol
         988           0 :         if (tau(icol,ilay) > 1.0e-8_wp) then
         989             :           !
         990             :           ! Toon et al. (JGR 1989) Eqs 26-27
         991             :           !
         992           0 :           Z = (lev_source_bot(icol)-lev_source_top(icol)) / (tau(icol,ilay)*(gamma1(icol,ilay)+gamma2(icol,ilay)))
         993           0 :           Zup_top        =  Z + lev_source_top(icol)
         994           0 :           Zup_bottom     =  Z + lev_source_bot(icol)
         995           0 :           Zdn_top        = -Z + lev_source_top(icol)
         996           0 :           Zdn_bottom     = -Z + lev_source_bot(icol)
         997           0 :           source_up(icol,ilay) = pi * (Zup_top    - rdif(icol,ilay) * Zdn_top    - tdif(icol,ilay) * Zup_bottom)
         998           0 :           source_dn(icol,ilay) = pi * (Zdn_bottom - rdif(icol,ilay) * Zup_bottom - tdif(icol,ilay) * Zdn_top)
         999             :         else
        1000           0 :           source_up(icol,ilay) = 0._wp
        1001           0 :           source_dn(icol,ilay) = 0._wp
        1002             :         end if
        1003             :       end do
        1004             :     end do
        1005           0 :     do icol = 1, ncol
        1006           0 :       source_sfc(icol) = pi * sfc_emis(icol) * sfc_src(icol)
        1007             :     end do
        1008           0 :   end subroutine lw_source_2str
        1009             :   ! -------------------------------------------------------------------------------------------------
        1010             :   !
        1011             :   !   Lower-level shortwave kernels
        1012             :   !
        1013             :   ! -------------------------------------------------------------------------------------------------
        1014             :   !
        1015             :   ! Two-stream solutions to diffuse reflectance and transmittance for a layer
        1016             :   !    with optical depth tau, single scattering albedo w0, and asymmetery parameter g.
        1017             :   ! Direct reflectance and transmittance used to compute direct beam source for diffuse radiation
        1018             :   !   in layers and at surface; report direct beam as a byproduct
        1019             :   ! Computing the direct-beam source for diffuse radiation at the same time as R and T for
        1020             :   !   direct radiation reduces memory traffic and use.
        1021             :   !
        1022             :   ! Equations are developed in Meador and Weaver, 1980,
        1023             :   !    doi:10.1175/1520-0469(1980)037<0630:TSATRT>2.0.CO;2
        1024             :   !
        1025             :   ! -------------------------------------------------------------------------------------------------
        1026    87194912 :   pure subroutine sw_dif_and_source(ncol, nlay, top_at_1, mu0, sfc_albedo, &
        1027    87194912 :                                     tau, w0, g,  &
        1028    87194912 :                                     Rdif, Tdif, source_dn, source_up, source_sfc, flux_dn_dir) bind (C, name="rte_sw_source_dir")
        1029             :     integer,                          intent(in   ) :: ncol, nlay
        1030             :     logical(wl),                      intent(in   ) :: top_at_1
        1031             :     real(wp), dimension(ncol       ), intent(in   ) :: sfc_albedo          ! surface albedo for direct radiation
        1032             :     real(wp), dimension(ncol,nlay  ), intent(in   ) :: tau, w0, g, mu0
        1033             :     real(wp), dimension(ncol,nlay  ), intent(  out) :: Rdif, Tdif, source_dn, source_up
        1034             :     real(wp), dimension(ncol       ), intent(  out) :: source_sfc ! Source function for upward radation at surface
        1035             :     real(wp), dimension(ncol,nlay+1), target, &
        1036             :                                       intent(inout) :: flux_dn_dir ! Direct beam flux
        1037             :
        1038             :     ! -----------------------
        1039             :     integer  :: i, j
        1040             :
        1041             :     ! Variables used in Meador and Weaver
        1042             :     real(wp) :: gamma1, gamma2, gamma3, gamma4, alpha1, alpha2
        1043             :
        1044             :
        1045             :     ! Ancillary variables
        1046             :     real(wp), parameter :: min_k = 1.e4_wp * epsilon(1._wp) ! Suggestion from Chiel van Heerwaarden
        1047             :     real(wp) :: k, exp_minusktau, k_mu, k_gamma3, k_gamma4
        1048             :     real(wp) :: RT_term, exp_minus2ktau
        1049             :     real(wp) :: Rdir, Tdir, Tnoscat
        1050    87194912 :     real(wp), pointer, dimension(:) :: dir_flux_inc, dir_flux_trans
        1051             :     integer  :: lay_index
        1052             :     real(wp) :: tau_s, w0_s, g_s, mu0_s
        1053             :     ! ---------------------------------
        1054             :
        1055  8283516640 :     do j = 1, nlay
        1056  8196321728 :       if(top_at_1) then
        1057  8196321728 :         lay_index      =  j
        1058  8196321728 :         dir_flux_inc   => flux_dn_dir(:,lay_index  )
        1059  8196321728 :         dir_flux_trans => flux_dn_dir(:,lay_index+1)
        1060             :       else
        1061           0 :         lay_index      =  nlay-j+1
        1062           0 :         dir_flux_inc   => flux_dn_dir(:,lay_index+1)
        1063           0 :         dir_flux_trans => flux_dn_dir(:,lay_index  )
        1064             :       end if
        1065             :
        1066 >13159*10^7 :       do i = 1, ncol
        1067             :         !
        1068             :         ! Scalars
        1069             :         !
        1070 >12331*10^7 :         tau_s = tau(i, lay_index)
        1071 >12331*10^7 :         w0_s  = w0 (i, lay_index)
        1072 >12331*10^7 :         g_s   = g  (i, lay_index)
        1073 >12331*10^7 :         mu0_s = mu0(i, lay_index)
        1074             :         !
        1075             :         ! Zdunkowski Practical Improved Flux Method "PIFM"
        1076             :         !  (Zdunkowski et al., 1980;  Contributions to Atmospheric Physics 53, 147-66)
        1077             :         !
        1078 >12331*10^7 :         gamma1 = (8._wp - w0_s * (5._wp + 3._wp * g_s)) * .25_wp
        1079 >12331*10^7 :         gamma2 =  3._wp *(w0_s * (1._wp -         g_s)) * .25_wp
        1080             :         !
        1081             :         ! Direct reflect and transmission
        1082             :         !
        1083             :         ! Eq 18;  k = SQRT(gamma1**2 - gamma2**2), limited below to avoid div by 0.
        1084             :         !   k = 0 for isotropic, conservative scattering; this lower limit on k
        1085             :         !   gives relative error with respect to conservative solution
        1086             :         !   of < 0.1% in Rdif down to tau = 10^-9
        1087 >12331*10^7 :         k = sqrt(max((gamma1 - gamma2) * (gamma1 + gamma2), min_k))
        1088 >12331*10^7 :         exp_minusktau = exp(-tau_s*k)
        1089 >12331*10^7 :         exp_minus2ktau = exp_minusktau * exp_minusktau
        1090             :
        1091             :         ! Refactored to avoid rounding errors when k, gamma1 are of very different magnitudes
        1092             :         RT_term = 1._wp / (k      * (1._wp + exp_minus2ktau)  + &
        1093 >12331*10^7 :                            gamma1 * (1._wp - exp_minus2ktau) )
        1094             :         ! Equation 25
        1095 >12331*10^7 :         Rdif(i,lay_index) = RT_term * gamma2 * (1._wp - exp_minus2ktau)
        1096             :
        1097             :         ! Equation 26
        1098 >12331*10^7 :         Tdif(i,lay_index) = RT_term * 2._wp * k * exp_minusktau
        1099             :
        1100             :         !
        1101             :         ! On a round earth, where mu0 can increase with depth in the atmosphere,
        1102             :         !   levels with mu0 <= 0 have no direct beam and hence no source for diffuse light
        1103             :         !
        1104 >13150*10^7 :         if(mu0_s > 0._wp) then
        1105 >12331*10^7 :           k_mu     = k * mu0_s
        1106             :           !
        1107             :           ! Equation 14, multiplying top and bottom by exp(-k*tau)
        1108             :           !   and rearranging to avoid div by 0.
        1109             :           !
        1110             :           RT_term =  w0_s * RT_term/merge(1._wp - k_mu*k_mu, &
        1111             :                                           epsilon(1._wp),    &
        1112 >12331*10^7 :                                           abs(1._wp - k_mu*k_mu) >= epsilon(1._wp))
        1113             :           !
        1114             :           ! Zdunkowski Practical Improved Flux Method "PIFM"
        1115             :           !  (Zdunkowski et al., 1980;  Contributions to Atmospheric Physics 53, 147-66)
        1116             :           !
        1117 >12331*10^7 :           gamma3 = (2._wp - 3._wp * mu0_s *         g_s ) * .25_wp
        1118 >12331*10^7 :           gamma4 =  1._wp - gamma3
        1119 >12331*10^7 :           alpha1 = gamma1 * gamma4 + gamma2 * gamma3           ! Eq. 16
        1120 >12331*10^7 :           alpha2 = gamma1 * gamma3 + gamma2 * gamma4           ! Eq. 17
        1121             :
        1122             :           !
        1123             :           ! Transmittance of direct, unscattered beam.
        1124             :           !
        1125 >12331*10^7 :           k_gamma3 = k * gamma3
        1126 >12331*10^7 :           k_gamma4 = k * gamma4
        1127 >12331*10^7 :           Tnoscat = exp(-tau_s/mu0_s)
        1128             :           Rdir = RT_term  *                                            &
        1129             :               ((1._wp - k_mu) * (alpha2 + k_gamma3)                  - &
        1130             :                (1._wp + k_mu) * (alpha2 - k_gamma3) * exp_minus2ktau - &
        1131 >12331*10^7 :                2.0_wp * (k_gamma3 - alpha2 * k_mu)  * exp_minusktau * Tnoscat)
        1132             :           !
        1133             :           ! Equation 15, multiplying top and bottom by exp(-k*tau),
        1134             :           !   multiplying through by exp(-tau/mu0) to
        1135             :           !   prefer underflow to overflow
        1136             :           ! Omitting direct transmittance
        1137             :           !
        1138             :           Tdir = -RT_term *                                                             &
        1139             :                 ((1._wp + k_mu) * (alpha1 + k_gamma4)                  * Tnoscat - &
        1140             :                  (1._wp - k_mu) * (alpha1 - k_gamma4) * exp_minus2ktau * Tnoscat - &
        1141 >12331*10^7 :                  2.0_wp * (k_gamma4 + alpha1 * k_mu)  * exp_minusktau)
        1142             :           ! Final check that energy is not spuriously created, by recognizing that
        1143             :           ! the beam can either be reflected, penetrate unscattered to the base of a layer,
        1144             :           ! or penetrate through but be scattered on the way - the rest is absorbed
        1145             :           ! Makes the equations safer in single precision. Credit: Robin Hogan, Peter Ukkonen
        1146 >12331*10^7 :           Rdir    = max(0.0_wp, min(Rdir, (1.0_wp - Tnoscat       ) ))
        1147 >12331*10^7 :           Tdir    = max(0.0_wp, min(Tdir, (1.0_wp - Tnoscat - Rdir) ))
        1148             :
        1149 >12331*10^7 :           source_up(i,lay_index) =    Rdir * dir_flux_inc(i)
        1150 >12331*10^7 :           source_dn(i,lay_index) =    Tdir * dir_flux_inc(i)
        1151 >12331*10^7 :           dir_flux_trans(i)      = Tnoscat * dir_flux_inc(i)
        1152             :         else
        1153           0 :           source_up(i,lay_index) = 0._wp
        1154           0 :           source_dn(i,lay_index) = 0._wp
        1155           0 :           dir_flux_trans(i)      = 0._wp
        1156             :         end if
        1157             :       end do
        1158             :     end do
        1159  1399006112 :     source_sfc(:) = dir_flux_trans(:)*sfc_albedo(:)
        1160             :
        1161    87194912 :   end subroutine sw_dif_and_source
        1162             : ! ---------------------------------------------------------------
        1163             : !
        1164             : ! Transport of diffuse radiation through a vertically layered atmosphere.
        1165             : !   Equations are after Shonk and Hogan 2008, doi:10.1175/2007JCLI1940.1 (SH08)
        1166             : !   This routine is shared by longwave and shortwave
        1167             : !
        1168             : ! -------------------------------------------------------------------------------------------------
        1169    87194912 : subroutine adding(ncol, nlay, top_at_1, &
        1170    87194912 :                   albedo_sfc,           &
        1171    87194912 :                   rdif, tdif,           &
        1172    87194912 :                   src_dn, src_up, src_sfc, &
        1173    87194912 :                   flux_up, flux_dn)
        1174             :   integer,                          intent(in   ) :: ncol, nlay
        1175             :   logical(wl),                      intent(in   ) :: top_at_1
        1176             :   real(wp), dimension(ncol       ), intent(in   ) :: albedo_sfc
        1177             :   real(wp), dimension(ncol,nlay  ), intent(in   ) :: rdif, tdif
        1178             :   real(wp), dimension(ncol,nlay  ), intent(in   ) :: src_dn, src_up
        1179             :   real(wp), dimension(ncol       ), intent(in   ) :: src_sfc
        1180             :   real(wp), dimension(ncol,nlay+1), intent(  out) :: flux_up
        1181             :   ! intent(inout) because top layer includes incident flux
        1182             :   real(wp), dimension(ncol,nlay+1), intent(inout) :: flux_dn
        1183             :   ! ------------------
        1184             :   integer :: ilev
        1185   174389824 :   real(wp), dimension(ncol,nlay+1)  :: albedo, &  ! reflectivity to diffuse radiation below this level
        1186             :                                                   ! alpha in SH08
        1187   174389824 :                                        src        ! source of diffuse upwelling radiation from emission or
        1188             :                                                   ! scattering of direct beam
        1189             :                                                   ! G in SH08
        1190   174389824 :   real(wp), dimension(ncol,nlay  )  :: denom      ! beta in SH08
        1191             :   ! ------------------
        1192             :   !
        1193             :   ! Indexing into arrays for upward and downward propagation depends on the vertical
        1194             :   !   orientation of the arrays (whether the domain top is at the first or last index)
        1195             :   ! We write the loops out explicitly so compilers will have no trouble optimizing them.
        1196             :   !
        1197    87194912 :   if(top_at_1) then
        1198    87194912 :     ilev = nlay + 1
        1199             :     ! Albedo of lowest level is the surface albedo...
        1200  1399006112 :     albedo(:,ilev)  = albedo_sfc(:)
        1201             :     ! ... and source of diffuse radiation is surface emission
        1202  1399006112 :     src(:,ilev) = src_sfc(:)
        1203             :
        1204             :     !
        1205             :     ! From bottom to top of atmosphere --
        1206             :     !   compute albedo and source of upward radiation
        1207             :     !
        1208  8283516640 :     do ilev = nlay, 1, -1
        1209 >13150*10^7 :       denom(:, ilev) = 1._wp/(1._wp - rdif(:,ilev)*albedo(:,ilev+1))                 ! Eq 10
        1210  8196321728 :       albedo(:,ilev) = rdif(:,ilev) + &
        1211 >13150*10^7 :                        tdif(:,ilev)*tdif(:,ilev) * albedo(:,ilev+1) * denom(:,ilev) ! Equation 9
        1212             :       !
        1213             :       ! Equation 11 -- source is emitted upward radiation at top of layer plus
        1214             :       !   radiation emitted at bottom of layer,
        1215             :       !   transmitted through the layer and reflected from layers below (tdiff*src*albedo)
        1216             :       !
        1217             :       src(:,ilev) =  src_up(:, ilev) + &
        1218             :                      tdif(:,ilev) * denom(:,ilev) *       &
        1219 >13159*10^7 :                        (src(:,ilev+1) + albedo(:,ilev+1)*src_dn(:,ilev))
        1220             :     end do
        1221             :
        1222             :     ! Eq 12, at the top of the domain upwelling diffuse is due to ...
        1223  1399006112 :     ilev = 1
        1224             :     flux_up(:,ilev) = flux_dn(:,ilev) * albedo(:,ilev) + & ! ... reflection of incident diffuse and
        1225  1399006112 :                       src(:,ilev)                          ! emission from below
        1226             :
        1227             :     !
        1228             :     ! From the top of the atmosphere downward -- compute fluxes
        1229             :     !
        1230  8283516640 :     do ilev = 2, nlay+1
        1231  8196321728 :       flux_dn(:,ilev) = (tdif(:,ilev-1)*flux_dn(:,ilev-1) + &  ! Equation 13
        1232  8196321728 :                          rdif(:,ilev-1)*src(:,ilev) +       &
        1233 >14789*10^7 :                          src_dn(:,ilev-1)) * denom(:,ilev-1)
        1234             :       flux_up(:,ilev) = flux_dn(:,ilev) * albedo(:,ilev) + & ! Equation 12
        1235 >13159*10^7 :                         src(:,ilev)
        1236             :     end do
        1237             :   else
        1238           0 :     ilev = 1
        1239             :     ! Albedo of lowest level is the surface albedo...
        1240           0 :     albedo(:,ilev)  = albedo_sfc(:)
        1241             :     ! ... and source of diffuse radiation is surface emission
        1242           0 :     src(:,ilev) = src_sfc(:)
        1243             :
        1244             :     !
        1245             :     ! From bottom to top of atmosphere --
        1246             :     !   compute albedo and source of upward radiation
        1247             :     !
        1248           0 :     do ilev = 1, nlay
        1249           0 :       denom(:, ilev  ) = 1._wp/(1._wp - rdif(:,ilev)*albedo(:,ilev))                ! Eq 10
        1250           0 :       albedo(:,ilev+1) = rdif(:,ilev) + &
        1251           0 :                          tdif(:,ilev)*tdif(:,ilev) * albedo(:,ilev) * denom(:,ilev) ! Equation 9
        1252             :       !
        1253             :       ! Equation 11 -- source is emitted upward radiation at top of layer plus
        1254             :       !   radiation emitted at bottom of layer,
        1255             :       !   transmitted through the layer and reflected from layers below (tdiff*src*albedo)
        1256             :       !
        1257             :       src(:,ilev+1) =  src_up(:, ilev) +  &
        1258             :                        tdif(:,ilev) * denom(:,ilev) *       &
        1259           0 :                        (src(:,ilev) + albedo(:,ilev)*src_dn(:,ilev))
        1260             :     end do
        1261             :
        1262             :     ! Eq 12, at the top of the domain upwelling diffuse is due to ...
        1263           0 :     ilev = nlay+1
        1264             :     flux_up(:,ilev) = flux_dn(:,ilev) * albedo(:,ilev) + & ! ... reflection of incident diffuse and
        1265           0 :                       src(:,ilev)                          ! scattering by the direct beam below
        1266             :
        1267             :     !
        1268             :     ! From the top of the atmosphere downward -- compute fluxes
        1269             :     !
        1270           0 :     do ilev = nlay, 1, -1
        1271           0 :       flux_dn(:,ilev) = (tdif(:,ilev)*flux_dn(:,ilev+1) + &  ! Equation 13
        1272           0 :                          rdif(:,ilev)*src(:,ilev) + &
        1273           0 :                          src_dn(:, ilev)) * denom(:,ilev)
        1274             :       flux_up(:,ilev) = flux_dn(:,ilev) * albedo(:,ilev) + & ! Equation 12
        1275           0 :                         src(:,ilev)
        1276             :
        1277             :     end do
        1278             :   end if
        1279    87194912 : end subroutine adding
        1280             : end module mo_rte_solver_kernels

Generated by: LCOV version 1.14