source: LMDZ6/trunk/libf/phylmd/ecrad/radsurf_intermediate.F90 @ 4582

! radsurf_intermediate.f90 - Derived type for intermediate radiative properties of surface
!
! (C) Copyright 2017- 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
!

module radsurf_intermediate

  use parkind1, only : jpim, jprb

  implicit none

  public

  !---------------------------------------------------------------------
  ! Derived type storing a description of radiative properties at the
  ! level of individual facets
  type surface_intermediate_type

    ! Surface "facet" properties, dimensioned (ng, nfacet, istartcol:iendcol)

    ! Longwave blackbody emission (W m-2) and emissivity from facets
    real(kind=jprb), allocatable, dimension(:,:,:) :: planck_facet
    real(kind=jprb), allocatable, dimension(:,:,:) :: lw_emissivity
    ! Shortwave direct and diffuse albedo
    real(kind=jprb), allocatable, dimension(:,:,:) :: sw_albedo_direct
    real(kind=jprb), allocatable, dimension(:,:,:) :: sw_albedo_diffuse

    ! Longwave blackbody emission (W m-2) from regions (volumes
    ! e.g. vegetation and urban canopies)
    real(kind=jprb), allocatable, dimension(:,:,:) :: planck_region
  
    ! Volumetric "region" properties of the canopy, e.g. vegetation or
    ! the space between buildings, dimensioned
    ! (ng,nregion,istartcol:iendcol)

    ! Shortwave reflectance and transmittance to diffuse radiation
    real(kind=jprb), allocatable, dimension(:,:,:) :: sw_ref_dif
    real(kind=jprb), allocatable, dimension(:,:,:) :: sw_tra_dif
    ! Shortwave reflectance and transmittance to direct radiation,
    ! where the transmittance has separate direct-to-diffuse
    ! transmittance (including scattering) and direct-to-direct
    ! transmittance (without scattering)
    real(kind=jprb), allocatable, dimension(:,:,:) :: sw_ref_dir
    real(kind=jprb), allocatable, dimension(:,:,:) :: sw_tra_dir_dif
    real(kind=jprb), allocatable, dimension(:,:,:) :: sw_tra_dir_dir
    ! Fraction of direct radiation at canyon top that is absorbed by
    ! the wall and the atmosphere
    real(kind=jprb), allocatable, dimension(:,:,:) :: sw_wall_abs_dir
    real(kind=jprb), allocatable, dimension(:,:,:) :: sw_air_abs_dir
    ! Ratio of diffuse absorption in a street canyon by the wall
    ! (rather than the air or vegetation in the canyon)
    real(kind=jprb), allocatable, dimension(:,:,:) :: sw_wall_abs_frac_dif
    real(kind=jprb), allocatable, dimension(:,:,:) :: lw_wall_abs_frac
    ! Shortwave direct and diffuse albedo at the top of a region
    real(kind=jprb), allocatable, dimension(:,:,:) :: sw_albedo_direct_reg
    real(kind=jprb), allocatable, dimension(:,:,:) :: sw_albedo_diffuse_reg
    ! Longwave reflectance, transmittance and source (note that source
    ! is the same up and down because the canopy temperature is
    ! assumed constant with height)
    real(kind=jprb), allocatable, dimension(:,:,:) :: lw_reflectance
    real(kind=jprb), allocatable, dimension(:,:,:) :: lw_transmittance
    real(kind=jprb), allocatable, dimension(:,:,:) :: lw_source
    ! Longwave emission and emissivity passed to the atmosphere scheme
    ! from the top of a region
    real(kind=jprb), allocatable, dimension(:,:,:) :: lw_emissivity_region
    real(kind=jprb), allocatable, dimension(:,:,:) :: lw_emission_region
    ! Total emission from wall and the air in an urban canopy, needed
    ! at the final partitioning stage to determine net fluxes into
    ! wall and air
    real(kind=jprb), allocatable, dimension(:,:,:) :: lw_total_wall_emission
    real(kind=jprb), allocatable, dimension(:,:,:) :: lw_total_canopy_emission

    ! Number of bands in which calculations are to be performed (can
    ! match either the spectral resolution of the atmosphere or that
    ! of the surface data)
    integer :: nswbands, nlwbands

    ! Do we represent gas radiative transfer in street/vegetation
    ! canopies?
    !logical :: do_canopy_gases_sw = .false.
    !logical :: do_canopy_gases_lw = .false.

    ! Do we use the full spectrum?
    !logical :: use_full_spectrum_sw = .false.
    !logical :: use_full_spectrum_lw = .true.

  contains

    procedure :: allocate    => allocate_surface_intermediate
    procedure :: deallocate  => deallocate_surface_intermediate
    procedure :: calc_boundary_conditions_sw
    procedure :: calc_boundary_conditions_lw
    procedure :: calc_boundary_conditions
    procedure :: partition_fluxes

  end type surface_intermediate_type

