source: LMDZ6/trunk/libf/phylmd/ecrad/radiation/radiation_spartacus_lw.F90 @ 5130

! radiation_spartacus_lw.F90 - SPARTACUS longwave solver
!
! (C) Copyright 2014- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
!
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.
!
! Author:  Robin Hogan
! Email:   r.j.hogan@ecmwf.int
!
! Modifications
!   2017-04-11  R. Hogan  Receive emission/albedo rather than planck/emissivity
!   2017-04-22  R. Hogan  Store surface fluxes at all g-points
!   2018-09-03  R. Hogan  Security via min_cloud_effective_size
!   2018-10-08  R. Hogan  Call calc_region_properties
!   2019-01-12  R. Hogan  Use inv_inhom_effective_size if allocated

module radiation_spartacus_lw

  public

contains

  ! Small routine for scaling cloud optical depth in the cloudy
  ! regions
#include "radiation_optical_depth_scaling.h"

  ! This module contains just one subroutine, the longwave solver
  ! using the Speedy Algorithm for Radiative Transfer through Cloud
  ! Sides (SPARTACUS), which can represent 3D effects using a matrix
  ! form of the two-stream equations.
  !
  ! Sections:
  !   1: Prepare general variables and arrays
  !   2: Prepare column-specific variables and arrays
  !   3: First loop over layers
  !     3.1: Layer-specific initialization
  !     3.2: Compute gamma variables
  !       3.2a: Clear-sky case
  !       3.2b: Cloudy case
  !     3.3: Compute reflection, transmission and emission
  !       3.3a: g-points with 3D effects
  !       3.3b: g-points without 3D effects
  !   4: Compute total sources and albedos
  !   5: Compute fluxes
  subroutine solver_spartacus_lw(nlev,istartcol,iendcol, &
       &  config, thermodynamics, cloud, & 
       &  od, ssa, g, od_cloud, ssa_cloud, g_cloud, planck_hl, &
       &  emission, albedo, &
       &  flux)

    use parkind1,                 only : jprb
    use yomhook,                  only : lhook, dr_hook, jphook

    use radiation_io,             only : nulout
    use radiation_config,         only : config_type, IPdfShapeGamma
    use radiation_thermodynamics, only : thermodynamics_type
    use radiation_cloud,          only : cloud_type
    use radiation_regions,        only : calc_region_properties
    use radiation_overlap,        only : calc_overlap_matrices
    use radiation_flux,           only : flux_type, indexed_sum
    use radiation_matrix
    use radiation_two_stream,     only : calc_two_stream_gammas_lw, &
         calc_reflectance_transmittance_lw, LwDiffusivityWP
    use radiation_lw_derivatives, only : calc_lw_derivatives_matrix
    use radiation_constants,      only : Pi, GasConstantDryAir, &
         AccelDueToGravity

    implicit none

    ! Inputs
    integer, intent(in) :: nlev               ! number of model levels
    integer, intent(in) :: istartcol, iendcol ! range of columns to process
    type(config_type), intent(in)        :: config
    type(thermodynamics_type),intent(in) :: thermodynamics
    type(cloud_type),   intent(in)       :: cloud

    ! Gas and aerosol optical depth of each layer at each longwave
    ! g-point
    real(jprb), intent(in), dimension(config%n_g_lw,nlev,istartcol:iendcol) :: od

    ! Gas and aerosol single-scattering albedo and asymmetry factor,
    ! only if longwave scattering by aerosols is to be represented
    real(jprb), intent(in), &
         &  dimension(config%n_g_lw_if_scattering,nlev,istartcol:iendcol) :: ssa, g

    ! Cloud and precipitation optical depth of each layer in each
    ! longwave band
    real(jprb), intent(in) :: od_cloud(config%n_bands_lw,nlev,istartcol:iendcol)

    ! Cloud and precipitation single-scattering albedo and asymmetry
    ! factor, only if longwave scattering by clouds is to be
    ! represented
    real(jprb), intent(in), dimension(config%n_bands_lw_if_scattering, &
         &                            nlev,istartcol:iendcol) :: ssa_cloud, g_cloud

    ! Planck function (emitted flux from a black body) at half levels
    real(jprb), intent(in), dimension(config%n_g_lw,nlev+1,istartcol:iendcol) &
         &  :: planck_hl

    ! Emission (Planck*emissivity) and albedo (1-emissivity) at the
    ! surface at each longwave g-point
    real(jprb), intent(in), dimension(config%n_g_lw, istartcol:iendcol) &
         &  :: emission, albedo

    ! Output
    type(flux_type),          intent(inout):: flux

    integer :: nreg, ng
    integer :: nregActive ! =1 in clear layer, =nreg in a cloudy layer
    integer :: jcol, jlev, jg, jreg, iband, jreg2
    integer :: ng3D ! Number of g-points with small enough gas optical
                    ! depth that 3D effects need to be represented

    ! Ratio of gas constant for dry air to acceleration due to gravity
    real(jprb), parameter &
         &  :: R_over_g = GasConstantDryAir / AccelDueToGravity

    ! Used in computing rates of lateral radiation transfer
    real(jprb), parameter :: four_over_pi = 4.0_jprb / Pi

    ! The tangent of the effective zenith angle for diffuse radiation
    ! that is appropriate for 3D transport
    real(jprb), parameter :: tan_diffuse_angle_3d = Pi * 0.5_jprb
    ! Incorrect value of tand(53) used by Hogan & Shonk (2013)
    !    real(jprb), parameter :: tan_diffuse_angle_3d = 1.32704482162041

    ! Optical depth, single scattering albedo and asymmetry factor in
    ! each region at each g-point
    real(jprb), dimension(config%n_g_lw, config%nregions) &
         &  :: od_region, ssa_region, g_region

    ! Scattering optical depths of gases and clouds
    real(jprb) :: scat_od, scat_od_cloud

    ! The area fractions of each region
    real(jprb) :: region_fracs(1:config%nregions,nlev,istartcol:iendcol)

    ! The scaling used for the optical depth in the cloudy regions
    real(jprb) :: od_scaling(2:config%nregions,nlev,istartcol:iendcol)

    ! The length of the interface between regions per unit area of
    ! gridbox, equal to L_diff^ab in Hogan and Shonk (2013). This is
    ! actually the effective length oriented to a photon with random
    ! azimuth angle. The three values are the edge length between
    ! regions a-b, b-c and a-c.
    real(jprb) :: edge_length(3)

    ! Element i,j gives the rate of 3D transfer of diffuse
    ! radiation from region i to region j, multiplied by the thickness
    ! of the layer in m
    real(jprb) :: transfer_rate(config%nregions,config%nregions)

    ! Directional overlap matrices defined at all layer interfaces
    ! including top-of-atmosphere and the surface
    real(jprb), dimension(config%nregions,config%nregions,nlev+1,istartcol:iendcol) &
         &  :: u_matrix, v_matrix

    ! Two-stream variables
    real(jprb), dimension(config%n_g_lw, config%nregions) &
         &  :: gamma1, gamma2

    ! Matrix Gamma multiplied by the layer thickness z1, so units
    ! metres.  After calling expm, this array contains the matrix
    ! exponential of the original.
    real(jprb), dimension(config%n_g_lw, 2*config%nregions, &
         &  2*config%nregions) :: Gamma_z1

    ! Diffuse reflection and transmission matrices of each layer
    real(jprb), dimension(config%n_g_lw, config%nregions, &
         &  config%nregions, nlev) :: reflectance, transmittance

    ! Clear-sky diffuse reflection and transmission matrices of each
    ! layer
    real(jprb), dimension(config%n_g_lw, nlev) :: ref_clear, trans_clear

    ! If the Planck term at the top of a layer in each region is the
    ! vector b0 then the following is vector [-b0; b0]*dz
    real(jprb), dimension(config%n_g_lw,2*config%nregions) :: planck_top

    ! The difference between the Planck terms at the bottom and top of
    ! a layer, doubled as with planck_top; in terms of b' in the
    ! paper, planck_diff is [-b'; b']*dz*dz
    real(jprb), dimension(config%n_g_lw,2*config%nregions) :: planck_diff

    ! Parts of the particular solution associated with the
    ! inhomogeneous ODE. In terms of quantities from the paper,
    ! solution0 is [c0;d0] and solution_diff is [c';d']*dz.
    real(jprb), dimension(config%n_g_lw,2*config%nregions) &
         &  :: solution0, solution_diff

    ! Used for computing the Planck emission per layer
    real(jprb), dimension(config%n_g_lw,config%nregions) &
         &  :: tmp_vectors

    ! The fluxes upwelling from the top of a layer (source_up) and
    ! downwelling from the bottom of the layer (source_dn) due to
    ! emission
    real(jprb), dimension(config%n_g_lw, config%nregions, nlev) &
         &  :: source_up, source_dn
    ! ...clear-sky equivalents
    real(jprb), dimension(config%n_g_lw, nlev) &
         &  :: source_up_clear, source_dn_clear

    ! Upwelling radiation just above a layer interface due to emission
    ! below that interface, where level index = 1 corresponds to the
    ! top-of-atmosphere
    real(jprb), dimension(config%n_g_lw, config%nregions, nlev+1) &
         &  :: total_source
    ! ...clear-sky equivalent
    real(jprb), dimension(config%n_g_lw, nlev+1) :: total_source_clear

    ! As total_source, but just below a layer interface
    real(jprb), dimension(config%n_g_lw, config%nregions) &
         &  :: total_source_below

    ! Total albedo of the atmosphere/surface just above a layer
    ! interface with respect to downwelling diffuse radiation at that
    ! interface, where level index = 1 corresponds to the
    ! top-of-atmosphere
    real(jprb), dimension(config%n_g_lw, config%nregions, &
         &  config%nregions, nlev+1) :: total_albedo
    ! ...clear-sky equivalent
    real(jprb), dimension(config%n_g_lw, nlev+1) :: total_albedo_clear

    ! As total_albedo, but just below a layer interface
    real(jprb), dimension(config%n_g_lw, config%nregions, config%nregions) &
         &  :: total_albedo_below

    ! The following is used to store matrices of the form I-A*B that
    ! are used on the denominator of some expressions
    real(jprb), dimension(config%n_g_lw, config%nregions, config%nregions) &
         &  :: denominator
    ! Clear-sky equivalent, but actually its reciprocal to replace
    ! some divisions by multiplications
    real(jprb), dimension(config%n_g_lw) :: inv_denom_scalar

    ! Layer depth (m)
    real(jprb) :: dz

    ! Upwelling and downwelling fluxes above/below layer interfaces
    real(jprb), dimension(config%n_g_lw, config%nregions) &
         &  :: flux_up_above, flux_dn_above, flux_dn_below
    ! Clear-sky upwelling and downwelling fluxes (which don't
    ! distinguish between whether they are above/below a layer
    ! interface)
    real(jprb), dimension(config%n_g_lw) :: flux_up_clear, flux_dn_clear

    ! Parameterization of internal flux distribution in a cloud that
    ! can lead to the flux just about to exit a cloud being different
    ! from the mean in-cloud value.  "lateral_od" is the typical
    ! absorption optical depth through the cloud in a horizontal
    ! direction.
    real(jprb), dimension(config%n_g_lw) :: lateral_od, sqrt_1_minus_ssa
    real(jprb), dimension(config%n_g_lw) :: side_emiss_thick, side_emiss

    ! In the optically thin limit, the flux near the edge is greater
    ! than that in the interior
    real(jprb), parameter :: side_emiss_thin = 1.4107

    real(jprb) :: aspect_ratio

    ! Keep a count of the number of calls to the two ways of computing
    ! reflectance/transmittance matrices
    integer :: n_calls_expm, n_calls_meador_weaver

    ! Identify clear-sky layers, with pseudo layers for outer space
    ! and below the ground, both treated as single-region clear skies
    logical :: is_clear_sky_layer(0:nlev+1)

    real(jphook) :: hook_handle

    if (lhook) call dr_hook('radiation_spartacus_lw:solver_spartacus_lw',0,hook_handle)

    ! --------------------------------------------------------
    ! Section 1: Prepare general variables and arrays
    ! --------------------------------------------------------

    ! Copy array dimensions to local variables for convenience
    nreg = config%nregions
    ng = config%n_g_lw

    ! Reset count of number of calls to the two ways to compute
    ! reflectance/transmission matrices
    n_calls_expm = 0
    n_calls_meador_weaver = 0

    ! Initialize dz to avoid compiler warning
    dz = 1.0_jprb

    ! Compute the wavelength-independent region fractions and
    ! optical-depth scalings
    call calc_region_properties(nlev,nreg,istartcol,iendcol, &
         &  config%i_cloud_pdf_shape == IPdfShapeGamma, &
         &  cloud%fraction, cloud%fractional_std, region_fracs, &
         &  od_scaling, config%cloud_fraction_threshold)

    ! Compute wavelength-independent overlap matrices u_matrix and v_matrix
    call calc_overlap_matrices(nlev, nreg, istartcol, iendcol, &
         &  region_fracs, cloud%overlap_param, &
         &  u_matrix, v_matrix, decorrelation_scaling=config%cloud_inhom_decorr_scaling, &
         &  cloud_fraction_threshold=config%cloud_fraction_threshold, &
         &  use_beta_overlap=config%use_beta_overlap, &
         &  cloud_cover=flux%cloud_cover_lw)

    if (config%iverbose >= 3) then
      write(nulout,'(a)',advance='no') '  Processing columns'
    end if

    ! Main loop over columns
    do jcol = istartcol, iendcol
      ! --------------------------------------------------------
      ! Section 2: Prepare column-specific variables and arrays
      ! --------------------------------------------------------

      if (config%iverbose >= 3) then
        write(nulout,'(a)',advance='no') '.'
      end if

      ! Define which layers contain cloud; assume that
      ! cloud%crop_cloud_fraction has already been called
      is_clear_sky_layer = .true.
      do jlev = 1,nlev
        if (cloud%fraction(jcol,jlev) > 0.0_jprb) then
          is_clear_sky_layer(jlev) = .false.
        end if
      end do

      ! --------------------------------------------------------
      ! Section 3: First loop over layers
      ! --------------------------------------------------------
      ! In this section the reflectance, transmittance and sources
      ! are computed for each layer
      do jlev = 1,nlev ! Start at top-of-atmosphere
        ! --------------------------------------------------------
        ! Section 3.1: Layer-specific initialization
        ! --------------------------------------------------------
        ! Array-wise assignments
        gamma1 = 0.0_jprb
        gamma2 = 0.0_jprb
        Gamma_z1= 0.0_jprb
        planck_top = 0.0_jprb
        planck_diff = 0.0_jprb
        transfer_rate(:,:) = 0.0_jprb
        solution0 = 0.0_jprb
        solution_diff = 0.0_jprb

        ! Reset clear-sky absorption/scattering properties in each
        ! region
        od_region = 0.0_jprb
        ssa_region = 0.0_jprb ! No Rayleigh scattering in longwave by default
        g_region = 0.0_jprb

        ! Set optical depth of clear-sky region (region 1) to the
        ! gas/aerosol optical depth
        od_region(1:ng,1) = od(1:ng, jlev, jcol)

        ! Can't deal with zero gas optical depths, but this is now
        ! dealt with in radiation_ifsrrtm
        ! where (od_region(:,1) < 1.0e-15_jprb) od_region(:,1) = 1.0e-15_jprb

        ! Set clear-sky single-scattering albedo and asymmetry
        ! factor, if non-zero
        if (config%do_lw_aerosol_scattering) then
          ssa_region(1:ng,1) = ssa(1:ng, jlev, jcol)
          g_region(1:ng,1)   = g(1:ng, jlev, jcol)
        end if

        ! --------------------------------------------------------
        ! Section 3.2: Compute gamma variables
        ! --------------------------------------------------------
        if (is_clear_sky_layer(jlev)) then
          ! --- Section 3.2a: Clear-sky case --------------------

          nregActive = 1   ! Only region 1 (clear-sky) is active

          ! Compute two-stream variables gamma1 and gamma2
          call calc_two_stream_gammas_lw(ng, &
               &  ssa_region(1:ng,1), g_region(1:ng,1), &
               &  gamma1(1:ng,1), gamma2(1:ng,1))

          if (config%use_expm_everywhere) then
            ! Treat cloud-free layer as 3D: the matrix exponential
            ! "expm" is used as much as possible (provided the
            ! optical depth is not too high). ng3D is initially
            ! set to the total number of g-points, then the
            ! clear-sky optical depths are searched until a value
            ! exceeds the threshold for treating 3D effects, and
            ! if found it is reset to that.  Note that we are
            ! assuming that the g-points have been reordered in
            ! approximate order of gas optical depth.
            ng3D = ng
            do jg = 1, ng
              if (od_region(jg,1) > config%max_gas_od_3D) then
                ng3D = jg-1
                exit
              end if
            end do
          else
            ! Otherwise treat cloud-free layer using the classical
            ! Meador & Weaver formulae for all g-points
            ng3D = 0
          end if

        else
          ! --- Section 3.2b: Cloudy case -----------------------

          ! Default number of g-points to treat with
          ! matrix-exponential scheme
          if (config%use_expm_everywhere) then
            ng3D = ng   ! All g-points
          else
            ng3D = 0    ! No g-points
          end if

          ! Do we compute 3D effects; note that if we only have one
          ! region and the sky is overcast then 3D calculations must
          ! be turned off as there will be only one region
          if (config%do_3d_effects .and. &
               &  allocated(cloud%inv_cloud_effective_size) .and. &
               &  .not. (nreg == 2 .and. cloud%fraction(jcol,jlev) &
               &  > 1.0_jprb-config%cloud_fraction_threshold)) then
            if (cloud%inv_cloud_effective_size(jcol,jlev) &
                 &  > 0.0_jprb) then
              ! 3D effects are only simulated if
              ! inv_cloud_effective_size is defined and greater
              ! than zero
              ng3D = ng

              ! The following is from the hydrostatic equation
              ! and ideal gas law: dz = dp * R * T / (p * g)
              dz = R_over_g &
                   &  * (thermodynamics%pressure_hl(jcol,jlev+1) &
                   &   - thermodynamics%pressure_hl(jcol,jlev)) &
                   &  * (thermodynamics%temperature_hl(jcol,jlev) &
                   &   + thermodynamics%temperature_hl(jcol,jlev+1)) &
                   &  / (thermodynamics%pressure_hl(jcol,jlev) &
                   &   + thermodynamics%pressure_hl(jcol,jlev+1))

              ! Compute cloud edge length per unit area of gridbox
              ! from rearranging Hogan & Shonk (2013) Eq. 45, but
              ! adding a factor of (1-frac) so that a region that
              ! fully occupies the gridbox (frac=1) has an edge of
              ! zero. We can therefore use the fraction of clear sky,
              ! region_fracs(1,jlev,jcol) for convenience instead. The
              ! pi on the denominator means that this is actually edge
              ! length with respect to a light ray with a random
              ! azimuthal direction.
              edge_length(1) = four_over_pi &
                   &  * region_fracs(1,jlev,jcol)*(1.0_jprb-region_fracs(1,jlev,jcol)) &
                   &  * min(cloud%inv_cloud_effective_size(jcol,jlev), &
                   &        1.0_jprb / config%min_cloud_effective_size)
              if (nreg > 2) then
                ! The corresponding edge length between the two cloudy
                ! regions is computed in the same way but treating the
                ! optically denser of the two regions (region 3) as
                ! the cloud; note that the fraction of this region may
                ! be less than that of the optically less dense of the
                ! two regions (region 2).  For increased flexibility,
                ! the user may specify the effective size of
                ! inhomogeneities separately from the cloud effective
                ! size.
                if (allocated(cloud%inv_inhom_effective_size)) then
                  edge_length(2) = four_over_pi &
                       &  * region_fracs(3,jlev,jcol)*(1.0_jprb-region_fracs(3,jlev,jcol)) &
                       &  * min(cloud%inv_inhom_effective_size(jcol,jlev), &
                       &        1.0_jprb / config%min_cloud_effective_size)
                else
                  edge_length(2) = four_over_pi &
                       &  * region_fracs(3,jlev,jcol)*(1.0_jprb-region_fracs(3,jlev,jcol)) &
                       &  * min(cloud%inv_cloud_effective_size(jcol,jlev), &
                       &        1.0_jprb / config%min_cloud_effective_size)
                end if
                ! In the case of three regions, some of the cloud
                ! boundary may go directly to the "thick" region
                if (config%clear_to_thick_fraction > 0.0_jprb) then
                  edge_length(3) = config%clear_to_thick_fraction &
                       &  * min(edge_length(1), edge_length(2))
                  edge_length(1) = edge_length(1) - edge_length(3)
                  edge_length(2) = edge_length(2) - edge_length(3)
                else 
                  edge_length(3) = 0.0_jprb
                end if
              end if

              do jreg = 1, nreg-1
                ! Compute lateral transfer rate from region jreg to
                ! jreg+1 following Hogan & Shonk (2013) Eq. 47, but
                ! multiplied by dz because the transfer rate is
                ! vertically integrated across the depth of the layer
                if (region_fracs(jreg,jlev,jcol) > epsilon(1.0_jprb)) then
                  transfer_rate(jreg,jreg+1) = dz &
                       &  * edge_length(jreg) &
                       &  * tan_diffuse_angle_3d / region_fracs(jreg,jlev,jcol)
                end if
                ! Compute transfer rate from region jreg+1 to jreg
                if (region_fracs(jreg+1,jlev,jcol) > epsilon(1.0_jprb)) then
                  transfer_rate(jreg+1,jreg) = dz &
                       &  * edge_length(jreg) &
                       &  * tan_diffuse_angle_3d / region_fracs(jreg+1,jlev,jcol)
                end if
              end do

              ! Compute transfer rates directly between regions 1 and
              ! 3
              if (edge_length(3) > 0.0_jprb) then
                if (region_fracs(1,jlev,jcol) > epsilon(1.0_jprb)) then
                  transfer_rate(1,3) = dz &
                       &  * edge_length(3) &
                       &  * tan_diffuse_angle_3d / region_fracs(1,jlev,jcol)
                end if
                if (region_fracs(3,jlev,jcol) > epsilon(1.0_jprb)) then
                  transfer_rate(3,1) = dz &
                       &  * edge_length(3) &
                       &  * tan_diffuse_angle_3d / region_fracs(3,jlev,jcol)
                end if
              end if

              ! Don't allow the transfer rate out of a region to be
              ! equivalent to a loss of exp(-10) through the layer
              where (transfer_rate > config%max_3d_transfer_rate) 
                transfer_rate = config%max_3d_transfer_rate
              end where
            end if ! Cloud has edge length required for 3D effects
          end if ! Include 3D effects

          ! In a cloudy layer the number of active regions equals
          ! the number of regions
          nregActive = nreg

          ! Compute scattering properties of the regions at each
          ! g-point, mapping from the cloud properties
          ! defined in each band.
          do jg = 1,ng
            ! Mapping from g-point to band
            iband = config%i_band_from_reordered_g_lw(jg)

            ! Scattering optical depth of clear-sky region
            scat_od = od_region(jg,1)*ssa_region(jg,1)

            ! Loop over cloudy regions
            do jreg = 2,nreg
              ! Add scaled cloud optical depth to clear-sky value
              od_region(jg,jreg) = od_region(jg,1) &
                   &  + od_cloud(iband,jlev,jcol)*od_scaling(jreg,jlev,jcol)
              if (config%do_lw_cloud_scattering) then
                ! Compute single-scattering albedo of gas-cloud
                ! combination
                scat_od_cloud = od_cloud(iband,jlev,jcol) &
                     &  * ssa_cloud(iband,jlev,jcol)*od_scaling(jreg,jlev,jcol)
                ssa_region(jg,jreg) = (scat_od+scat_od_cloud) &
                     &  / od_region(jg,jreg)
                if (scat_od + scat_od_cloud > 0.0_jprb) then
                  ! Compute asymmetry factor of gas-cloud
                  ! combination
                  g_region(jg,jreg) &
                       &  = (scat_od*g_region(jg,1) &
                       &  + scat_od_cloud*g_cloud(iband,jlev,jcol)) &
                       &  / (scat_od + scat_od_cloud)
                  ! else g_region is already set to zero
                end if
                ! else ssa_region is already set to zero
              end if

              ! Apply maximum cloud optical depth for stability in the
              ! 3D case
              if (od_region(jg,jreg) > config%max_cloud_od) then
                od_region(jg,jreg) = config%max_cloud_od
              end if

            end do

            ! Calculate two-stream variables gamma1 and gamma2 of
            ! all regions at once
            call calc_two_stream_gammas_lw(nreg, &
                 &  ssa_region(jg,:), g_region(jg,:), &
                 &  gamma1(jg,:), gamma2(jg,:))

            ! Loop is in order of g-points with typically
            ! increasing optical depth: if optical depth of
            ! clear-sky region exceeds a threshold then turn off
            ! 3D effects for any further g-points
            if (ng3D == ng &
                 &  .and. od_region(jg,1) > config%max_gas_od_3D) then
              ng3D = jg-1
            end if
          end do ! Loop over g points
        end if ! Cloudy level

        ! --------------------------------------------------------
        ! Section 3.3: Compute reflection, transmission and emission
        ! --------------------------------------------------------
        if (ng3D > 0) then
          ! --- Section 3.3a: g-points with 3D effects ----------

          ! 3D effects need to be represented in "ng3D" of the g
          ! points.  This is done by creating ng3D square matrices
          ! each of dimension 2*nreg by 2*nreg, computing the matrix
          ! exponential, then computing the various
          ! transmission/reflectance matrices from that.
          do jreg = 1,nregActive
            ! Write the diagonal elements of -Gamma1*z1
            Gamma_z1(1:ng3D,jreg,jreg) &
                 &  = od_region(1:ng3D,jreg)*gamma1(1:ng3D,jreg)
            ! Write the diagonal elements of +Gamma2*z1
            Gamma_z1(1:ng3D,jreg+nreg,jreg) &
                 &  = od_region(1:ng3D,jreg)*gamma2(1:ng3D,jreg)

            ! Write the vectors corresponding to the inhomogeneous
            ! parts of the matrix ODE
            planck_top(1:ng3D,nreg+jreg) = od_region(1:ng3D,jreg) &
                 &  *(1.0_jprb-ssa_region(1:ng3D,jreg))*region_fracs(jreg,jlev,jcol) &
                 &  *planck_hl(1:ng3D,jlev,jcol)*LwDiffusivityWP
            planck_top(1:ng3D,jreg) = -planck_top(1:ng3D,nreg+jreg)
            planck_diff(1:ng3D,nreg+jreg) = od_region(1:ng3D,jreg) &
                 &  * (1.0_jprb-ssa_region(1:ng3D,jreg))*region_fracs(jreg,jlev,jcol) &
                 &  * (planck_hl(1:ng3D,jlev+1,jcol) &
                 &  -planck_hl(1:ng3D,jlev,jcol))*LwDiffusivityWP
            planck_diff(1:ng3D,jreg) = -planck_diff(1:ng3D,nreg+jreg)
          end do

          if (nregActive < nreg) then
            ! To avoid NaNs in some situations we need to put
            ! non-zeros in the cloudy-sky regions even if the cloud
            ! fraction is zero
            do jreg = nregActive+1,nreg
              Gamma_z1(1:ng3D,jreg,jreg) = Gamma_z1(1:ng3D,1,1)
              Gamma_z1(1:ng3D,nreg+jreg,jreg) = Gamma_z1(1:ng3D,nreg+1,1)
            end do
          end if

          ! Parameterization for the effective emissivity of the side
          ! of the cloud
          if (config%do_lw_side_emissivity &
             & .and. region_fracs(1,jlev,jcol) > 0.0_jprb .and. region_fracs(2,jlev,jcol) > 0.0_jprb &
             & .and. config%do_3d_effects &
             & .and. cloud%inv_cloud_effective_size(jcol,jlev) > 0.0_jprb) then
            aspect_ratio = 1.0_jprb / (min(cloud%inv_cloud_effective_size(jcol,jlev), &
                 &                         1.0_jprb / config%min_cloud_effective_size) &
                 &                     * region_fracs(1,jlev,jcol) * dz)
            lateral_od(1:ng3D) = (aspect_ratio / (nreg-1.0_jprb)) &
                 &  * sum(od_region(1:ng3D,2:nreg)*(1.0_jprb-ssa_region(1:ng3D,2:nreg)),2)
            sqrt_1_minus_ssa(1:ng3D) = sqrt(1.0_jprb - ssa_region(1:ng3D,2))
            side_emiss_thick(1:ng3D) = 2.0_jprb * sqrt_1_minus_ssa(1:ng3D) / &
                 & (sqrt_1_minus_ssa(1:ng3D) &
                 &  + sqrt(1.0_jprb-ssa_region(1:ng3D,2)*g_region(1:ng3D,2)))
            side_emiss(1:ng3D) = (side_emiss_thin - side_emiss_thick(1:ng3D)) &
                 &  / (lateral_od(1:ng3D) + 1.0_jprb) & 
                 &  + side_emiss_thick(1:ng3D)
          else
            side_emiss(1:ng3D) = 1.0_jprb
          end if

          do jreg = 1,nregActive-1
            ! Add the terms assocated with 3D transport to Gamma1*z1.
            ! First the rate from region jreg to jreg+1
            Gamma_z1(1:ng3D,jreg,jreg) = Gamma_z1(1:ng3D,jreg,jreg) &
                 &  + transfer_rate(jreg,jreg+1)
            Gamma_z1(1:ng3D,jreg+1,jreg) = -transfer_rate(jreg,jreg+1)

            ! Next the rate from region jreg+1 to jreg
            if (jreg > 1) then
              ! Flow between one cloudy region and another
              Gamma_z1(1:ng3D,jreg+1,jreg+1) = Gamma_z1(1:ng3D,jreg+1,jreg+1) &
                   &  + transfer_rate(jreg+1,jreg)
              Gamma_z1(1:ng3D,jreg,jreg+1) = -transfer_rate(jreg+1,jreg)
            else
              ! Only for the lateral transfer between cloud and clear
              ! skies do we account for the effective emissivity of
              ! the side of the cloud
              Gamma_z1(1:ng3D,jreg+1,jreg+1) = Gamma_z1(1:ng3D,jreg+1,jreg+1) &
                   &  + side_emiss(1:ng3D) * transfer_rate(jreg+1,jreg)
              Gamma_z1(1:ng3D,jreg,jreg+1) = -side_emiss(1:ng3D)*transfer_rate(jreg+1,jreg)
            end if
          end do

          ! Possible flow between regions a and c
          if (edge_length(3) > 0.0_jprb) then
            Gamma_z1(1:ng3D,1,1) = Gamma_z1(1:ng3D,1,1) &
                 &  + transfer_rate(1,3)
            Gamma_z1(1:ng3D,3,1) = -transfer_rate(1,3)
            Gamma_z1(1:ng3D,3,3) = Gamma_z1(1:ng3D,3,3) &
                 &  + side_emiss(1:ng3D) * transfer_rate(3,1)
            Gamma_z1(1:ng3D,1,3) = -side_emiss(1:ng3D)*transfer_rate(3,1)
          end if

          ! Copy Gamma1*z1
          Gamma_z1(1:ng3D,nreg+1:2*nreg,nreg+1:2*nreg) &
               &  = -Gamma_z1(1:ng3D,1:nreg,1:nreg)
          ! Copy Gamma2*z1
          Gamma_z1(1:ng3D,1:nreg,nreg+1:2*nreg) &
               &  = -Gamma_z1(1:ng3D,nreg+1:2*nreg,1:nreg)

          ! Compute the parts of the particular solution
          solution_diff(1:ng3D,1:2*nreg) &
               &  = solve_vec(ng,ng3D,2*nreg,Gamma_z1,planck_diff)
          solution_diff(1:ng3D,1:2*nreg) = - solution_diff(1:ng3D,1:2*nreg) 
          solution0(1:ng3D,1:2*nreg) = solve_vec(ng,ng3D,2*nreg,Gamma_z1, &
               &  solution_diff-planck_top)

          ! Compute the matrix exponential of Gamma_z1, returning the
          ! result in-place
          call expm(ng, ng3D, 2*nreg, Gamma_z1, IMatrixPatternDense)

          ! Update count of expm calls
          n_calls_expm = n_calls_expm + ng3D

          ! Diffuse reflectance matrix
          reflectance(1:ng3D,:,:,jlev) = -solve_mat(ng,ng3D,nreg, &
               &  Gamma_z1(1:ng3D,1:nreg,1:nreg), &
               &  Gamma_z1(1:ng3D,1:nreg,nreg+1:2*nreg))
          ! Diffuse transmission matrix
          transmittance(1:ng3D,:,:,jlev) = mat_x_mat(ng,ng3D,nreg, &
               &  Gamma_z1(1:ng3D,nreg+1:2*nreg,1:nreg), &
               &  reflectance(1:ng3D,:,:,jlev)) &
               &  + Gamma_z1(1:ng3D,nreg+1:2*nreg,nreg+1:2*nreg)

          ! Upward and downward sources due to emission within the
          ! layer
          tmp_vectors(1:ng3D,:) = solution0(1:ng3D,1:nreg) &
               &  + solution_diff(1:ng3D,1:nreg) &
               &  - mat_x_vec(ng,ng3D,nreg,Gamma_z1(:,1:nreg,nreg+1:2*nreg), &
               &  solution0(:,nreg+1:2*nreg))
          source_up(1:ng3D,:,jlev) = solution0(1:ng3D,1:nreg) &
               &  - solve_vec(ng,ng3D,nreg,Gamma_z1(:,1:nreg,1:nreg), &
               &  tmp_vectors)
          tmp_vectors(1:ng3D,:) = source_up(1:ng3D,:,jlev) &
               &  - solution0(1:ng3D,1:nreg)
          source_dn(1:ng3D,:,jlev) = mat_x_vec(ng,ng3D,nreg,&
               &  Gamma_z1(:,nreg+1:2*nreg,1:nreg),tmp_vectors) &
               &  + solution0(1:ng3D,nreg+1:2*nreg) &
               &  - mat_x_vec(ng,ng3D,nreg, &
               &  Gamma_z1(:,nreg+1:2*nreg,nreg+1:2*nreg), &
               &  solution0(:,nreg+1:2*nreg)) &
               &  + solution_diff(1:ng3D,nreg+1:2*nreg)
        end if ! we are treating 3D effects for some g points

        ! --- Section 3.3b: g-points without 3D effects ----------

        ! Compute reflectance, transmittance and upward and downward
        ! sources for clear skies, using the Meador-Weaver formulas
        ! for reflectance and transmittance, and equivalent solutions
        ! to the coupled ODEs for the sources.
        call calc_reflectance_transmittance_lw(ng, &
             &  od_region(1:ng,1), &
             &  gamma1(1:ng,1), gamma2(1:ng,1), &
             &  planck_hl(1:ng,jlev,jcol), planck_hl(1:ng,jlev+1,jcol), &
             &  ref_clear(1:ng,jlev), trans_clear(1:ng,jlev), &
             &  source_up_clear(1:ng,jlev), source_dn_clear(1:ng,jlev))

        n_calls_meador_weaver = n_calls_meador_weaver + ng

        if (ng3D < ng) then
          ! Some of the g points are to be treated using the
          ! conventional plane-parallel method.  First zero the
          ! relevant parts of the arrays
          reflectance  (ng3D+1:ng,:,:,jlev) = 0.0_jprb
          transmittance(ng3D+1:ng,:,:,jlev) = 0.0_jprb
          source_up    (ng3D+1:ng,:,jlev) = 0.0_jprb
          source_dn    (ng3D+1:ng,:,jlev) = 0.0_jprb

          ! Since there is no lateral transport, the clear-sky parts
          ! of the arrays can be copied from the clear-sky arrays
          reflectance(ng3D+1:ng,1,1,jlev) = ref_clear(ng3D+1:ng,jlev)
          transmittance(ng3D+1:ng,1,1,jlev) = trans_clear(ng3D+1:ng,jlev)
          source_up(ng3D+1:ng,1,jlev) = region_fracs(1,jlev,jcol)*source_up_clear(ng3D+1:ng,jlev)
          source_dn(ng3D+1:ng,1,jlev) = region_fracs(1,jlev,jcol)*source_dn_clear(ng3D+1:ng,jlev)

          ! Compute reflectance, transmittance and upward and downward
          ! sources, for each cloudy region, using the Meador-Weaver
          ! formulas for reflectance and transmittance, and equivalent
          ! solutions to the coupled ODEs for the sources.
          do jreg = 2, nregActive
            call calc_reflectance_transmittance_lw(ng-ng3D, &
                 &  od_region(ng3D+1:ng,jreg), &
                 &  gamma1(ng3D+1:ng,jreg), gamma2(ng3D+1:ng,jreg), &
                 &  region_fracs(jreg,jlev,jcol)*planck_hl(ng3D+1:ng,jlev,jcol), &
                 &  region_fracs(jreg,jlev,jcol)*planck_hl(ng3D+1:ng,jlev+1,jcol), &
                 &  reflectance(ng3D+1:ng,jreg,jreg,jlev), &
                 &  transmittance(ng3D+1:ng,jreg,jreg,jlev), &
                 &  source_up(ng3D+1:ng,jreg,jlev), &
                 &  source_dn(ng3D+1:ng,jreg,jlev))
          end do
          n_calls_meador_weaver &
               &  = n_calls_meador_weaver + (ng-ng3D)*(nregActive-1)
        end if

      end do ! Loop over levels

      ! --------------------------------------------------------
      ! Section 4: Compute total sources and albedos
      ! --------------------------------------------------------

      total_albedo(:,:,:,:) = 0.0_jprb
      total_source(:,:,:) = 0.0_jprb

      if (config%do_clear) then
        total_albedo_clear(:,:) = 0.0_jprb
        total_source_clear(:,:) = 0.0_jprb
      end if

      ! Calculate the upwelling radiation emitted by the surface, and
      ! copy the surface albedo into total_albedo 
      do jreg = 1,nreg
        do jg = 1,ng
          ! region_fracs(jreg,nlev,jcol) is the fraction of each
          ! region in the lowest model level
          total_source(jg,jreg,nlev+1) = region_fracs(jreg,nlev,jcol)*emission(jg,jcol)
          total_albedo(jg,jreg,jreg,nlev+1) = albedo(jg,jcol)
        end do
      end do
      ! Equivalent surface values for computing clear-sky fluxes 
      if (config%do_clear) then
        do jg = 1,ng
          total_source_clear(jg,nlev+1) = emission(jg,jcol)
        end do
        ! In the case of surface albedo there is no dependence on
        ! cloud fraction so we can copy the all-sky value
        total_albedo_clear(1:ng,nlev+1) = total_albedo(1:ng,1,1,nlev+1)
      end if

      ! Loop back up through the atmosphere computing the total albedo
      ! and the total upwelling due to emission below each level
      do jlev = nlev,1,-1
        if (config%do_clear) then
          ! Use adding method for clear-sky arrays; note that there
          ! is no need to consider "above" and "below" quantities
          ! since with no cloud overlap to worry about, these are
          ! the same
          inv_denom_scalar(:) = 1.0_jprb &
               &  / (1.0_jprb - total_albedo_clear(:,jlev+1)*ref_clear(:,jlev))
          total_albedo_clear(:,jlev) = ref_clear(:,jlev) &
               &  + trans_clear(:,jlev)*trans_clear(:,jlev)*total_albedo_clear(:,jlev+1) &
               &  * inv_denom_scalar(:)
          total_source_clear(:,jlev) = source_up_clear(:,jlev) &
               &  + trans_clear(:,jlev)*(total_source_clear(:,jlev+1) &
               &  + total_albedo_clear(:,jlev+1)*source_dn_clear(:,jlev)) &
               &  * inv_denom_scalar(:)
        end if

        if (is_clear_sky_layer(jlev)) then
          ! Clear-sky layer: use scalar adding method
          inv_denom_scalar(:) = 1.0_jprb &
               &  / (1.0_jprb - total_albedo(:,1,1,jlev+1)*reflectance(:,1,1,jlev))
          total_albedo_below = 0.0_jprb
          total_albedo_below(:,1,1) = reflectance(:,1,1,jlev) &
               &  + transmittance(:,1,1,jlev)  * transmittance(:,1,1,jlev) &
               &  * total_albedo(:,1,1,jlev+1) * inv_denom_scalar(:)
          total_source_below = 0.0_jprb
          total_source_below(:,1) = source_up(:,1,jlev) &
               &  + transmittance(:,1,1,jlev)*(total_source(:,1,jlev+1) &
               &  + total_albedo(:,1,1,jlev+1)*source_dn(:,1,jlev)) &
               &  * inv_denom_scalar(:)
        else if (config%do_3d_effects .or. &
             &   config%do_3d_lw_multilayer_effects) then
          ! Cloudy layer: use matrix adding method
          denominator = identity_minus_mat_x_mat(ng,ng,nreg,&
               &  total_albedo(:,:,:,jlev+1), reflectance(:,:,:,jlev))
          total_albedo_below = reflectance(:,:,:,jlev) &
               &  + mat_x_mat(ng,ng,nreg,transmittance(:,:,:,jlev), &
               &  solve_mat(ng,ng,nreg,denominator, &
               &  mat_x_mat(ng,ng,nreg,total_albedo(:,:,:,jlev+1), &
               &  transmittance(:,:,:,jlev))))
          total_source_below = source_up(:,:,jlev) & 
               &  + mat_x_vec(ng,ng,nreg,transmittance(:,:,:,jlev), &
               &  solve_vec(ng,ng,nreg,denominator, &
               &  total_source(:,:,jlev+1) + mat_x_vec(ng,ng,nreg, &
               &  total_albedo(:,:,:,jlev+1),source_dn(:,:,jlev))))
        else 
          ! Cloudy layer for which reflectance, transmittance and
          ! total_albedo matrices are diagonal
          total_albedo_below = 0.0_jprb
          total_source_below = 0.0_jprb
          do jreg = 1,nreg
            inv_denom_scalar(:) = 1.0_jprb / (1.0_jprb &
                 &  - total_albedo(:,jreg,jreg,jlev+1)*reflectance(:,jreg,jreg,jlev))
            total_albedo_below(:,jreg,jreg) = reflectance(:,jreg,jreg,jlev) &
                 &  + transmittance(:,jreg,jreg,jlev)*transmittance(:,jreg,jreg,jlev) &
                 &  * total_albedo(:,jreg,jreg,jlev+1) &
                 &  * inv_denom_scalar(:)
            total_source_below(:,jreg) = source_up(:,jreg,jlev) &
                 &  + transmittance(:,jreg,jreg,jlev)*(total_source(:,jreg,jlev+1) &
                 &  + total_albedo(:,jreg,jreg,jlev+1)*source_dn(:,jreg,jlev)) &
                 &  * inv_denom_scalar(:)
          end do

        end if

        ! Account for cloud overlap when converting albedo and
        ! source below a layer interface to the equivalent values
        ! just above
        if (is_clear_sky_layer(jlev) .and. is_clear_sky_layer(jlev-1)) then
          ! If both layers are cloud free, this is trivial...
          total_albedo(:,:,:,jlev) = 0.0_jprb
          total_albedo(:,1,1,jlev) = total_albedo_below(:,1,1)
          total_source(:,:,jlev) = 0.0_jprb
          total_source(:,1,jlev) = total_source_below(:,1)
        else 
          total_source(:,:,jlev) = singlemat_x_vec(ng,ng,nreg,&
               &  u_matrix(:,:,jlev,jcol), total_source_below)

