source: LMDZ6/trunk/libf/phylmd/ecrad/radiation_tripleclouds_lw.F90 @ 4388

! radiation_tripleclouds_lw.F90 - Longwave "Tripleclouds" solver
!
! (C) Copyright 2016- 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-28  R. Hogan  Receive emission/albedo rather than planck/emissivity
!   2017-04-22  R. Hogan  Store surface fluxes at all g-points
!   2017-10-23  R. Hogan  Renamed single-character variables
!   2018-10-08  R. Hogan  Call calc_region_properties

module radiation_tripleclouds_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
  ! "Tripleclouds" solver in which cloud inhomogeneity is treated by
  ! dividing each model level into three regions, one clear and two
  ! cloudy (with differing optical depth). This approach was described
  ! by Shonk and Hogan (2008).

  subroutine solver_tripleclouds_lw(nlev,istartcol,iendcol, &
       &  config, 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

!    use radiation_io, only             : nulout
    use radiation_config, only         : config_type, IPdfShapeGamma
    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, add_indexed_sum
    use radiation_matrix, only         : singlemat_x_vec
    use radiation_two_stream, only     : calc_two_stream_gammas_lw, &
         &                               calc_reflectance_transmittance_lw, &
         &                               calc_no_scattering_transmittance_lw
    use radiation_lw_derivatives, only : calc_lw_derivatives_region

    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(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
    ! and at the surface at each longwave g-point
    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

    ! Optical depth, single scattering albedo and asymmetry factor in
    ! each g-point including gas, aerosol and clouds
    real(jprb), dimension(config%n_g_lw) :: od_total, ssa_total, g_total

    ! Modified optical depth after Tripleclouds scaling to represent
    ! cloud inhomogeneity
    real(jprb), dimension(config%n_g_lw) :: od_cloud_new

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

    ! Local constants
    integer, parameter :: nregions = 3

    ! In a clear-sky layer this will be 1, otherwise equal to nregions
    integer :: nreg

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

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

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

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

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

    ! Emission by a layer into the upwelling or downwelling diffuse
    ! streams
    real(jprb), dimension(config%n_g_lw, nregions, nlev) &
         &  :: Sup, Sdn

    ! ...clear-sky equivalent
    real(jprb), dimension(config%n_g_lw, nlev) &
         &  :: Sup_clear, Sdn_clear

    ! 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, nregions, nlev+1) :: total_albedo

    ! 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, nregions, nlev+1) :: total_source

    ! ...equivalent values for clear-skies
    real(jprb), dimension(config%n_g_lw, nlev+1) :: total_albedo_clear, total_source_clear

    ! Total albedo and source of the atmosphere just below a layer interface
    real(jprb), dimension(config%n_g_lw, nregions) &
         &  :: total_albedo_below, total_source_below

    ! Downwelling flux below and above an interface between
    ! layers into a plane perpendicular to the direction of the sun
    real(jprb), dimension(config%n_g_lw, nregions) &
         &  :: flux_dn, flux_dn_below, flux_up

    ! ...clear-sky equivalent (no distinction between "above/below")
    real(jprb), dimension(config%n_g_lw) &
         &  :: flux_dn_clear, flux_up_clear

    ! Clear-sky equivalent, but actually its reciprocal to replace
    ! some divisions by multiplications
    real(jprb), dimension(config%n_g_lw, nregions) :: inv_denom

    ! 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)

    integer :: jcol, jlev, jg, jreg, jreg2, ng

    real(jprb) :: hook_handle

    if (lhook) call dr_hook('radiation_tripleclouds_lw:solver_tripleclouds_lw',0,hook_handle)

    ! --------------------------------------------------------
    ! Section 1: Prepare general variables and arrays
    ! --------------------------------------------------------
    ! Copy array dimensions to local variables for convenience
    ng   = config%n_g_lw

    ! Compute the wavelength-independent region fractions and
    ! optical-depth scalings
    call calc_region_properties(nlev,nregions,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,nregions,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)

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

      ! 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: Loop over layers to compute reflectance and transmittance
      ! --------------------------------------------------------
      ! In this section the reflectance, transmittance and sources
      ! are computed for each layer
      do jlev = 1,nlev ! Start at top-of-atmosphere

        ! Array-wise assignments
        gamma1 = 0.0_jprb
        gamma2 = 0.0_jprb

        nreg = nregions
        if (is_clear_sky_layer(jlev)) then
          nreg = 1
          reflectance(:,2:,jlev)   = 0.0_jprb
          transmittance(:,2:,jlev)   = 0.0_jprb
          Sup(:,2:,jlev) = 0.0_jprb
          Sdn(:,2:,jlev) = 0.0_jprb
        end if
        do jreg = 1,nreg
          if (jreg == 1) then
            ! Clear-sky calculation
            if (.not. config%do_lw_aerosol_scattering) then
              call calc_no_scattering_transmittance_lw(ng, od(:,jlev,jcol), &
                   &  planck_hl(:,jlev,jcol), planck_hl(:,jlev+1, jcol), &
                   &  transmittance(:,1,jlev), Sup(:,1,jlev), Sdn(:,1,jlev))
              reflectance(:,1,jlev) = 0.0_jprb
            else
              call calc_two_stream_gammas_lw(ng, &
                   &  ssa(:,jlev,jcol), g(:,jlev,jcol), gamma1, gamma2)
              call calc_reflectance_transmittance_lw(ng, &
                   &  od(:,jlev,jcol), gamma1, gamma2, &
                   &  planck_hl(:,jlev,jcol), planck_hl(:,jlev+1,jcol), &
                   &  reflectance(:,1,jlev), transmittance(:,1,jlev), &
                   &  Sup(:,1,jlev), Sdn(:,1,jlev))
            end if
          else
            ! Cloudy sky
            ! Add scaled cloud optical depth to clear-sky value
            od_cloud_new = od_cloud(config%i_band_from_reordered_g_lw,jlev,jcol) &
                 &       * od_scaling(jreg,jlev,jcol)
            od_total = od(:,jlev,jcol) + od_cloud_new

            if (config%do_lw_cloud_scattering) then
              ssa_total = 0.0_jprb
              g_total   = 0.0_jprb              
              if (config%do_lw_aerosol_scattering) then
                where (od_total > 0.0_jprb)
                  ssa_total = (ssa(:,jlev,jcol)*od(:,jlev,jcol) &
                       &     + ssa_cloud(config%i_band_from_reordered_g_lw,jlev,jcol) &
                       &     *  od_cloud_new) & 
                       &     / od_total
                end where
                where (ssa_total > 0.0_jprb .and. od_total > 0.0_jprb)
                  g_total = (g(:,jlev,jcol)*ssa(:,jlev,jcol)*od(:,jlev,jcol) &
                       &     +   g_cloud(config%i_band_from_reordered_g_lw,jlev,jcol) &
                       &     * ssa_cloud(config%i_band_from_reordered_g_lw,jlev,jcol) &
                       &     *  od_cloud_new) &
                       &     / (ssa_total*od_total)
                end where
              else
                where (od_total > 0.0_jprb)
                  ssa_total = ssa_cloud(config%i_band_from_reordered_g_lw,jlev,jcol) &
                       &     * od_cloud_new / od_total
                end where
                where (ssa_total > 0.0_jprb .and. od_total > 0.0_jprb)
                  g_total = g_cloud(config%i_band_from_reordered_g_lw,jlev,jcol) &
                       &     * ssa_cloud(config%i_band_from_reordered_g_lw,jlev,jcol) &
                       &     *  od_cloud_new / (ssa_total*od_total)
                end where
              end if
              call calc_two_stream_gammas_lw(ng, &
                   &  ssa_total, g_total, gamma1, gamma2)
              call calc_reflectance_transmittance_lw(ng, &
                   &  od_total, gamma1, gamma2, &
                   &  planck_hl(:,jlev,jcol), planck_hl(:,jlev+1,jcol), &
                   &  reflectance(:,jreg,jlev), transmittance(:,jreg,jlev), &
                   &  Sup(:,jreg,jlev), Sdn(:,jreg,jlev))
            else
              ! No-scattering case: use simpler functions for
              ! transmission and emission
              call calc_no_scattering_transmittance_lw(ng, od_total, &
                   &  planck_hl(:,jlev,jcol), planck_hl(:,jlev+1, jcol), &
                   &  transmittance(:,jreg,jlev), Sup(:,jreg,jlev), Sdn(:,jreg,jlev))
              reflectance(:,jreg,jlev) = 0.0_jprb
            end if
          end if
        end do

        ! Copy over the clear-sky emission
        Sup_clear(:,jlev) = Sup(:,1,jlev)
        Sdn_clear(:,jlev) = Sdn(:,1,jlev)
        if (.not. is_clear_sky_layer(jlev)) then
          ! Emission is scaled by the size of each region
          do jreg = 1,nregions
            Sup(:,jreg,jlev) = region_fracs(jreg,jlev,jcol) * Sup(:,jreg,jlev)
            Sdn(:,jreg,jlev) = region_fracs(jreg,jlev,jcol) * Sdn(:,jreg,jlev)
          end do
        end if

      end do ! Loop over levels

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

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

      ! Calculate the upwelling radiation emitted by the surface, and
      ! copy the surface albedo into total_albedo 
      do jreg = 1,nregions
        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,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,nlev+1)
      end if

      ! Work up from the surface computing the total albedo of the
      ! atmosphere and the total upwelling due to emission below each
      ! level below using the adding method
      do jlev = nlev,1,-1

        total_albedo_below        = 0.0_jprb

        if (config%do_clear) then
          ! For clear-skies there is no need to consider "above" and
          ! "below" quantities since with no cloud overlap to worry
          ! about, these are the same
          inv_denom(:,1) = 1.0_jprb &
               &  / (1.0_jprb - total_albedo_clear(:,jlev+1)*reflectance(:,1,jlev))
          total_albedo_clear(:,jlev) = reflectance(:,1,jlev) &
               &  + transmittance(:,1,jlev)*transmittance(:,1,jlev)*total_albedo_clear(:,jlev+1) &
               &  * inv_denom(:,1)
          total_source_clear(:,jlev) = Sup_clear(:,jlev) &
               &  + transmittance(:,1,jlev)*(total_source_clear(:,jlev+1) &
               &  + total_albedo_clear(:,jlev+1)*Sdn_clear(:,jlev)) &
               &  * inv_denom(:,1)
        end if

        if (is_clear_sky_layer(jlev)) then
          inv_denom(:,1) = 1.0_jprb &
               &  / (1.0_jprb - total_albedo(:,1,jlev+1)*reflectance(:,1,jlev))
          total_albedo_below = 0.0_jprb
          total_albedo_below(:,1) = reflectance(:,1,jlev) &
               &  + transmittance(:,1,jlev)*transmittance(:,1,jlev)*total_albedo(:,1,jlev+1) &
               &  * inv_denom(:,1)
          total_source_below = 0.0_jprb
          total_source_below(:,1) = Sup(:,1,jlev) &
               &  + transmittance(:,1,jlev)*(total_source(:,1,jlev+1) &
               &  + total_albedo(:,1,jlev+1)*Sdn(:,1,jlev)) &
               &  * inv_denom(:,1)
        else
          inv_denom = 1.0_jprb / (1.0_jprb - total_albedo(:,:,jlev+1)*reflectance(:,:,jlev))
          total_albedo_below = reflectance(:,:,jlev) &
               &  + transmittance(:,:,jlev)*transmittance(:,:,jlev)*total_albedo(:,:,jlev+1) &
               &  * inv_denom
          total_source_below = Sup(:,:,jlev) &
               &  + transmittance(:,:,jlev)*(total_source(:,:,jlev+1) &
               &  + total_albedo(:,:,jlev+1)*Sdn(:,:,jlev)) &
               &  * inv_denom
        end if

        ! Account for cloud overlap when converting albedo below a
        ! layer interface to the equivalent values just above
        if (is_clear_sky_layer(jlev) .and. is_clear_sky_layer(jlev-1)) then
          total_albedo(:,:,jlev) = total_albedo_below(:,:)
          total_source(:,:,jlev) = total_source_below(:,:)
        else
          total_source(:,:,jlev) = singlemat_x_vec(ng,ng,nregions,&
               &  u_matrix(:,:,jlev,jcol), total_source_below)
          ! Use overlap matrix and 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,nregions
            do jreg2 = 1,nregions
              total_albedo(:,jreg,jlev) &
                   &  = total_albedo(:,jreg,jlev) &
                   &  + total_albedo_below(:,jreg2) &
                   &  * v_matrix(jreg2,jreg,jlev,jcol)

            end do
          end do

        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 diffuse downwelling
      flux_dn = 0.0_jprb

      if (config%do_clear) then
        flux_dn_clear = 0.0_jprb
      end if

      ! Store the TOA broadband fluxes
      flux%lw_up(jcol,1) = sum(total_source(:,:,1))
      flux%lw_dn(jcol,1) = 0.0_jprb
      if (config%do_clear) then
        flux%lw_up_clear(jcol,1) = sum(total_source_clear(:,1))
        flux%lw_dn_clear(jcol,1) = 0.0_jprb
      end if

      ! Save the spectral fluxes if required
      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
          flux_dn_clear = (transmittance(:,1,jlev)*flux_dn_clear &
               &  + reflectance(:,1,jlev)*total_source_clear(:,jlev+1) + Sdn_clear(:,jlev) ) &
               &  / (1.0_jprb - reflectance(:,1,jlev)*total_albedo_clear(:,jlev+1))
          flux_up_clear = total_source_clear(:,jlev+1) &
               &        + flux_dn_clear*total_albedo_clear(:,jlev+1)
        end if

        if (is_clear_sky_layer(jlev)) then
          flux_dn(:,1) = (transmittance(:,1,jlev)*flux_dn(:,1) &
               &       + reflectance(:,1,jlev)*total_source(:,1,jlev+1) + Sdn(:,1,jlev) ) &
               &  / (1.0_jprb - reflectance(:,1,jlev)*total_albedo(:,1,jlev+1))
          flux_dn(:,2:)  = 0.0_jprb
          flux_up(:,1) = total_source(:,1,jlev+1) + flux_dn(:,1)*total_albedo(:,1,jlev+1)
          flux_up(:,2:)  = 0.0_jprb
        else
          flux_dn = (transmittance(:,:,jlev)*flux_dn &
               &     + reflectance(:,:,jlev)*total_source(:,:,jlev+1) + Sdn(:,:,jlev) ) &
               &  / (1.0_jprb - reflectance(:,:,jlev)*total_albedo(:,:,jlev+1))
          flux_up = total_source(:,:,jlev+1) + flux_dn*total_albedo(:,:,jlev+1)
        end if

        if (.not. (is_clear_sky_layer(jlev) &
             &    .and. is_clear_sky_layer(jlev+1))) then
          ! Account for overlap rules in translating fluxes just above
          ! a layer interface to the values just below
          flux_dn_below = singlemat_x_vec(ng,ng,nregions, &
               &  v_matrix(:,:,jlev+1,jcol), flux_dn)
          flux_dn = flux_dn_below
        end if ! Otherwise the fluxes in each region are the same so
               ! nothing to do

        ! Store the broadband fluxes
        flux%lw_up(jcol,jlev+1) = sum(sum(flux_up,1))
        flux%lw_dn(jcol,jlev+1) = sum(sum(flux_dn,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

        ! Save the spectral fluxes if required
        if (config%do_save_spectral_flux) then
          call indexed_sum(sum(flux_up,2), &
               &           config%i_spec_from_reordered_g_lw, &
               &           flux%lw_up_band(:,jcol,jlev+1))
          call indexed_sum(sum(flux_dn,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,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 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_region(ng, nlev, nregions, jcol, transmittance, &
             &  u_matrix(:,:,:,jcol), sum(flux_up,2), flux%lw_derivatives)
      end if

    end do ! Loop over columns

    if (lhook) call dr_hook('radiation_tripleclouds_lw:solver_tripleclouds_lw',1,hook_handle)

  end subroutine solver_tripleclouds_lw

end module radiation_tripleclouds_lw