contains

  !---------------------------------------------------------------------
  subroutine allocate_surface_intermediate(this, istartcol, iendcol, &
       &                                   config, surface)
    
    use radiation_config,   only : config_type
    use radsurf_properties, only : surface_type

    class(surface_intermediate_type), intent(inout) :: this
    integer(kind=jpim),               intent(in)    :: istartcol, iendcol
    type(config_type),                intent(in)    :: config
    type(surface_type),               intent(in)    :: surface

    call this%deallocate

    !this%use_full_spectrum_sw = config%use_canopy_full_spectrum_sw
    !this%use_full_spectrum_lw = config%use_canopy_full_spectrum_lw
    ! Assume that canopy gases will only be used if we do surface
    ! calculations at the full spectral resolution
    !this%do_canopy_gases_sw   = config%do_canopy_gases_sw
    !this%do_canopy_gases_lw   = config%do_canopy_gases_lw

    if (config%use_canopy_full_spectrum_sw) then
      ! Calculations will be performed at the same spectral resolution
      ! as in the atmosphere
      this%nswbands = config%n_g_sw
    else
      ! Calculations will be performed at the same spectral resolution
      ! as the input surface-property data
      this%nswbands = surface%nalbedobands
    end if

    if (config%use_canopy_full_spectrum_lw) then
      ! Calculations will be performed at the same spectral resolution
      ! as in the atmosphere
      this%nlwbands = config%n_g_lw
    else
      ! Calculations will be performed at the same spectral resolution
      ! as the input surface-property data
      this%nlwbands = surface%nemissbands
    end if

    if (config%do_lw) then
      allocate(this%planck_facet     (this%nlwbands,surface%nfacet,istartcol:iendcol))
      allocate(this%lw_emissivity    (this%nlwbands,surface%nfacet,istartcol:iendcol))
    end if
    if (config%do_sw) then
      allocate(this%sw_albedo_direct (this%nswbands,surface%nfacet,istartcol:iendcol))
      allocate(this%sw_albedo_diffuse(this%nswbands,surface%nfacet,istartcol:iendcol))
    end if

    if (surface%nregion > 0) then
      if (config%do_lw) then
        allocate(this%planck_region   (this%nlwbands,surface%nregion,istartcol:iendcol))
        allocate(this%lw_reflectance  (this%nlwbands,surface%nregion,istartcol:iendcol))
        allocate(this%lw_transmittance(this%nlwbands,surface%nregion,istartcol:iendcol))
        allocate(this%lw_source       (this%nlwbands,surface%nregion,istartcol:iendcol))
        allocate(this%lw_emission_region(this%nlwbands,surface%nregion,istartcol:iendcol))
        allocate(this%lw_emissivity_region(this%nlwbands,surface%nregion,istartcol:iendcol))
        allocate(this%lw_wall_abs_frac(this%nlwbands,surface%nregion,istartcol:iendcol))
        allocate(this%lw_total_wall_emission(this%nlwbands,surface%nregion,istartcol:iendcol))
        allocate(this%lw_total_canopy_emission(this%nlwbands,surface%nregion,istartcol:iendcol))
      end if
      if (config%do_sw) then
        allocate(this%sw_ref_dif    (this%nswbands,surface%nregion,istartcol:iendcol))
        allocate(this%sw_tra_dif    (this%nswbands,surface%nregion,istartcol:iendcol))
        allocate(this%sw_ref_dir    (this%nswbands,surface%nregion,istartcol:iendcol))
        allocate(this%sw_tra_dir_dif(this%nswbands,surface%nregion,istartcol:iendcol))
        allocate(this%sw_tra_dir_dir(this%nswbands,surface%nregion,istartcol:iendcol))
        allocate(this%sw_wall_abs_dir(this%nswbands,surface%nregion,istartcol:iendcol))
        allocate(this%sw_air_abs_dir(this%nswbands,surface%nregion,istartcol:iendcol))
        allocate(this%sw_wall_abs_frac_dif(this%nswbands,surface%nregion,istartcol:iendcol))
        allocate(this%sw_albedo_diffuse_reg(this%nswbands,surface%nregion,istartcol:iendcol))
        allocate(this%sw_albedo_direct_reg(this%nswbands,surface%nregion,istartcol:iendcol))
      end if
    end if

  end subroutine allocate_surface_intermediate

  !---------------------------------------------------------------------
  subroutine deallocate_surface_intermediate(this)

    class(surface_intermediate_type), intent(inout) :: this

    if (allocated(this%planck_facet)) then
      deallocate(this%planck_facet)
    end if
    if (allocated(this%planck_region)) then
      deallocate(this%planck_region)
    end if
    if (allocated(this%lw_emissivity)) then
      deallocate(this%lw_emissivity)
    end if
    if (allocated(this%sw_albedo_direct)) then
      deallocate(this%sw_albedo_direct)
    end if
    if (allocated(this%sw_albedo_diffuse)) then
      deallocate(this%sw_albedo_diffuse)
    end if

    if (allocated(this%lw_reflectance)) then
      deallocate(this%lw_reflectance)
    end if
    if (allocated(this%lw_transmittance)) then
      deallocate(this%lw_transmittance)
    end if
    if (allocated(this%lw_source)) then
      deallocate(this%lw_source)
    end if
    if (allocated(this%lw_emissivity_region)) then
      deallocate(this%lw_emissivity_region)
    end if
    if (allocated(this%lw_emission_region)) then
      deallocate(this%lw_emission_region)
    end if
    if (allocated(this%lw_wall_abs_frac)) then
      deallocate(this%lw_wall_abs_frac)
    end if
    if (allocated(this%lw_total_wall_emission)) then
      deallocate(this%lw_total_wall_emission)
    end if
    if (allocated(this%lw_total_canopy_emission)) then
      deallocate(this%lw_total_canopy_emission)
    end if

    if (allocated(this%sw_ref_dif)) then
      deallocate(this%sw_ref_dif)
    end if
    if (allocated(this%sw_tra_dif)) then
      deallocate(this%sw_tra_dif)
    end if
    if (allocated(this%sw_ref_dir)) then
      deallocate(this%sw_ref_dir)
    end if
    if (allocated(this%sw_tra_dir_dif)) then
      deallocate(this%sw_tra_dir_dif)
    end if
    if (allocated(this%sw_tra_dir_dir)) then
      deallocate(this%sw_tra_dir_dir)
    end if
    if (allocated(this%sw_wall_abs_frac_dif)) then
      deallocate(this%sw_wall_abs_frac_dif)
    end if
    if (allocated(this%sw_wall_abs_dir)) then
      deallocate(this%sw_wall_abs_dir)
    end if
    if (allocated(this%sw_air_abs_dir)) then
      deallocate(this%sw_air_abs_dir)
    end if
    if (allocated(this%sw_albedo_direct_reg)) then
      deallocate(this%sw_albedo_direct_reg)
    end if
    if (allocated(this%sw_albedo_diffuse_reg)) then
      deallocate(this%sw_albedo_diffuse_reg)
    end if

  end subroutine deallocate_surface_intermediate



  !---------------------------------------------------------------------
  ! Use the detailed physical properties of the surface in "surface",
  ! and optionally the atmospheric optical properties of the lowest
  ! atmospheric level, to work out the shortwave direct/diffuse albedo
  ! that is presented to the rest of the radiation scheme, and stored
  ! in "single_level". Also, store the necessary information so that
  ! after the atmospheric radiation scheme has been run, the net
  ! fluxes at each.
  subroutine calc_boundary_conditions_sw(this, istartcol, iendcol, &   ! in
       &  config, surface, &                    ! in
       &  single_level, &                                 ! out
       &  ext_sw_air, ssa_sw_air, g_sw_air)               ! in, optional

    use yomhook,                only : lhook, dr_hook

    use radiation_io,           only : nulerr, radiation_abort
    use radiation_config,       only : config_type
    use radiation_single_level, only : single_level_type
    use radsurf_properties,     only : surface_type, ITileFlat,ITileVegetation,ITileUrban3D
    use radiation_two_stream,   only : calc_two_stream_gammas_sw, &
         &                             calc_reflectance_transmittance_sw, &
         &                             calc_reflectance_transmittance_z_sw 

    use radiation_constants,    only: Pi, GasConstantDryAir, AccelDueToGravity

    class(surface_intermediate_type), intent(inout) :: this
    integer(kind=jpim),               intent(in)    :: istartcol, iendcol
    type(config_type),                intent(in)    :: config
    type(surface_type),               intent(in)    :: surface
    type(single_level_type),          intent(inout) :: single_level

    ! Input properties of the air in the lowest model level:
    ! extinction coefficient (m-1), single-scattering albedo and
    ! asymmetry factor
    real(kind=jprb), intent(in), dimension(config%n_g_sw,istartcol:iendcol), optional &
         &  :: ext_sw_air, ssa_sw_air, g_sw_air

    ! Shortwave region properties
    real(kind=jprb), dimension(this%nswbands) :: od_sw_region, ssa_sw_region, g_sw_region
    ! Extinction coefficient better for urban areas
    real(kind=jprb), dimension(this%nswbands) :: ext_sw_region
    ! Optical depth of air
    real(kind=jprb), dimension(this%nswbands) :: od_sw_air

    real(kind=jprb)    :: tile_fraction, cos_sza

    ! Exchange coefficients (m-1) for direct and diffuse radiation
    ! into building walls
    real(kind=jprb)    :: fdiff, fdir

    ! One minus the building fraction
    real(kind=jprb)    :: canyon_fraction

    ! Tangent of solar zenith angle
    real(kind=jprb)    :: tan_sza

    ! Two-stream coefficients
    real(jprb), dimension(this%nswbands) :: gamma1_sw, gamma2_sw
    real(jprb), dimension(this%nswbands) :: gamma0_sw, gamma3_sw, gamma4_sw

    real(jprb), dimension(this%nswbands) :: inv_denominator_sw

    integer(kind=jpim) :: jcol,jtile,iregion

    ! Mapping from albedo bands to reordered shortwave g points
    integer(kind=jpim) :: i_albedo_from_g(config%n_g_sw)

    ! Indices to different facets
    integer(kind=jpim) :: iground, iwall, iroof

    real(jprb) :: hook_handle

    if (lhook) call dr_hook('radsurf_intermediate:calc_boundary_conditions_sw',0,hook_handle)

    if (present(ext_sw_air)) then
      if (.not. present(ssa_sw_air) .or. .not. present(g_sw_air)) then
        write(nulerr,'(a)') '*** Error: if ext_sw_air is provided then ssa_sw_air and g_sw_air must also be provided'
        call radiation_abort()
      end if
    end if

    if (size(single_level%sw_albedo,2) /= this%nswbands) then
      write(nulerr,'(a,i0,a,i0)') '*** Error: single-level albedo has ', &
           &  size(single_level%sw_albedo,2), ' bands, needs ', this%nswbands
      call radiation_abort()
    end if

    if (config%use_canopy_full_spectrum_sw) then
      ! Put shortwave albedo on g-point grid and permute
      i_albedo_from_g = config%i_albedo_from_band_sw(config%i_band_from_reordered_g_sw)
      do jcol = istartcol,iendcol
        this%sw_albedo_diffuse(:,:,jcol) = surface%sw_albedo(jcol,i_albedo_from_g,:)
        if (allocated(surface%sw_albedo_direct)) then
          this%sw_albedo_direct(:,:,jcol) = surface%sw_albedo_direct(jcol,i_albedo_from_g,:)
        else
          this%sw_albedo_direct = this%sw_albedo_diffuse
        end if
      end do
    else
      do jcol = istartcol,iendcol
        ! No change to spectral resolution: simply permute 
        this%sw_albedo_diffuse(:,:,jcol) = surface%sw_albedo(jcol,:,:)
        if (allocated(surface%sw_albedo_direct)) then
          this%sw_albedo_direct(:,:,jcol) = surface%sw_albedo_direct(jcol,:,:)
        else
          this%sw_albedo_direct = this%sw_albedo_diffuse
        end if
      end do
    end if

    if (surface%is_simple) then
      ! We have a "traditional" representation: one flat tile
      single_level%sw_albedo_direct(istartcol:iendcol,:) &
           &  = transpose(this%sw_albedo_direct(:,1,istartcol:iendcol))
      single_level%sw_albedo       (istartcol:iendcol,:) &
           &  = transpose(this%sw_albedo_diffuse(:,1,istartcol:iendcol))
    else
      ! More complex description of surface

      ! Firstly initialize outputs to zero
      single_level%sw_albedo_direct(istartcol:iendcol,:) = 0.0_jprb
      single_level%sw_albedo       (istartcol:iendcol,:) = 0.0_jprb

      ! Loop over column and tile
      do jcol = istartcol,iendcol
        cos_sza = single_level%cos_sza(jcol)

        do jtile = 1,surface%ntile
          tile_fraction = surface%tile_fraction(jcol,jtile)

          if (tile_fraction > 0.0_jprb) then

            select case (surface%i_representation(jtile))
            case (ITileFlat)
              ! SIMPLE FLAT TILE

              ! Add the contribution from this simple flat tile to the
              ! accumulated values for the column
              iground = surface%i_ground_facet(jcol,jtile)
              single_level%sw_albedo_direct(jcol,:) = single_level%sw_albedo_direct(jcol,:) &
                   &  + tile_fraction * this%sw_albedo_direct(:,iground,jcol)
              single_level%sw_albedo(jcol,:) = single_level%sw_albedo(jcol,:) &
                   &  + tile_fraction * this%sw_albedo_diffuse(:,iground,jcol)

            case (ITileVegetation)
              ! HORIZONTALLY HOMOGENEOUS VEGETATION CANOPY WITH
              ! SELLERS-LIKE FORMULATION

              iground = surface%i_ground_facet(jcol,jtile)
              iregion = surface%i_region_1(jcol,jtile)

              ! Shortwave calculation
              if (present(ext_sw_air)) then
                od_sw_air = surface%canopy_depth(jcol,jtile)*ext_sw_air(:,jcol)
                od_sw_region = od_sw_air + surface%vegetation_optical_depth(jcol,jtile)
                ssa_sw_region= (ssa_sw_air(:,jcol)*od_sw_air &
                     &          + surface%vegetation_optical_depth(jcol,jtile) &
                     &            * surface%vegetation_sw_albedo(jcol,:,jtile)) &
                     &       / od_sw_region
                ! Assume asymmetry factor of vegetation is zero
                g_sw_region = g_sw_air(:,jcol)*ssa_sw_air(:,jcol)*od_sw_air &
                     &  / (ssa_sw_region*od_sw_region)
              else
                od_sw_region = surface%vegetation_optical_depth(jcol,jtile)
                ssa_sw_region= surface%vegetation_sw_albedo(jcol,:,jtile)
                g_sw_region  = 0.0_jprb
              end if
              call calc_two_stream_gammas_sw(this%nswbands, &
                   &  cos_sza, ssa_sw_region, g_sw_region, &
                   &  gamma1_sw, gamma2_sw, gamma3_sw)
              call calc_reflectance_transmittance_sw(this%nswbands, &
                   &  cos_sza, od_sw_region, ssa_sw_region, &
                   &  gamma1_sw, gamma2_sw, gamma3_sw, &
                   &  this%sw_ref_dif(:,iregion,jcol), this%sw_tra_dif(:,iregion,jcol), &
                   &  this%sw_ref_dir(:,iregion,jcol), this%sw_tra_dir_dif(:,iregion,jcol), &
                   &  this%sw_tra_dir_dir(:,iregion,jcol))

              ! Shortwave adding method for a single layer
              inv_denominator_sw(:) = 1.0_jprb &
                   &  / (1.0_jprb - this%sw_albedo_diffuse(:,iground,jcol) & 
                   &               *this%sw_ref_dif(:,iregion,jcol))
              this%sw_albedo_diffuse_reg(:,iregion,jcol) = this%sw_ref_dif(:,iregion,jcol) &
                   &  + this%sw_tra_dif(:,iregion,jcol)**2 &
                   &  * this%sw_albedo_diffuse(:,iground,jcol) * inv_denominator_sw(:)
              this%sw_albedo_direct_reg(:,iregion,jcol) = this%sw_ref_dir(:,iregion,jcol) &
                   &     + (this%sw_tra_dir_dir(:,iregion,jcol)*this%sw_albedo_direct (:,iground,jcol) &
                   &       +this%sw_tra_dir_dif(:,iregion,jcol)*this%sw_albedo_diffuse(:,iground,jcol)) &
                   &     * this%sw_tra_dif(:,iregion,jcol) * inv_denominator_sw(:)
              single_level%sw_albedo(jcol,:) = single_level%sw_albedo(jcol,:) &
                   &  + tile_fraction * this%sw_albedo_diffuse_reg(:,iregion,jcol)
              single_level%sw_albedo_direct(jcol,:) = single_level%sw_albedo_direct(jcol,:) &
                   &  + tile_fraction * this%sw_albedo_direct_reg(:,iregion,jcol)

            case (ITileUrban3D)
              ! URBAN CANOPY WITH NO VEGETATION, TREATED USING
              ! SPARTACUS METHODOLOGY

              iground  = surface%i_ground_facet(jcol,jtile)
              iwall    = surface%i_wall_facet(jcol,jtile)
              iroof    = surface%i_roof_facet(jcol,jtile)
              iregion = surface%i_region_1(jcol,jtile)

              canyon_fraction = 1.0_jprb - surface%building_fraction(jcol,jtile)

              fdiff = 0.5_jprb * surface%building_normalized_perimeter(jcol,jtile) / canyon_fraction

              tan_sza = sqrt(1.0_jprb / (cos_sza*cos_sza) - 1.0_jprb)
              fdir  = surface%building_normalized_perimeter(jcol,jtile) * tan_sza / (Pi * canyon_fraction)

              if (present(ext_sw_air)) then
                ext_sw_region = ext_sw_air(:,jcol)
                ssa_sw_region = ssa_sw_air(:,jcol)
                g_sw_region   = g_sw_air(:,jcol)
              else
                ext_sw_region = 0.0_jprb
                ssa_sw_region = 0.0_jprb
                g_sw_region   = 0.0_jprb
              end if

              ! Get gammas for atmosphere only
              call calc_two_stream_gammas_sw(this%nswbands, &
                   &  cos_sza, ssa_sw_region, g_sw_region, &
                   &  gamma1_sw, gamma2_sw, gamma3_sw)

              ! At this point gamma1_sw-gamma2_sw is the rate of
              ! absorption per unit optical depth of the air in the
              ! canyon
              this%sw_wall_abs_frac_dif(:,iregion,jcol) &
                   &  = fdiff * (1.0_jprb - this%sw_albedo_diffuse(:,iwall,jcol))
              this%sw_wall_abs_frac_dif(:,iregion,jcol) = this%sw_wall_abs_frac_dif(:,iregion,jcol) &
                   &  / max(1.0e-8_jprb,ext_sw_region*(gamma1_sw-gamma2_sw) &
                   &                    + this%sw_wall_abs_frac_dif(:,iregion,jcol))

              gamma4_sw = 1.0_jprb - gamma3_sw
              gamma0_sw = ext_sw_region / cos_sza + fdir
              gamma1_sw = ext_sw_region * gamma1_sw &
                   &    + fdiff * (1.0_jprb - 0.5_jprb*this%sw_albedo_diffuse(:,iwall,jcol))
              gamma2_sw = ext_sw_region * gamma2_sw &
                   &    + fdiff * 0.5_jprb*this%sw_albedo_diffuse(:,iwall,jcol)
              gamma3_sw = ext_sw_region * ssa_sw_region * gamma3_sw &
                   &    + 0.5_jprb * fdir * this%sw_albedo_direct(:,iwall,jcol)
              gamma4_sw = ext_sw_region * ssa_sw_region * gamma4_sw &
                   &    + 0.5_jprb * fdir * this%sw_albedo_direct(:,iwall,jcol)

              call calc_reflectance_transmittance_z_sw(this%nswbands, &
                   &  cos_sza, surface%canopy_depth(jcol,jtile), gamma0_sw, &
                   &  gamma1_sw, gamma2_sw, gamma3_sw, gamma3_sw, &
                   &  this%sw_ref_dif(:,iregion,jcol), this%sw_tra_dif(:,iregion,jcol), &
                   &  this%sw_ref_dir(:,iregion,jcol), this%sw_tra_dir_dif(:,iregion,jcol), &
                   &  this%sw_tra_dir_dir(:,iregion,jcol))

              ! Compute fraction of direct at canyon top that is absorbed by wall and air
              this%sw_wall_abs_dir(:,iregion,jcol) &
                   &  = (1.0_jprb - this%sw_tra_dir_dir(:,iregion,jcol)) &
                   &  * fdir * (1.0_jprb - this%sw_albedo_direct(:,iwall,jcol)) * cos_sza &
                   &  /  max(1.0e-8_jprb, fdir*cos_sza + ext_sw_region)
              this%sw_air_abs_dir(:,iregion,jcol)&
                   &  = (1.0_jprb - this%sw_tra_dir_dir(:,iregion,jcol)) &
                   &  * ext_sw_region * (1.0_jprb - ssa_sw_region) &
                   &  /  max(1.0e-8_jprb, fdir*cos_sza + ext_sw_region)

              ! Add roof component
              single_level%sw_albedo(jcol,:) = single_level%sw_albedo(jcol,:) &
                   &  + tile_fraction * surface%building_fraction(jcol,jtile) &
                   &  * this%sw_albedo_diffuse(:,iroof,jcol)
              single_level%sw_albedo_direct(jcol,:) = single_level%sw_albedo_direct(jcol,:) &
                   &  + tile_fraction * surface%building_fraction(jcol,jtile) &
                   &  * this%sw_albedo_direct(:,iroof,jcol)

              ! Add canyon component using shortwave adding method for a single layer
              inv_denominator_sw(:) = 1.0_jprb &
                   &  / (1.0_jprb - this%sw_albedo_diffuse(:,iground,jcol) & 
                   &               *this%sw_ref_dif(:,iregion,jcol))
              this%sw_albedo_diffuse_reg(:,iregion,jcol) = this%sw_ref_dif(:,iregion,jcol) &
                   &  + this%sw_tra_dif(:,iregion,jcol)**2 &
                   &  * this%sw_albedo_diffuse(:,iground,jcol) * inv_denominator_sw(:)
              this%sw_albedo_direct_reg(:,iregion,jcol) = this%sw_ref_dir(:,iregion,jcol) &
                   &     + (this%sw_tra_dir_dir(:,iregion,jcol)*this%sw_albedo_direct (:,iground,jcol) &
                   &       +this%sw_tra_dir_dif(:,iregion,jcol)*this%sw_albedo_diffuse(:,iground,jcol)) &
                   &     * this%sw_tra_dif(:,iregion,jcol) * inv_denominator_sw(:)
              single_level%sw_albedo(jcol,:) = single_level%sw_albedo(jcol,:) &
                   &  + tile_fraction * canyon_fraction * this%sw_albedo_diffuse_reg(:,iregion,jcol)
              single_level%sw_albedo_direct(jcol,:) = single_level%sw_albedo_direct(jcol,:) &
                   &  + tile_fraction * canyon_fraction * this%sw_albedo_direct_reg(:,iregion,jcol)

            end select
          end if
        end do
      end do
    end if

    if (lhook) call dr_hook('radsurf_intermediate:calc_boundary_conditions_sw',1,hook_handle)

  end subroutine calc_boundary_conditions_sw


  !---------------------------------------------------------------------
  ! As calc_boundary_conditions_sw, but for the longwave.
  subroutine calc_boundary_conditions_lw(this, istartcol, iendcol, &   ! in
       &  config, surface, &                    ! in
       &  single_level, &                                 ! out
       &  ext_lw_air, ssa_lw_air, g_lw_air)               ! in, optional

    use yomhook,                only : lhook, dr_hook

    use radiation_io,           only : nulerr, radiation_abort
    use radiation_config,       only : config_type
    use radiation_single_level, only : single_level_type
    use radsurf_properties,     only : surface_type, ITileFlat,ITileVegetation,ITileUrban3D
    use radiation_two_stream,   only : calc_two_stream_gammas_lw, &
         &                             calc_reflectance_transmittance_isothermal_lw, &
         &                             LwDiffusivity
    use radiation_ifs_rrtm,     only : planck_function
    use radiation_constants,    only: Pi, GasConstantDryAir, AccelDueToGravity, StefanBoltzmann

    class(surface_intermediate_type), intent(inout) :: this
    integer(kind=jpim),               intent(in)    :: istartcol, iendcol
    type(config_type),                intent(in)    :: config
    type(surface_type),               intent(in)    :: surface
    type(single_level_type),          intent(inout) :: single_level

    ! Input properties of the air in the lowest model level:
    ! extinction coefficient (m-1), single-scattering albedo and
    ! asymmetry factor
    real(kind=jprb), intent(in), dimension(config%n_g_lw,istartcol:iendcol), optional &
         &  :: ext_lw_air, ssa_lw_air, g_lw_air

    ! Longwave region properties
    real(kind=jprb), dimension(this%nlwbands) :: od_lw_region, ssa_lw_region, g_lw_region

    ! Optical depth of air
    real(kind=jprb), dimension(this%nlwbands) :: od_lw_air

    ! Effective Planck function of urban canopy as a weighted average
    ! of wall and air emission
    real(kind=jprb), dimension(this%nlwbands) :: planck_canopy

    ! Vegetation emissivity using local spectral representation
    real(kind=jprb), dimension(this%nlwbands) :: vegetation_lw_emissivity

    real(kind=jprb)    :: tile_fraction

    ! Exchange coefficient (m-1) for diffuse radiation into building
    ! walls, multiplied by canyon depth (m) to get a dimensionless
    ! analogue for optical depth
    real(kind=jprb)    :: od_lw_wall

    ! One minus the building fraction
    real(kind=jprb)    :: canyon_fraction

    ! Two-stream coefficients
    real(jprb), dimension(this%nlwbands) :: gamma1_lw, gamma2_lw

    real(jprb), dimension(this%nlwbands) :: inv_denominator_lw

    integer(kind=jpim) :: jcol,jtile,jfacet,iregion

    ! Indices to different facets
    integer(kind=jpim) :: iground, iwall, iroof

    real(jprb) :: hook_handle

    if (lhook) call dr_hook('radsurf_intermediate:calc_boundary_conditions_lw',0,hook_handle)

    if (present(ssa_lw_air)) then
      if (.not. present(g_lw_air)) then
        write(nulerr,'(a)') '*** Error: if ssa_lw_air is provided then g_lw_air must also be provided'
        call radiation_abort()
      end if
    end if

    if (size(single_level%lw_emissivity,2) /= this%nlwbands) then
      write(nulerr,'(a,i0,a,i0)') '*** Error: single-level emissivity has ', &
           &  size(single_level%lw_emissivity,2), ' bands, needs ', this%nlwbands
      call radiation_abort()
    end if

    ! FIX: This section ought to check what tiles are present in each column
    if (config%use_canopy_full_spectrum_lw) then
      ! Put longwave emissivity on g-point grid and permute      
      do jcol = istartcol,iendcol
        this%lw_emissivity(:,:,jcol) = surface%lw_emissivity(jcol,config%i_emiss_from_band_lw( &
             &                           config%i_band_from_reordered_g_lw),:)
      end do
      ! Compute Planck function at each facet
      do jfacet = 1,surface%nfacet
        do jcol = istartcol,iendcol
          ! FIX this function assumes contiguous data for planck_facet
          ! so copies into a temporary
          call planck_function(config, &
               &  surface%skin_temperature(jcol,jfacet), &
               &  this%planck_facet(:,jfacet,jcol))
        end do
      end do
      do jtile = 1,surface%ntile
        do jcol = istartcol,iendcol
          iregion = surface%i_region_1(jcol,jtile)
          if (iregion > 0) then
            call planck_function(config, &
                 &  surface%canopy_temperature(jcol,jtile), &
                 &  this%planck_region(:,iregion,jcol:jcol))
          end if
        end do
      end do
    else
      if (surface%nemissbands /= 1) then
        write(nulerr,'(a)') '*** Error: insufficient information to compute Planck function in emissivity bands'
        call radiation_abort()
      end if
      do jcol = istartcol,iendcol
        ! No change to spectral resolution: simply permute
        this%lw_emissivity(:,:,jcol) = surface%lw_emissivity(jcol,:,:)
      end do
      ! Broadband calculation: Stefan-Boltzmann law
      do jfacet = 1,surface%nfacet
        this%planck_facet(1,jfacet,istartcol:iendcol) &
             &  = StefanBoltzmann*surface%skin_temperature(istartcol:iendcol,jfacet)**4
      end do
      do jtile = 1,surface%ntile
        do jcol = istartcol,iendcol
          iregion = surface%i_region_1(jcol,jtile)
          if (iregion > 0) then
            this%planck_region(1,iregion,jcol) &
                 &  = StefanBoltzmann*surface%canopy_temperature(jcol,jtile)
          end if
        end do
      end do
    end if

    if (surface%is_simple) then
      ! We have a "traditional" representation: one flat tile
      single_level%lw_emissivity(istartcol:iendcol,:) &
           &  = transpose(this%lw_emissivity(:,1,istartcol:iendcol))
      single_level%lw_emission(istartcol:iendcol,:) &
           &  = transpose(this%planck_facet(:,1,istartcol:iendcol) &
           &              * this%lw_emissivity(:,1,istartcol:iendcol))

    else
      ! More complex description of surface

      ! Firstly initialize outputs to zero
      single_level%lw_emissivity(istartcol:iendcol,:) = 0.0_jprb
      single_level%lw_emission  (istartcol:iendcol,:) = 0.0_jprb

      ! Loop over column and tile
      do jcol = istartcol,iendcol

        do jtile = 1,surface%ntile
          tile_fraction = surface%tile_fraction(jcol,jtile)

          if (tile_fraction > 0.0_jprb) then

            select case (surface%i_representation(jtile))
            case (ITileFlat)
              ! SIMPLE FLAT TILE

              ! Add the contribution from this simple flat tile to the
              ! accumulated values for the column
              iground = surface%i_ground_facet(jcol,jtile)
              single_level%lw_emissivity(jcol,:) = single_level%lw_emissivity(jcol,:) &
                   &  + tile_fraction*this%lw_emissivity(:,iground,jcol)
              single_level%lw_emission(jcol,:) = single_level%lw_emission(jcol,:) + tile_fraction &
                   &  * (this%planck_facet(:,iground,jcol)*this%lw_emissivity(:,iground,jcol))

            case (ITileVegetation)
              ! HORIZONTALLY HOMOGENEOUS VEGETATION CANOPY

              iground  = surface%i_ground_facet(jcol,jtile)
              iregion = surface%i_region_1(jcol,jtile)

              ! Convert vegetation emissivity to local spectral representation
              if (config%use_canopy_full_spectrum_lw) then
                vegetation_lw_emissivity = surface%vegetation_lw_emissivity(jcol, &
                     &  config%i_emiss_from_band_lw(config%i_band_from_reordered_g_lw), &
                     &  jtile)
              else
                vegetation_lw_emissivity = surface%vegetation_lw_emissivity(jcol,:,jtile)
              end if

              if (present(ext_lw_air)) then
                od_lw_air = surface%canopy_depth(jcol,jtile)*ext_lw_air(:,jcol)
                od_lw_region = od_lw_air + surface%vegetation_optical_depth(jcol,jtile)
                if (present(ssa_lw_air)) then
                  ssa_lw_region= (ssa_lw_air(:,jcol)*od_lw_air &
                       &          + surface%vegetation_optical_depth(jcol,jtile) &
                       &            * (1.0_jprb-vegetation_lw_emissivity)) &
                       &       / od_lw_region
                  ! Assume asymmetry factor of vegetation is zero
                  g_lw_region = g_lw_air(:,jcol)*ssa_lw_air(:,jcol)*od_lw_air &
                       &  / (ssa_lw_region * od_lw_region)
                else
                  ! No longwave scattering properties for air
                  ssa_lw_region= surface%vegetation_optical_depth(jcol,jtile) &
                       &            * (1.0_jprb-vegetation_lw_emissivity) &
                       &       / od_lw_region
                  ! Assume asymmetry factor of vegetation is zero
                  g_lw_region = 0.0_jprb

                end if
              else
                ! No gases
                od_lw_region = surface%vegetation_optical_depth(jcol,jtile)
                ssa_lw_region= 1.0_jprb-vegetation_lw_emissivity
                g_lw_region  = 0.0_jprb
              end if

              call calc_two_stream_gammas_lw(this%nlwbands, &
                   &  ssa_lw_region, g_lw_region, &
                   &  gamma1_lw, gamma2_lw)
              call calc_reflectance_transmittance_isothermal_lw(this%nlwbands, &
                   &  od_lw_region, gamma1_lw, gamma2_lw, this%planck_region(:,iregion,jcol), &
                   &  this%lw_reflectance(:,iregion,jcol), this%lw_transmittance(:,iregion,jcol), &
                   &  this%lw_source(:,iregion,jcol))
              
              ! Longwave adding method for a single layer
              inv_denominator_lw(:) = 1.0_jprb &
                   &  / (1.0_jprb - (1.0_jprb - this%lw_emissivity(:,iground,jcol)) &
                   &               *this%lw_reflectance(:,iregion,jcol))
              single_level%lw_emissivity(jcol,:) = single_level%lw_emissivity(jcol,:) &
                   &  + tile_fraction * (1.0_jprb - (this%lw_reflectance(:,iregion,jcol) &
                   &  + this%lw_transmittance(:,iregion,jcol)**2 &
                   &  * (1.0_jprb - this%lw_emissivity(:,iground,jcol)) * inv_denominator_lw(:)))
              single_level%lw_emission(jcol,:) = single_level%lw_emission(jcol,:) + tile_fraction &
                   &  * (this%lw_source(:,iregion,jcol) * (1.0_jprb + inv_denominator_lw(:) &
                   &       * (1.0_jprb - this%lw_emissivity(:,iground,jcol)) &
                   &       * this%lw_transmittance(:,iregion,jcol)) &
                   &     +this%planck_facet(:,iground,jcol)*this%lw_emissivity(:,iground,jcol) &
                   &     *this%lw_transmittance(:,iregion,jcol)*inv_denominator_lw(:))

            case (ITileUrban3D)
              ! URBAN CANOPY WITH NO VEGETATION, TREATED USING
              ! SPARTACUS METHODOLOGY

              iground  = surface%i_ground_facet(jcol,jtile)
              iwall    = surface%i_wall_facet(jcol,jtile)
              iroof    = surface%i_roof_facet(jcol,jtile)
              iregion  = surface%i_region_1(jcol,jtile)

              canyon_fraction = 1.0_jprb - surface%building_fraction(jcol,jtile)

              ! fdiff woould be 0.5 * building_normalized_perimeter /
              ! canyon_fraction, but to get equivalent zenith optical
              ! depth we multiply by the building height and divide by
              ! the longwave diffusivity angle