!          if (config%do_3d_effects .or. config%do_3d_lw_multilayer_effects) then
          if (config%do_3d_lw_multilayer_effects) then
            ! Use the overlap matrices u_matrix and v_matrix
            total_albedo(:,:,:,jlev) = singlemat_x_mat(ng,ng,nreg,&
                 &  u_matrix(:,:,jlev,jcol), &
                 &  mat_x_singlemat(ng,ng,nreg,total_albedo_below,&
                 &  v_matrix(:,:,jlev,jcol)))
          else
            total_albedo(:,:,:,jlev) = 0.0_jprb
            ! "total_albedo" is diagonal and we wish to exclude
            ! anomalous horizontal transport described by Shonk &
            ! Hogan (2008).  Therefore, the operation we perform is
            ! essentially diag(total_albedo) = matmul(transpose(v_matrix),
            ! diag(total_albedo_below)).
            do jreg = 1,nreg
              do jreg2 = 1,nreg
                total_albedo(:,jreg,jreg,jlev) &
                     &  = total_albedo(:,jreg,jreg,jlev) &
                     &  + total_albedo_below(:,jreg2,jreg2) &
                     &  * v_matrix(jreg2,jreg,jlev,jcol)
              end do
            end do
          end if
        end if

      end do ! Reverse loop over levels

      ! --------------------------------------------------------
      ! Section 5: Compute fluxes
      ! --------------------------------------------------------

      ! Top-of-atmosphere fluxes into the regions of the top-most
      ! layer: zero since we assume no extraterrestrial longwave
      ! radiation (the shortwave scheme accounts for solar radiation
      ! in the "longwave" part of the spectrum)
      flux_dn_below = 0.0_jprb
      flux%lw_dn(jcol,1) = 0.0_jprb
      if (config%do_clear) then
        flux_dn_clear = 0.0_jprb
        flux%lw_dn_clear(jcol,1) = 0.0_jprb
      end if

      ! Store the outgoing longwave radiation at top-of-atmosphere
      flux%lw_up(jcol,1) = sum(sum(total_source(:,:,1),1))
      if (config%do_clear) then
        flux%lw_up_clear(jcol,1) = sum(total_source_clear(:,1))
      end if

      if (config%do_save_spectral_flux) then
        call indexed_sum(sum(total_source(:,:,1),2), &
             &           config%i_spec_from_reordered_g_lw, &
             &           flux%lw_up_band(:,jcol,1))
        flux%lw_dn_band(:,jcol,1) = 0.0_jprb
        if (config%do_clear) then
          call indexed_sum(total_source_clear(:,1), &
               &           config%i_spec_from_reordered_g_lw, &
               &           flux%lw_up_clear_band(:,jcol,1))
          flux%lw_dn_clear_band(:,jcol,1) = 0.0_jprb
        end if
      end if

      ! Final loop back down through the atmosphere to compute fluxes
      do jlev = 1,nlev
        if (config%do_clear) then
          ! Scalar operations for clear-sky fluxes
          flux_dn_clear(:) = (trans_clear(:,jlev)*flux_dn_clear(:) &
               + ref_clear(:,jlev)*total_source_clear(:,jlev+1) &
               + source_dn_clear(:,jlev)) &
               / (1.0_jprb - ref_clear(:,jlev)*total_albedo_clear(:,jlev+1))
          flux_up_clear(:) = total_source_clear(:,jlev+1) &
               + total_albedo_clear(:,jlev+1)*flux_dn_clear(:)
        end if

        if (is_clear_sky_layer(jlev)) then
          ! Scalar operations for clear-sky layer
          flux_dn_above(:,1) = (transmittance(:,1,1,jlev)*flux_dn_below(:,1) &
               &  + reflectance(:,1,1,jlev)*total_source(:,1,jlev+1) &
               &  + source_dn(:,1,jlev)) &
               &  / (1.0_jprb - reflectance(:,1,1,jlev)*total_albedo(:,1,1,jlev+1))
          flux_dn_above(:,2:nreg) = 0.0_jprb
          flux_up_above(:,1) = total_source(:,1,jlev+1) &
               &  + total_albedo(:,1,1,jlev+1)*flux_dn_above(:,1)
          flux_up_above(:,2:nreg) = 0.0_jprb
        else if (config%do_3d_effects .or. config%do_3d_lw_multilayer_effects) then
          ! Matrix operations for cloudy layer
          denominator = identity_minus_mat_x_mat(ng,ng,nreg,reflectance(:,:,:,jlev), &
               &  total_albedo(:,:,:,jlev+1))
          flux_dn_above = solve_vec(ng,ng,nreg,denominator, &
               &  mat_x_vec(ng,ng,nreg,transmittance(:,:,:,jlev),flux_dn_below) &
               &  + mat_x_vec(ng,ng,nreg,reflectance(:,:,:,jlev), &
               &  total_source(:,:,jlev+1)) &
               &  + source_dn(:,:,jlev))
          flux_up_above = mat_x_vec(ng,ng,nreg,total_albedo(:,:,:,jlev+1), &
               &  flux_dn_above) + total_source(:,:,jlev+1)
        else
          do jreg = 1,nreg
            ! Scalar operations for all regions, requiring that
            ! reflectance, transmittance and total_albedo are diagonal
            flux_dn_above(:,jreg) = (transmittance(:,jreg,jreg,jlev)*flux_dn_below(:,jreg) &
                 &  + reflectance(:,jreg,jreg,jlev)*total_source(:,jreg,jlev+1) &
                 &  + source_dn(:,jreg,jlev)) &
                 &  / (1.0_jprb - reflectance(:,jreg,jreg,jlev) &
                 &              * total_albedo(:,jreg,jreg,jlev+1))
            flux_up_above(:,jreg) = total_source(:,jreg,jlev+1) &
                 &  + total_albedo(:,jreg,jreg,jlev+1)*flux_dn_above(:,jreg)
          end do
        end if

        ! Account for overlap rules in translating fluxes just above
        ! a layer interface to the values just below
        if (is_clear_sky_layer(jlev) .and. is_clear_sky_layer(jlev+1)) then
          flux_dn_below = flux_dn_above
        else
          flux_dn_below = singlemat_x_vec(ng,ng,nreg,v_matrix(:,:,jlev+1,jcol), &
               &    flux_dn_above)
        end if

        ! Store the broadband fluxes
        flux%lw_up(jcol,jlev+1) = sum(sum(flux_up_above,1))
        flux%lw_dn(jcol,jlev+1) = sum(sum(flux_dn_above,1))
        if (config%do_clear) then
          flux%lw_up_clear(jcol,jlev+1) = sum(flux_up_clear)
          flux%lw_dn_clear(jcol,jlev+1) = sum(flux_dn_clear)
        end if

        if (config%do_save_spectral_flux) then
          call indexed_sum(sum(flux_up_above,2), &
               &           config%i_spec_from_reordered_g_lw, &
               &           flux%lw_up_band(:,jcol,jlev+1))
          call indexed_sum(sum(flux_dn_above,2), &
               &           config%i_spec_from_reordered_g_lw, &
               &           flux%lw_dn_band(:,jcol,jlev+1))
          if (config%do_clear) then
            call indexed_sum(flux_up_clear, &
                 &           config%i_spec_from_reordered_g_lw, &
                 &           flux%lw_up_clear_band(:,jcol,jlev+1))
            call indexed_sum(flux_dn_clear, &
                 &           config%i_spec_from_reordered_g_lw, &
                 &           flux%lw_dn_clear_band(:,jcol,jlev+1))
          end if
        end if

      end do ! Final loop over levels

      ! Store surface spectral downwelling fluxes, which at this point
      ! are at the surface
      flux%lw_dn_surf_g(:,jcol) = sum(flux_dn_above,2)
      if (config%do_clear) then
        flux%lw_dn_surf_clear_g(:,jcol) = flux_dn_clear
      end if

      ! Compute the longwave derivatives needed by Hogan and Bozzo
      ! (2015) approximate radiation update scheme
      if (config%do_lw_derivatives) then
        ! Note that at this point flux_up_above contains the spectral
        ! fluxes into the regions of the lowest layer; we sum over
        ! regions first to provide a simple spectral flux upwelling
        ! from the surface
        call calc_lw_derivatives_matrix(ng, nlev, nreg, jcol, transmittance, &
             &  u_matrix(:,:,:,jcol), sum(flux_up_above,2), flux%lw_derivatives)
      end if
        
    end do ! Loop over columns

    if (config%iverbose >= 3) then
      write(nulout,*)
    end if

    ! Report number of calls to each method of solving single-layer
    ! two-stream equations
    if (config%iverbose >= 4) then
      write(nulout,'(a,i0)') '  Matrix-exponential calls: ', n_calls_expm
      write(nulout,'(a,i0)') '  Meador-Weaver calls: ', n_calls_meador_weaver
    end if

    if (lhook) call dr_hook('radiation_spartacus_lw:solver_spartacus_lw',1,hook_handle)

  end subroutine solver_spartacus_lw

end module radiation_spartacus_lw