!              od_lw_wall = 0.5_jprb * surface%building_normalized_perimeter(jcol,jtile) &
!                   &     * surface%canopy_depth(jcol,jtile) / (LwDiffusivity * canyon_fraction)

              ! Or first compute H/W
              od_lw_wall = 0.5_jprb * surface%building_normalized_perimeter(jcol,jtile) &
                   &     * surface%canopy_depth(jcol,jtile) / canyon_fraction
              ! And then use the Harman et al. (2004) formula for
              ! street-to-sky transmittance T=sqrt[(H/W)^2+1]-H/W, to
              ! compute equivalent optical depth knowing that the
              ! two-stream scheme treatment will be T=exp(-D*od)
              od_lw_wall = -log(sqrt(od_lw_wall*od_lw_wall + 1) - od_lw_wall) / LwDiffusivity

              if (present(ext_lw_air)) then
                od_lw_air     = ext_lw_air(:,jcol)*surface%canopy_depth(jcol,jtile) !*10.0_jprb
                od_lw_region  = od_lw_air + od_lw_wall
                if (present(ssa_lw_air)) then
                  ssa_lw_region = (od_lw_air * ssa_lw_air(:,jcol) &
                       &    + od_lw_wall * (1.0_jprb - this%lw_emissivity(:,iwall,jcol))) &
                       &               / max(od_lw_region,1.0e-6_jprb)
                  ! We assume that any scattering off the walls is
                  ! equally likely to be up or down, so the effective
                  ! asymmetry factor is zero
                  g_lw_region = (od_lw_air * ssa_lw_air(:,jcol)*g_lw_air(:,jcol)) &
                       &      / max(od_lw_region*ssa_lw_region,1.0e-6_jprb)
                  ! Effective Planck function of the canopy is the
                  ! weighted average of wall and air emission
                  this%lw_total_wall_emission(:,iregion,jcol) = LwDiffusivity &
                       &  * od_lw_wall*this%lw_emissivity(:,iwall,jcol)*this%planck_facet(:,iwall,jcol)
                  this%lw_total_canopy_emission(:,iregion,jcol) = LwDiffusivity &
                       &  * od_lw_air*(1.0_jprb-ssa_lw_air(:,jcol))*this%planck_region(:,iregion,jcol)
                  planck_canopy = (this%lw_total_wall_emission(:,iregion,jcol) &
                       &           +this%lw_total_canopy_emission(:,iregion,jcol)) &
                       &        / max(od_lw_region*(1.0_jprb-ssa_lw_region)*LwDiffusivity,1.0e-6_jprb)
                else
                  ssa_lw_region = od_lw_wall * (1.0_jprb - this%lw_emissivity(:,iwall,jcol)) &
                       &        / max(od_lw_region,1.0e-6_jprb)
                  ! We assume that any scattering off the walls is
                  ! equally likely to be up or down, so the effective
                  ! asymmetry factor is zero
                  g_lw_region = 0.0_jprb
                  ! Effective Planck function of the canopy is the
                  ! weighted average of wall and air emission
                  this%lw_total_wall_emission(:,iregion,jcol) = LwDiffusivity &
                       &  * od_lw_wall*this%lw_emissivity(:,iwall,jcol)*this%planck_facet(:,iwall,jcol)
                  this%lw_total_canopy_emission(:,iregion,jcol) = LwDiffusivity &
                       &  * od_lw_air*this%planck_region(:,iregion,jcol)
                  planck_canopy = (this%lw_total_wall_emission(:,iregion,jcol) &
                       &           +this%lw_total_canopy_emission(:,iregion,jcol)) &
                       &        / max(od_lw_region*(1.0_jprb-ssa_lw_region)*LwDiffusivity,1.0e-6_jprb)
                end if

                ! Compute fraction of canyon absorption by the wall 
                this%lw_wall_abs_frac(:,iregion,jcol) = od_lw_wall * this%lw_emissivity(:,iwall,jcol) &
                     &    / max(od_lw_region*(1.0_jprb-ssa_lw_region),1.0e-6_jprb)


              else
                od_lw_region  = od_lw_wall
                ssa_lw_region = 1.0_jprb - this%lw_emissivity(:,iwall,jcol);
                g_lw_region   = 0.0_jprb

                ! All absorption and emission is by the wall
                this%lw_wall_abs_frac(:,iregion,jcol) = 1.0_jprb
                this%lw_total_wall_emission(:,iregion,jcol) = LwDiffusivity &
                     &  * od_lw_wall*this%lw_emissivity(:,iwall,jcol)*this%planck_facet(:,iwall,jcol)
                this%lw_total_canopy_emission(:,iregion,jcol) = 0.0_jprb
                planck_canopy = this%planck_facet(:,iwall,jcol)
              end if

              call calc_two_stream_gammas_lw(this%nlwbands, &
                   &  ssa_lw_region, g_lw_region, &
                   &  gamma1_lw, gamma2_lw)
              call calc_reflectance_transmittance_isothermal_lw(this%nlwbands, &
                   &  od_lw_region, gamma1_lw, gamma2_lw, planck_canopy, &
                   &  this%lw_reflectance(:,iregion,jcol), this%lw_transmittance(:,iregion,jcol), &
                   &  this%lw_source(:,iregion,jcol))

              ! Add roof component
              single_level%lw_emissivity(jcol,:) = single_level%lw_emissivity(jcol,:) &
                   &  + tile_fraction * surface%building_fraction(jcol,jtile) &
                   &  * this%lw_emissivity(:,iroof,jcol)
              single_level%lw_emission(jcol,:) = single_level%lw_emission(jcol,:) &
                   &  + tile_fraction * surface%building_fraction(jcol,jtile) &
                   &  * this%lw_emissivity(:,iroof,jcol)*this%planck_facet(:,iroof,jcol)
              
              ! Add canyon component: longwave adding method for a single layer
              inv_denominator_lw(:) = 1.0_jprb &
                   &  / (1.0_jprb - (1.0_jprb - this%lw_emissivity(:,iground,jcol)) &
                   &               *this%lw_reflectance(:,iregion,jcol))
              this%lw_emissivity_region(:,iregion,jcol) &
                   &  = (1.0_jprb - (this%lw_reflectance(:,iregion,jcol) &
                   &  + this%lw_transmittance(:,iregion,jcol)**2 &
                   &  * (1.0_jprb - this%lw_emissivity(:,iground,jcol)) * inv_denominator_lw(:)))
              single_level%lw_emissivity(jcol,:) = single_level%lw_emissivity(jcol,:) &
                   &  + tile_fraction * canyon_fraction * this%lw_emissivity_region(:,iregion,jcol) 
              this%lw_emission_region(:,iregion,jcol) &
                   &  = this%lw_source(:,iregion,jcol) * (1.0_jprb + inv_denominator_lw(:) &
                   &       * (1.0_jprb - this%lw_emissivity(:,iground,jcol)) &
                   &       * this%lw_transmittance(:,iregion,jcol)) &
                   &     +this%planck_facet(:,iground,jcol)*this%lw_emissivity(:,iground,jcol) &
                   &     *this%lw_transmittance(:,iregion,jcol)*inv_denominator_lw(:)
              single_level%lw_emission(jcol,:) = single_level%lw_emission(jcol,:) &
                   &  + tile_fraction * canyon_fraction * this%lw_emission_region(:,iregion,jcol)

            end select
          end if
        end do
      end do
    end if

    if (lhook) call dr_hook('radsurf_intermediate:calc_boundary_conditions_lw',1,hook_handle)

  end subroutine calc_boundary_conditions_lw


  !---------------------------------------------------------------------
  ! Compute both shortwave and longwave boundary conditions neglecting
  ! gases within vegetation and urban canopies
  subroutine calc_boundary_conditions_vacuum(this, istartcol, iendcol, &   ! in
       &                              config, surface,          &   ! in
       &                              single_level)                 ! out

    use radiation_config,       only : config_type
    use radiation_single_level, only : single_level_type
    use radsurf_properties,     only : surface_type

    class(surface_intermediate_type), intent(inout) :: this
    integer(kind=jpim),               intent(in)    :: istartcol, iendcol
    type(config_type),                intent(in)    :: config
    type(surface_type),               intent(in)    :: surface
    type(single_level_type),          intent(inout) :: single_level

    call this%calc_boundary_conditions_sw(istartcol, iendcol, &
         &                                config, surface, single_level)
    call this%calc_boundary_conditions_lw(istartcol, iendcol, &
         &                                config, surface, single_level)

  end subroutine calc_boundary_conditions_vacuum


  !---------------------------------------------------------------------
  ! Compute both shortwave and longwave boundary conditions
  subroutine calc_boundary_conditions(this, istartcol, iendcol, &   ! in
       &                              config, surface,          &   ! in
       &                              thermodynamics, gas,      &   ! in
       &                              single_level)                 ! out

    use radiation_config,       only : config_type
    use radiation_thermodynamics,only: thermodynamics_type
    use radiation_gas,          only : gas_type
    use radiation_single_level, only : single_level_type
    use radiation_ifs_rrtm,     only : gas_optics
    use radsurf_properties,     only : surface_type
    use radiation_constants,    only : GasConstantDryAir, &
         &                             AccelDueToGravity

    class(surface_intermediate_type), intent(inout) :: this
    integer(kind=jpim),               intent(in)    :: istartcol, iendcol
    type(config_type),                intent(in)    :: config
    type(surface_type),               intent(in)    :: surface
    type(gas_type),                   intent(in)    :: gas
    type(thermodynamics_type),        intent(in)    :: thermodynamics
    type(single_level_type),          intent(inout) :: single_level

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

    ! Number of leves to request gas optical depths for; we only need
    ! one but gas optics scheme assumes more
    integer, parameter :: NGasLevels = 1

    ! Canopy longwave optical depth
    real(jprb), dimension(config%n_g_lw,NGasLevels,istartcol:iendcol) :: od_lw

    ! Canopy longwave extinction coefficient
    real(jprb), dimension(config%n_g_lw,istartcol:iendcol) :: ext_lw

    ! Canopy shortwave optical depth, single scattering albedo and
    ! asymmetry factor of gases and aerosols at each shortwave g-point
    real(jprb), dimension(config%n_g_sw,NGasLevels,istartcol:iendcol) :: od_sw, ssa_sw, g_sw

    ! Thickness of lowest model level
    real(jprb) :: layer_depth

    ! Index to final level of full model grid
    integer :: iendlev

    ! Column loop counter, number of columns
    integer :: jcol, ncol

    ncol    = ubound(thermodynamics%pressure_hl,1)
    iendlev = ubound(thermodynamics%pressure_hl,2)-1

    call gas_optics(ncol, NGasLevels, istartcol, iendcol, &
         &  config, single_level, thermodynamics, gas, & 
         &  od_lw, od_sw, ssa_sw)

    ! Scale optical depths to extinction
    do jcol = istartcol,iendcol
      layer_depth = R_over_g &
             &  * (thermodynamics%pressure_hl(jcol,iendlev+1) &
             &     - thermodynamics%pressure_hl(jcol,iendlev)) &
             &  * (thermodynamics%temperature_hl(jcol,iendlev) &
             &     + thermodynamics%temperature_hl(jcol,iendlev+1)) &
             &  / (thermodynamics%pressure_hl(jcol,iendlev) &
             &     + thermodynamics%pressure_hl(jcol,iendlev+1))
      !if (config%do_sw) then
      !  ext_sw(jcol) = ext_sw(jcol) / layer_depth
      !end if
      if (config%do_lw) then
        ext_lw(:,jcol) = od_lw(:,NGasLevels,jcol) / layer_depth
      end if
    end do

    call this%calc_boundary_conditions_sw(istartcol, iendcol, &
         &                                config, surface, single_level)
    if (config%do_canopy_gases_lw) then
      call this%calc_boundary_conditions_lw(istartcol, iendcol, &
           &                                config, surface, single_level, &
           &                                ext_lw_air=ext_lw)
    else
      call this%calc_boundary_conditions_lw(istartcol, iendcol, &
           &                                config, surface, single_level)
    end if

  end subroutine calc_boundary_conditions


  !---------------------------------------------------------------------
  subroutine partition_fluxes(this, istartcol, iendcol, config, surface, flux, &
       &                      surface_flux)

    use yomhook,           only : lhook, dr_hook

    use radiation_flux,    only : flux_type
    use radsurf_flux,      only : surface_flux_type
    use radiation_config,  only : config_type
    use radsurf_properties,only : surface_type, ITileFlat,ITileVegetation,ITileUrban3D

    class(surface_intermediate_type), intent(in)    :: this
    integer(kind=jpim),               intent(in)    :: istartcol, iendcol
    type(config_type),                intent(in)    :: config
    type(surface_type),               intent(in)    :: surface
    type(flux_type),                  intent(in)    :: flux
    type(surface_flux_type),          intent(inout) :: surface_flux

    real(kind=jprb), dimension(config%n_g_lw)       :: lw_dn_g, lw_up_g
    real(kind=jprb), dimension(config%n_canopy_bands_sw) &
         &  :: sw_dn_diffuse_g, sw_dn_direct_g, sw_up_g, sw_abs_g
    real(kind=jprb), dimension(config%n_canopy_bands_lw) :: lw_abs_g

    ! Ratio of street planar area to wall frontal area
    real(kind=jprb)    :: wall_scaling

    real(kind=jprb)    :: tile_fraction
    integer(kind=jpim) :: iground, iroof, iwall, iregion
    integer(kind=jpim) :: isurf ! index to lowest flux level (=nlev+1)
    integer(kind=jpim) :: jcol, jtile

    real(kind=jprb) :: hook_handle

    if (lhook) call dr_hook('radsurf_intermediate:partition_fluxes',0,hook_handle)

    if (.not. surface%is_simple) then
      isurf = size(flux%lw_dn,2)

      if (config%do_lw) then
        surface_flux%lw_dn_facet(istartcol:iendcol,:) = 0.0_jprb
        surface_flux%lw_up_facet(istartcol:iendcol,:) = 0.0_jprb
        surface_flux%lw_abs_canopy(istartcol:iendcol,:)=0.0_jprb
      end if
      if (config%do_sw) then
        surface_flux%sw_dn_facet(istartcol:iendcol,:) = 0.0_jprb
        surface_flux%sw_dn_direct_facet(istartcol:iendcol,:) = 0.0_jprb
        surface_flux%sw_up_facet(istartcol:iendcol,:) = 0.0_jprb
        surface_flux%sw_abs_canopy(istartcol:iendcol,:)=0.0_jprb
      end if

      ! Loop over column and tile
      do jcol = istartcol,iendcol
        do jtile = 1,surface%ntile
          tile_fraction = surface%tile_fraction(jcol,jtile)

          if (tile_fraction > 0.0_jprb) then

            select case (surface%i_representation(jtile))
            case (ITileFlat)
              iground = surface%i_ground_facet(jcol,jtile)
              if (config%do_lw) then
                surface_flux%lw_dn_facet(jcol,iground) = flux%lw_dn(jcol,isurf)
                surface_flux%lw_up_facet(jcol,iground) &
                     &  = sum(this%lw_emissivity(:,iground,jcol)*this%planck_facet(:,iground,jcol) &
                     &        + (1.0_jprb - this%lw_emissivity(:,iground,jcol)) &
                     &        * flux%lw_dn_surf_canopy(:,jcol))
              end if
              if (config%do_sw) then
                surface_flux%sw_dn_facet(jcol,iground)        = flux%sw_dn(jcol,isurf)
                surface_flux%sw_dn_direct_facet(jcol,iground) = flux%sw_dn_direct(jcol,isurf)
                surface_flux%sw_up_facet(jcol,iground) &
                     & = sum(this%sw_albedo_diffuse(:,iground,jcol) &
                     &       * flux%sw_dn_diffuse_surf_canopy(:,jcol) &
                     &       + this%sw_albedo_direct (:,iground,jcol) &
                     &       * flux%sw_dn_direct_surf_canopy (:,jcol))
              end if

            case (ITileVegetation)
              iground  = surface%i_ground_facet(jcol,jtile)
              iregion = surface%i_region_1(jcol,jtile)
              ! Adding method in longwave and shortwave
              if (config%do_lw) then
                ! Surface downwelling fluxes at each g point
                lw_dn_g = (this%lw_transmittance(:,iregion,jcol)*flux%lw_dn_surf_canopy(:,jcol) &
                     &    +this%lw_reflectance  (:,iregion,jcol) &
                     &     *this%lw_emissivity(:,iground,jcol)*this%planck_facet(:,iground,jcol) &
                     &    +this%lw_source(:,iregion,jcol)) &
                     &  / (1.0_jprb - (1.0_jprb - this%lw_emissivity(:,iground,jcol)) &
                     &               *this%lw_reflectance(:,iregion,jcol))
                surface_flux%lw_dn_facet(jcol,iground) = sum(lw_dn_g)
                surface_flux%lw_up_facet(jcol,iground) &
                     &  = sum((1.0_jprb-this%lw_emissivity(:,iground,jcol))*lw_dn_g &
                     &       +this%lw_emissivity(:,iground,jcol)*this%planck_facet(:,iground,jcol))
                surface_flux%lw_abs_canopy(jcol,jtile) = flux%lw_dn(jcol,isurf) &
                     &  - flux%lw_up(jcol,isurf) - surface_flux%lw_dn_facet(jcol,iground) &
                     &  + surface_flux%lw_up_facet(jcol,iground)

              end if
              if (config%do_sw) then
                sw_dn_direct_g = this%sw_tra_dir_dir(:,iregion,jcol) &
                     &  * flux%sw_dn_direct_surf_canopy(:,jcol)
                ! Note that the following is initially just the
                ! upwelling due to scattering of the direct beam
                sw_up_g = sw_dn_direct_g * this%sw_albedo_direct(:,iground,jcol)
                sw_dn_diffuse_g &
                     &  = (this%sw_tra_dif(:,iregion,jcol)*flux%sw_dn_diffuse_surf_canopy(:,jcol) &
                     &    +this%sw_ref_dif(:,iregion,jcol)*sw_up_g &
                     &    +this%sw_tra_dir_dif(:,iregion,jcol)*flux%sw_dn_direct_surf_canopy(:,jcol)) &
                     &  / (1.0_jprb - this%sw_albedo_diffuse(:,iground,jcol) & 
                     &               *this%sw_ref_dif(:,iregion,jcol))
                sw_up_g = sw_up_g + sw_dn_diffuse_g * this%sw_albedo_diffuse(:,iground,jcol)
                surface_flux%sw_dn_direct_facet(jcol,iground) = sum(sw_dn_direct_g)
                surface_flux%sw_dn_facet(jcol,iground) = surface_flux%sw_dn_direct_facet(jcol,iground) &
                     &  + sum(sw_dn_diffuse_g)
                surface_flux%sw_up_facet(jcol,iground) = sum(sw_up_g)
                surface_flux%sw_abs_canopy(jcol,jtile) &
                     &  = flux%sw_dn(jcol,isurf) - flux%sw_up(jcol,isurf) &
                     &  - surface_flux%sw_dn_facet(jcol,iground) &
                     &  + surface_flux%sw_up_facet(jcol,iground)
              end if
            case (ITileUrban3D)
              iground = surface%i_ground_facet(jcol,jtile)
              iroof   = surface%i_roof_facet(jcol,jtile)
              iwall   = surface%i_wall_facet(jcol,jtile)
              iregion = surface%i_region_1(jcol,jtile)

              ! We want wall fluxes per unit area of wall, not per
              ! unit area of street
              wall_scaling = (1.0_jprb - surface%building_fraction(jcol,jtile)) &
                   &  / max(1.0e-4_jprb, surface%building_normalized_perimeter(jcol,jtile) &
                   &                     * surface%canopy_depth(jcol,jtile))

              if (config%do_sw) then
                ! Roof fluxes
                surface_flux%sw_dn_facet(jcol,iroof)        = flux%sw_dn(jcol,isurf)
                surface_flux%sw_dn_direct_facet(jcol,iroof) = flux%sw_dn_direct(jcol,isurf)
                surface_flux%sw_up_facet(jcol,iroof) &
                     & = sum(this%sw_albedo_diffuse(:,iroof,jcol) &
                     &       * flux%sw_dn_diffuse_surf_canopy(:,jcol) &
                     &       + this%sw_albedo_direct (:,iroof,jcol) &
                     &       * flux%sw_dn_direct_surf_canopy (:,jcol))

                ! Ground fluxes
                sw_dn_direct_g = this%sw_tra_dir_dir(:,iregion,jcol) &
                     &         * flux%sw_dn_direct_surf_canopy(:,jcol)
                ! Note that the following is initially just the
                ! upwelling due to scattering of the direct beam
                sw_up_g = sw_dn_direct_g * this%sw_albedo_direct(:,iground,jcol)
                sw_dn_diffuse_g &
                     &  = (this%sw_tra_dif(:,iregion,jcol)*flux%sw_dn_diffuse_surf_canopy(:,jcol) &
                     &    +this%sw_ref_dif(:,iregion,jcol)*sw_up_g &
                     &    +this%sw_tra_dir_dif(:,iregion,jcol)*flux%sw_dn_direct_surf_canopy(:,jcol)) &
                     &  / (1.0_jprb - this%sw_albedo_diffuse(:,iground,jcol) & 
                     &               *this%sw_ref_dif(:,iregion,jcol))
                sw_up_g = sw_up_g + sw_dn_diffuse_g * this%sw_albedo_diffuse(:,iground,jcol)
                surface_flux%sw_dn_direct_facet(jcol,iground) = sum(sw_dn_direct_g)
                surface_flux%sw_dn_facet(jcol,iground) = surface_flux%sw_dn_direct_facet(jcol,iground) &
                     &  + sum(sw_dn_diffuse_g)
                surface_flux%sw_up_facet(jcol,iground) = sum(sw_up_g)

                ! Wall fluxes
                sw_abs_g = flux%sw_dn_direct_surf_canopy(:,jcol)*this%sw_wall_abs_dir(:,iregion,jcol)
                surface_flux%sw_dn_direct_facet(jcol,iwall) &
                     &  = wall_scaling * sum(sw_abs_g / (1.0_jprb - this%sw_albedo_direct(:,iwall,jcol)))
                ! Initially just the direct reflection
                surface_flux%sw_up_facet(jcol,iwall) = wall_scaling &
                     &  * sum(sw_abs_g * this%sw_albedo_direct(:,iwall,jcol) &
                     &        / (1.0_jprb - this%sw_albedo_direct(:,iwall,jcol)))
                ! Initially just the direct absorption
                surface_flux%sw_abs_canopy(jcol,jtile) &
                     &  = sum(flux%sw_dn_direct_surf_canopy(:,jcol)*this%sw_air_abs_dir(:,iregion,jcol))

                ! Diffuse absorption within canopy
                sw_abs_g = flux%sw_dn_direct_surf_canopy(:,jcol) &
                     &   * (1.0-this%sw_albedo_direct_reg (:,iregion,jcol)) &
                     &   + flux%sw_dn_diffuse_surf_canopy(:,jcol) &
                     &   * (1.0-this%sw_albedo_diffuse_reg(:,iregion,jcol)) &
                     &   - sw_dn_direct_g - sw_dn_diffuse_g + sw_up_g - sw_abs_g
                ! Add diffuse absorption
                surface_flux%sw_abs_canopy(jcol,jtile) = surface_flux%sw_abs_canopy(jcol,jtile) &
                     &   + sum(sw_abs_g * (1.0_jprb - this%sw_wall_abs_frac_dif(:,iregion,jcol)))
                ! Add diffuse reflection
                surface_flux%sw_up_facet(jcol,iwall) = surface_flux%sw_up_facet(jcol,iwall) &
                     &   +  wall_scaling * sum(sw_abs_g * this%sw_wall_abs_frac_dif(:,iregion,jcol) &
                     &         * this%sw_albedo_diffuse(:,iwall,jcol) &
                     &         / (1.0_jprb - this%sw_albedo_diffuse(:,iwall,jcol)))
                ! Add diffuse into wall
                surface_flux%sw_dn_facet(jcol,iwall) = surface_flux%sw_dn_direct_facet(jcol,iwall) &
                     &   + wall_scaling * sum(sw_abs_g * this%sw_wall_abs_frac_dif(:,iregion,jcol) &
                     &         / (1.0_jprb - this%sw_albedo_diffuse(:,iwall,jcol)))
!!$
!!$                sw_direct_abs_g  = flux%sw_dn_direct_surf_g (:,jcol) - sw_dn_direct_g
!!$                sw_diffuse_abs_g = flux%sw_dn_diffuse_surf_g(:,jcol)*(1.0_jprb - this%sw_albedo_diffuse_reg(:,iregion,jcol)) &
!!$                     &   + flux%sw_dn_direct_surf_g (:,jcol)*(1.0_jprb - this%sw_albedo_direct_reg (:,iregion,jcol)) &
!!$                     &   - sw_dn_direct_g - sw_dn_diffuse_g + sw_up_g
!!$
!!$
!!$                surface_flux%sw_abs_canopy(jcol,jtile) = sum((1.0_jprb - this%sw_wall_abs_frac(:,iregion,jcol))*sw_abs_g)
!!$                surface_flux%sw_dn_facet(jcol,iwall) = sum(this%sw_wall_abs_frac(:,iregion,jcol)*sw_abs_g &
!!$                     &                           / (1.0_jprb - this%sw_albedo_diffuse(:,iwall,jcol)))
!!$                surface_flux%sw_up_facet(jcol,iwall) = sum(this%sw_albedo_diffuse(:,iwall,jcol) &
!!$                     &                             *this%sw_wall_abs_frac(:,iregion,jcol)*sw_abs_g &
!!$                     &                           / (1.0_jprb - this%sw_albedo_diffuse(:,iwall,jcol)))
              end if

              if (config%do_lw) then
                surface_flux%lw_dn_facet(jcol,iroof) = flux%lw_dn(jcol,isurf)
                surface_flux%lw_up_facet(jcol,iroof) &
                     &  = sum(this%lw_emissivity(:,iroof,jcol)*this%planck_facet(:,iroof,jcol) &
                     &  +  (1.0_jprb - this%lw_emissivity(:,iroof,jcol))*flux%lw_dn_surf_canopy(:,jcol))
                lw_dn_g = (this%lw_transmittance(:,iregion,jcol)*flux%lw_dn_surf_canopy(:,jcol) &
                     &    +this%lw_reflectance  (:,iregion,jcol) &
                     &     *this%lw_emissivity(:,iground,jcol)*this%planck_facet(:,iground,jcol) &
                     &    +this%lw_source(:,iregion,jcol)) &
                     &  / (1.0_jprb - (1.0_jprb - this%lw_emissivity(:,iground,jcol)) &
                     &               *this%lw_reflectance(:,iregion,jcol))
                lw_up_g = (1.0_jprb-this%lw_emissivity(:,iground,jcol))*lw_dn_g &
                     &       +this%lw_emissivity(:,iground,jcol)*this%planck_facet(:,iground,jcol)
                surface_flux%lw_dn_facet(jcol,iground) = sum(lw_dn_g)
                surface_flux%lw_up_facet(jcol,iground) = sum(lw_up_g)
!                lw_abs_g = flux%lw_dn(jcol,isurf)-flux%lw_up(jcol,isurf) &
!                     &  - surface_flux%lw_dn_facet(jcol,iground)+surface_flux%lw_up_facet(jcol,iground)
!                lw_abs_g = flux%lw_dn_surf_canopy(:,jcol) &
!                     &  * (1.0_jprb-this%lw_emissivity_region(:,iregion,jcol)) &
!                     &  - this%lw_emission_region(:,iregion,jcol) &
!                     &  - lw_dn_g*(1.0_jprb-this%lw_emissivity(:,iground,jcol)) &
!                     &  + this%lw_emissivity(:,iground,jcol)*this%planck_facet(:,iground,jcol)
!                lw_abs_g = flux%lw_dn_surf_canopy(:,jcol) &
!                     &  * this%lw_emissivity_region(:,iregion,jcol) &
!                     &  - this%lw_emission_region(:,iregion,jcol) &
!                     &  - lw_dn_g*this%lw_emissivity(:,iground,jcol) &
!                     &  + this%lw_emissivity(:,iground,jcol)*this%planck_facet(:,iground,jcol)
                lw_abs_g = (flux%lw_dn_surf_canopy(:,jcol) + lw_up_g) &
                     &   * (1.0_jprb - this%lw_reflectance(:,iregion,jcol) &
                     &               - this%lw_transmittance(:,iregion,jcol)) &
                     &   + this%lw_total_wall_emission(:,iregion,jcol) &
                     &   + this%lw_total_canopy_emission(:,iregion,jcol) &
                     &   - 2.0_jprb * this%lw_source(:,iregion,jcol)

!                surface_flux%lw_up_facet(jcol,iwall) = sum(this%planck_facet(:,iwall,jcol) &
!                     &  + wall_scaling*(1.0_jprb-this%lw_emissivity(:,iwall,jcol)) &
!                     &    *this%lw_wall_abs_frac(:,iregion,jcol)*lw_abs_g &
!                     &    / this%lw_emissivity(:,iwall,jcol))
!                surface_flux%lw_dn_facet(jcol,iwall) = sum((this%planck_facet(:,iwall,jcol) &
!                     &  + this%lw_total_wall_emission(:,iregion,jcol)) / this%lw_emissivity(:,iwall,jcol))

!                surface_flux%lw_dn_facet(jcol,iwall) = wall_scaling*sum((this%lw_total_wall_emission(:,iregion,jcol) &
!                     &  + this%lw_wall_abs_frac(:,iregion,jcol)*lw_abs_g) &
!                     &  / this%lw_emissivity(:,iwall,jcol))
!                surface_flux%lw_up_facet(jcol,iwall) = surface_flux%lw_dn_facet(jcol,iwall) &
!                     &  - wall_scaling*sum(this%lw_wall_abs_frac(:,iregion,jcol)*lw_abs_g)

                surface_flux%lw_dn_facet(jcol,iwall) &
                     &  = wall_scaling*sum(this%lw_wall_abs_frac(:,iregion,jcol)*lw_abs_g &
                     &                     / this%lw_emissivity(:,iwall,jcol))
                surface_flux%lw_up_facet(jcol,iwall) = surface_flux%lw_dn_facet(jcol,iwall) &
                     &  + wall_scaling * sum(this%lw_total_wall_emission(:,iregion,jcol) &
                     &                       - this%lw_wall_abs_frac(:,iregion,jcol)*lw_abs_g)

                surface_flux%lw_abs_canopy(jcol,jtile) &
                     &  = sum(lw_abs_g*(1.0_jprb-this%lw_wall_abs_frac(:,iregion,jcol)) &
                     &        -this%lw_total_canopy_emission(:,iregion,jcol))
              end if
            end select

          end if

        end do
      end do
    end if

    if (lhook) call dr_hook('radsurf_intermediate:partition_fluxes',1,hook_handle)

  end subroutine partition_fluxes

end module radsurf_intermediate