scale_atm_phy_rd_dgm_simple.F90

Go to the documentation of this file.

!> module FElib / Atmosphere / Physics radiation
!!
!! @par Description
!!      A module to provide simplified radiation schemes based on Vallis et al. (2018), 
!!      which are a gray radiation scheme and a two-band longwave and one-band shortwave radiation scheme
!!  
!! @author Yuta Kawai, Team SCALE
!!
!! @par Reference
!!  - Vallis, G.K., Colyer, G., Geen, R., Gerber, E., Jucker, M., Maher, P.,
!!    Paterson, A., Pietschnig, M., Penn, J., and Thomson, S.I. 2018:
!!    Isca, v1.0: a framework for the global modelling of the atmospheres of Earth and other planets at varying levels of complexity.
!!    Geosci. Model Dev., 11, 843–859.
!!  - Frierson, D.M.W., Held, I.M., and Zurita-Gotor, P. 2006:
!!    A gray-radiation aquaplanet moist GCM. Part I: Static stability and eddy scale.
!!    J. Atmos. Sci., 63, 2548–2566.
!!  - Byrne, M.P., and O’Gorman, P.A. 2013:
!!    Land–ocean warming contrast over a wide range of climates:
!!    convective quasi-equilibrium theory and idealized simulations.
!!    J. Climate, 26, 4000–4016.
!!  - Geen, R., Czaja, A., and Haigh, J.D. 2016:
!!    The effects of increasing humidity on heat transport by extratropical waves.
!!    Geophys. Res. Lett., 43, 8314–8321.
!!
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_atm_phy_rd_dgm_simple
  !-----------------------------------------------------------------------------
  !
  !++ Used modules
  !
  use scale_precision
  use scale_io
  use scale_prc
  use scale_prof
 
  use scale_atmos_phy_rd_common, only: &
    i_up, i_dn, i_lw, i_sw
  
  use scale_element_base, only: &
    elementbase1d, elementbase2d, elementbase3d
  use scale_element_hexahedral, only: hexahedralelement
  use scale_localmesh_2d, only: localmesh2d  
  use scale_localmesh_3d, only: localmesh3d
  use scale_mesh_base3d, only: meshbase3d
  use scale_localmeshfield_base, only: localmeshfield3d
  use scale_meshfield_base, only: meshfield3d
 
  !-----------------------------------------------------------------------------
  implicit none
  private  
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedure
  !
 
 
  integer, parameter :: N_LW_BND_MAX = 2
 
  !> Derived type to represent a simplified radiation scheme based on Vallis et al. (2018)
  type, public :: atmphyradsimple
    integer :: optdep_type                  !< Type of optical depth calculation
 
    !- Parameters for optical depth calculation with a gray-radiation scheme
    real(rp) :: optdep_vallis_eq5_params(4) !< Parameters associated with optical depth calculation (Eq.5 in Vallis et al. (2018))
                                            !! A, mu, B, C
 
    !- Parameters for optical depth calculation with a two-band longwave and one-band shortwave radiation scheme
    real(rp) :: optdep_vallis_eq6a_params(2) !< Parameters associated with optical depth calculation (Eq.6a in Vallis et al. (2018))
                                             !! a_SW, c_SW
    real(rp) :: optdep_vallis_eq7a_params(4) !< Parameters associated with optical depth calculation (Eq.7a in Vallis et al. (2018))
                                             !! a_LW, b_LW, c_LW, d_LW
    real(rp) :: optdep_vallis_eq7b_params(4) !< Parameters associated with optical depth calculation (Eq.7b in Vallis et al. (2018))
                                             !! a_win, b_win, c_win, d_win
    integer :: optdep_vallis_eq6a_kadd          !< Number of layers above the model top in the SW optical depth calculation.
    real(rp) :: optdep_vallis_eq6a_qv_above_mot !< Specific humidity above the model top used in the SW optical depth calculation.
 
    real(rp) :: difffactor            !< Diffusivity factor for the two-stream approximation
    real(rp) :: frac_lw(n_lw_bnd_max) !< The fraction of the longwave spectrum at each band
    real(rp) :: pco2                  !< CO2 concentration in ppmv
 
  contains
    procedure, public :: init  => atm_phy_rd_dgm_simple_init
    procedure, public :: final => atm_phy_rd_dgm_simple_final
    procedure, public :: calculate_rad_flux  => atm_phy_rd_dgm_simple_flux
    procedure, private :: calc_optical_thick
  end type atmphyradsimple
  
  !-----------------------------------------------------------------------------
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !  
  integer, parameter :: optdep_type_v2018_gray = 1         !< An idealized radiation scheme based on Eq.5.
                                                           !! This is gray in infrared so that a single optical thickness is defined for the entire longwave spectrum, 
                                                           !! which includes a parameterization of long-wave absorption by CO2.
 
  integer, parameter :: optdep_type_v2018_lwbnd2swbnd1 = 2 !< An idealized radiation scheme based on Eqs. 6 and 7.
                                                           !! This scheme has two infrared bands and one shortwave band, as described in Green et al. (2016),
                                                           !! which provides an intermediate complexity between gray radiation and more sophisticated radiation schemes.
                                                           !! All bands were originally parameterized by fitting to data from SBDART for a range of atmospheric profiles. 
                                                           !! Vallis et al. (2018) include the CO2 absorption in each band and changes the functional form of the non-window optical depth.
                                                           !! The upward flux of shortwave radiation is assumed to be transparent.
 
contains
 
  !> Initialize an object to represent a gray-radiation scheme
!OCL SERIAL
  subroutine atm_phy_rd_dgm_simple_init( this )
    implicit none
    class(atmphyradsimple), intent(inout) :: this
 
    character(len=H_MID) :: simple_rd_type                  = 'Gray'  !< Type of radiation scheme. 'Gray', 'LWbnd2+SWbnd1'
    character(len=H_MID) :: v2018_gray_optdepth_params_type = 'V2018' !< How to specify the parameters for optical depth calculation (Eq.5 in Vallis et al. (2018)). 'BO2013', 'V2018', 'USER'
    real(rp) :: v2018_gray_optdepth_lw_params(4)                      !< Array of parameters for optical depth calculation with LW radiation (Eq.5 in Vallis et al. (2018)). A, mu, B, C
 
    character(len=H_MID) :: v2018_lwbnd2swbnd1_optdepth_params_type = 'V2018' !< How to specify the parameters for optical depth calculation (Eq.6 in Vallis et al. (2018)). 'V2018', 'USER'
    real(rp) :: v2018_lwbnd2swbnd1_optdepth_sw_params(2)              !< Array of parameters for optical depth calculation with SW band (< 4 μm) (Eq.6 in Vallis et al. (2018)). a_SW, c_SW
    real(rp) :: v2018_lwbnd2swbnd1_optdepth_lw_params(4)              !< Array of parameters for optical depth calculation with LW no-window band (> 4 μm) (Eq.7a in Vallis et al. (2018)). a_LW, b_LW, c_LW, d_LW
    real(rp) :: v2018_lwbnd2swbnd1_optdepth_lw_window_params(4)       !< Array of parameters for optical depth calculation with LW window band (8-14 μm) (Eq.7b in Vallis et al. (2018)). a_win, b_win, c_win, d_win
    real(rp) :: v2018_lwbnd2swbnd1_lw_window_frac                     !< Fraction of for LW window band in the longwave spectrum. 
    integer :: v2018_lwbnd2swbnd1_optdepth_sw_kadd          = 4       !< Number of layers above the model top in the SW optical depth calculation.
    real(rp) :: v2018_lwbnd2swbnd1_optdepth_sw_qv_above_mot = 1e-3_rp !< Optical thickness for shortwave radiation above the model top.
 
    real(rp) :: diffusivityfactor = 1.0_rp !< Diffusivity factor for the two-stream approximation, which is typically set to 1.66 for a plane-parallel atmosphere with isotropic scattering.
                                           !! However, we set it to 1.0 in this idealized radiation scheme, following the formulation in Vallis et al. (2018).
 
    real(rp) :: pco2 = 360.0_rp            !< CO2 concentration in ppmv
 
    namelist / param_atmos_phy_rd_dgm_simple / &
      simple_rd_type,                               &
      !-
      v2018_gray_optdepth_params_type,              &
      v2018_gray_optdepth_lw_params,                &
      !-
      v2018_lwbnd2swbnd1_optdepth_params_type,      &
      v2018_lwbnd2swbnd1_optdepth_sw_params,        &
      v2018_lwbnd2swbnd1_optdepth_lw_params,        &
      v2018_lwbnd2swbnd1_optdepth_lw_window_params, &
      v2018_lwbnd2swbnd1_lw_window_frac,            &
      v2018_lwbnd2swbnd1_optdepth_sw_kadd,          &
      v2018_lwbnd2swbnd1_optdepth_sw_qv_above_mot,  &
      !-
      diffusivityfactor,                            &
      pco2
        
    integer :: ierr
    !--------------------------------------------------------------------
 
    log_newline
    log_info("ATMOS_PHY_RD_dgm_simple_setup",*) 'Setup'
 
    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_atmos_phy_rd_dgm_simple,iostat=ierr)
    if( ierr < 0 ) then !--- missing
       log_info("ATMOS_PHY_RD_dgm_simple_setup",*) 'Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       log_error("ATMOS_PHY_RD_dgm_simple_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_RD_DGM_SIMPLE. Check!'
       call prc_abort
    endif
    log_nml(param_atmos_phy_rd_dgm_simple)
 
    select case(trim(simple_rd_type))
    case('Gray')
      this%optdep_type = optdep_type_v2018_gray
      this%FRAC_LW = [ 1.0_rp, 0.0_rp ]
    case('LWbnd2+SWbnd1')
      this%optdep_type = optdep_type_v2018_lwbnd2swbnd1
    case default
      log_error("ATMOS_PHY_RD_dgm_simple_setup",*) 'SIMPLE_RD_TYPE is invalid. Check!', trim(simple_rd_type)
      call prc_abort
    end select
 
    this%pCO2 = pco2
    this%diffFactor = diffusivityfactor
 
    if ( this%optdep_type == optdep_type_v2018_gray ) then
      select case(trim(v2018_gray_optdepth_params_type))
      case('BO2013')
        this%OPTDEP_VALLIS_EQ5_PARAMS = [ 0.8678_rp, 1.0_rp, 1997.9_rp, 0.0_rp ]
      case('V2018')
        this%OPTDEP_VALLIS_EQ5_PARAMS = [ 0.1627_rp, 1.0_rp, 1997.9_rp, 0.17_rp ]
      case('USER')
        this%OPTDEP_VALLIS_EQ5_PARAMS(:) = v2018_gray_optdepth_lw_params(:)
      case default
        log_error("ATMOS_PHY_RD_dgm_simple_setup",*) 'V2018_GRAY_OPTDEPTH_PARAMS_TYPE is invalid. Check!', trim(v2018_gray_optdepth_params_type)
        call prc_abort
      end select
    end if
 
    if ( this%optdep_type == optdep_type_v2018_lwbnd2swbnd1 ) then
      select case(trim(v2018_lwbnd2swbnd1_optdepth_params_type))
      case('V2018') 
        ! Based on Vallis et al. (2018), the default values for these coefficients were fitted to output from Santa Barbara DISORT Atmospheric Radiative Transfer 60 (SBDART).1
        this%OPTDEP_VALLIS_EQ6a_PARAMS = [ 5.96e-2_rp, 2.9e-3_rp ]
        this%OPTDEP_VALLIS_EQ7a_PARAMS = [ 0.1_rp, 23.8_rp, 254.0_rp, 0.2023_rp ]
        this%OPTDEP_VALLIS_EQ7b_PARAMS = [ 0.215_rp, 1.4711e2_rp, 1.0814e4_rp, 0.0954_rp  ]
        this%FRAC_LW = [ 0.6268_rp, 0.3732_rp ]
        this%OPTDEP_VALLIS_EQ6a_kadd = v2018_lwbnd2swbnd1_optdepth_sw_kadd
        this%OPTDEP_VALLIS_EQ6a_QV_above_MOT = v2018_lwbnd2swbnd1_optdepth_sw_qv_above_mot
      case('USER')
        this%OPTDEP_VALLIS_EQ6a_PARAMS(:) = v2018_lwbnd2swbnd1_optdepth_sw_params(:)
        this%OPTDEP_VALLIS_EQ7a_PARAMS(:) = v2018_lwbnd2swbnd1_optdepth_lw_params(:)
        this%OPTDEP_VALLIS_EQ7b_PARAMS(:) = v2018_lwbnd2swbnd1_optdepth_lw_window_params(:)
        this%FRAC_LW = [ 1.0_rp - v2018_lwbnd2swbnd1_lw_window_frac, v2018_lwbnd2swbnd1_lw_window_frac ]
        this%OPTDEP_VALLIS_EQ6a_kadd = v2018_lwbnd2swbnd1_optdepth_sw_kadd
        this%OPTDEP_VALLIS_EQ6a_QV_above_MOT = v2018_lwbnd2swbnd1_optdepth_sw_qv_above_mot
      case default
        log_error("ATMOS_PHY_RD_dgm_simple_setup",*) 'V2018_LWbnd2SWbnd1_OPTDEPTH_PARAMS_TYPE is invalid. Check!', trim(v2018_lwbnd2swbnd1_optdepth_params_type)
        call prc_abort
      end select
    end if
    return
  end subroutine atm_phy_rd_dgm_simple_init
 
  !> Finalize an object to represent a gray-radiation scheme
!OCL SERIAL
  subroutine atm_phy_rd_dgm_simple_final(this)
    implicit none
    class(atmphyradsimple), intent(inout) :: this
    !--------------------------------------------------------------------
    return
  end subroutine atm_phy_rd_dgm_simple_final
 
  !> Calculate radiative fluxes assuming only absorption and emission of radiation, without scattering.
  !!
!OCL SERIAL
  subroutine atm_phy_rd_dgm_simple_flux( this, &
    flux, flux_top, sflx_up, sflx_dn,                & ! (out)
    solins, pres, temp, dens, qv, sfc_temp, sfc_alb, & ! (in)
    lcmesh, elem3d, lcmesh2d, elem2d )                 ! (in)
    use scale_const, only: &
      pi => const_pi, &
      stb => const_stb
    implicit none
    class(atmphyradsimple), intent(inout) :: this
    class(localmesh3d), intent(in) :: lcmesh
    class(elementbase3d), intent(in) :: elem3d
    class(localmesh2d), intent(in) :: lcmesh2d
    class(elementbase2d), intent(in) :: elem2d
    real(rp), intent(out) :: flux(elem3d%nnode_h1d**2,elem3d%nnode_v,lcmesh%ne2d,lcmesh%nez,2,2)
    real(rp), intent(out) :: flux_top(elem2d%np,lcmesh2d%ne,2,2)
    real(rp), intent(out) :: sflx_up(elem2d%np,lcmesh2d%ne,2)
    real(rp), intent(out) :: sflx_dn(elem2d%np,lcmesh2d%ne,2)
    real(rp), intent(in) :: pres(elem3d%nnode_v,lcmesh%nez,elem3d%nnode_h1d**2,lcmesh%ne2d)
    real(rp), intent(in) :: temp(elem3d%nnode_v,lcmesh%nez,elem3d%nnode_h1d**2,lcmesh%ne2d)
    real(rp), intent(in) :: dens(elem3d%nnode_v,lcmesh%nez,elem3d%nnode_h1d**2,lcmesh%ne2d)
    real(rp), intent(in) :: qv(elem3d%nnode_v,lcmesh%nez,elem3d%nnode_h1d**2,lcmesh%ne2d)
    real(rp), intent(in) :: solins(elem2d%np,lcmesh2d%nea)
    real(rp), intent(in) :: sfc_alb(elem2d%np,lcmesh2d%nea)
    real(rp), intent(in) :: sfc_temp(elem2d%np,lcmesh2d%nea)
 
    integer :: ke, ke_z, ke_h
    integer :: p, p_z, p_h
 
    real(rp) :: dtau_lw(elem3d%nnode_v-1,lcmesh%nez,n_lw_bnd_max)
    real(rp) :: dtau_sw(elem3d%nnode_v-1,lcmesh%nez)
    real(rp) :: dtau_sw_kadd                                       !< Optical thickness for shortwave radiation above the model top. 
                                                                   ! This is used to calculate the downward shortwave flux at the model top.
 
    real(rp) :: trans_lw(elem3d%nnode_v-1,lcmesh%nez,n_lw_bnd_max) !< Transmission across the layer
    real(rp) :: trans_sw(elem3d%nnode_v-1,lcmesh%nez)              !< Transmission across the layer
    real(rp) :: co2(elem3d%nnode_v,lcmesh%nez)
    real(rp) :: src(elem3d%nnode_v-1,lcmesh%nez)                   !< Source term across the layer
    real(rp) :: temp_
    real(rp) :: d
 
    real(rp) :: flux_up_lw(elem3d%nnode_v,lcmesh%nez,n_lw_bnd_max)
    real(rp) :: flux_dn_lw(elem3d%nnode_v,lcmesh%nez,n_lw_bnd_max)
    real(rp) :: flux_up_lw_tot(elem3d%nnode_v,lcmesh%nez)
    real(rp) :: flux_dn_lw_tot(elem3d%nnode_v,lcmesh%nez)
    real(rp) :: flux_up_sw(elem3d%nnode_v,lcmesh%nez)
    real(rp) :: flux_dn_sw(elem3d%nnode_v,lcmesh%nez)
 
    integer :: n_lw_bnd
    integer :: bnd_i
 
    real(rp) :: r_lw(n_lw_bnd_max) !< Fraction of the longwave spectrum in each band
    !--------------------------------------------------------------------
 
    d = this%diffFactor
    r_lw(:) = this%FRAC_LW(:)
 
    select case( this%optdep_type )
    case( optdep_type_v2018_gray )
      n_lw_bnd = 1
    case( optdep_type_v2018_lwbnd2swbnd1 )
      n_lw_bnd = 2
    end select
 
    !$omp parallel do collapse(2) &
    !$omp private( dtau_lw, dtau_sw, dtau_sw_kadd, trans_lw, trans_sw, flux_up_lw, flux_dn_lw, flux_up_sw, flux_dn_sw, &
    !$omp          flux_up_lw_tot, flux_dn_lw_tot, CO2, temp_, Src )
    do ke_h=1, lcmesh%Ne2D
    do p_h=1, elem3d%Nnode_h1D**2
      co2(:,:) = this%pCO2
      call this%calc_optical_thick( dtau_lw, dtau_sw, dtau_sw_kadd,  & ! (out)
        pres(:,:,p_h,ke_h), qv(:,:,p_h,ke_h), co2(:,:),              & ! (in)
        elem3d%Nnode_v, lcmesh%NeZ                                   ) ! (in)
 
      trans_lw(:,:,1:n_lw_bnd) = exp(- d * dtau_lw(:,:,1:n_lw_bnd))
      trans_sw(:,:) = exp(- dtau_sw(:,:))
 
      !- Calculate the source term for longwave radiation
      do ke_z=1, lcmesh%NeZ
      do p_z=1, elem3d%Nnode_v-1
        temp_ = 0.5_rp * ( temp(p_z,ke_z,p_h,ke_h) + temp(p_z+1,ke_z,p_h,ke_h) )
        src(p_z,ke_z) = temp_**4 * stb
      end do
      end do
 
      !- Calculate downward radiative fluxes for longwave and shortwave radiation
 
      flux_dn_lw_tot(:,:) = 0.0_rp
      do bnd_i =1, n_lw_bnd
        flux_dn_lw(elem3d%Nnode_v,lcmesh%NeZ,bnd_i) = 0.0_rp
        do ke_z=lcmesh%NeZ, 1, -1
          do p_z=elem3d%Nnode_v-1, 1, -1
            flux_dn_lw(p_z,ke_z,bnd_i) = flux_dn_lw(p_z+1,ke_z,bnd_i) * trans_lw(p_z,ke_z,bnd_i) &
                                       + r_lw(bnd_i) * src(p_z,ke_z) * ( 1.0_rp - trans_lw(p_z,ke_z,bnd_i) )
          end do
          if ( ke_z > 1 ) then
            flux_dn_lw(elem3d%Nnode_v,ke_z-1,bnd_i) = flux_dn_lw(1,ke_z,bnd_i)
          end if
        end do
 
        flux_dn_lw_tot(:,:) = flux_dn_lw_tot(:,:) + flux_dn_lw(:,:,bnd_i)
      end do
 
      flux_dn_sw(elem3d%Nnode_v,lcmesh%NeZ) = solins(p_h,ke_h) * exp(-dtau_sw_kadd)
      do ke_z=lcmesh%NeZ, 1, -1
      do p_z=elem3d%Nnode_v-1, 1, -1
        flux_dn_sw(p_z,ke_z) = flux_dn_sw(p_z+1,ke_z) * trans_sw(p_z,ke_z)
      end do
        if ( ke_z > 1 ) then
          flux_dn_sw(elem3d%Nnode_v,ke_z-1) = flux_dn_sw(1,ke_z)
        end if      
      end do
 
      !- Calculate upward radiative fluxes for longwave and shortwave radiation
 
      flux_up_lw_tot(:,:) = 0.0_rp
      do bnd_i =1, n_lw_bnd
        flux_up_lw(1,1,bnd_i) = r_lw(bnd_i) * sfc_temp(p_h,ke_h)**4 * stb
        do ke_z=1, lcmesh%NeZ
          do p_z=1, elem3d%Nnode_v-1
            flux_up_lw(p_z+1,ke_z,bnd_i) = flux_up_lw(p_z,ke_z,bnd_i) * trans_lw(p_z,ke_z,bnd_i) &
                                         + r_lw(bnd_i) * src(p_z,ke_z) * ( 1.0_rp - trans_lw(p_z,ke_z,bnd_i) )
          end do
          if ( ke_z < lcmesh%NeZ ) then
            flux_up_lw(1,ke_z+1,bnd_i) = flux_up_lw(elem3d%Nnode_v,ke_z,bnd_i)
          end if
        end do
 
        flux_up_lw_tot(:,:) = flux_up_lw_tot(:,:) + flux_up_lw(:,:,bnd_i)
      end do
 
      ! * Note that the upward flux of shortwave is transparent *
 
      flux_up_sw(1,1) = sfc_alb(p_h,ke_h) * flux_dn_sw(1,1)
      do ke_z=1, lcmesh%NeZ
        do p_z=1, elem3d%Nnode_v-1
          flux_up_sw(p_z+1,ke_z) = flux_up_sw(p_z,ke_z)  ! * trans_sw(p_z,ke_z)
        end do
        if ( ke_z < lcmesh%NeZ ) then
          flux_up_sw(1,ke_z+1) = flux_up_sw(elem3d%Nnode_v,ke_z)
        end if
      end do
 
      !- Store the calculated fluxes into the output arrays
 
      do ke_z=1, lcmesh%NeZ
      do p_z=1, elem3d%Nnode_v
        flux(p_h,p_z,ke_h,ke_z,i_lw,i_up) = flux_up_lw_tot(p_z,ke_z)
        flux(p_h,p_z,ke_h,ke_z,i_lw,i_dn) = flux_dn_lw_tot(p_z,ke_z)
 
        flux(p_h,p_z,ke_h,ke_z,i_sw,i_up) = flux_up_sw(p_z,ke_z)
        flux(p_h,p_z,ke_h,ke_z,i_sw,i_dn) = flux_dn_sw(p_z,ke_z)
      end do
      end do
 
      flux_top(p_h,ke_h,i_lw,i_up) = flux_up_lw_tot(elem3d%Nnode_v,lcmesh%NeZ)
      flux_top(p_h,ke_h,i_lw,i_dn) = flux_dn_lw_tot(elem3d%Nnode_v,lcmesh%NeZ)
      flux_top(p_h,ke_h,i_sw,i_up) = flux_up_sw(elem3d%Nnode_v,lcmesh%NeZ)
      flux_top(p_h,ke_h,i_sw,i_dn) = flux_dn_sw(elem3d%Nnode_v,lcmesh%NeZ)
      sflx_dn(p_h,ke_h,i_lw) = flux_dn_lw_tot(1,1)
      sflx_dn(p_h,ke_h,i_sw) = flux_dn_sw(1,1)
      sflx_up(p_h,ke_h,i_lw) = flux_up_lw_tot(1,1)
      sflx_up(p_h,ke_h,i_sw) = flux_up_sw(1,1)
    end do
    end do
 
    return
  end subroutine atm_phy_rd_dgm_simple_flux
 
!-- private -----
  !> Calculate optical thickness for longwave and shortwave radiation
!OCL SERIAL
  subroutine calc_optical_thick( this, dtau_lw, dtau_sw, dtau_sw_kadd, &
    pres, qv, CO2, Nnode_v, NeZ )
    implicit none
    class(atmphyradsimple), intent(in) :: this
    integer, intent(in) :: nnode_v
    integer, intent(in) :: nez
    real(rp), intent(out) :: dtau_lw(nnode_v-1,nez,n_lw_bnd_max)
    real(rp), intent(out) :: dtau_sw(nnode_v-1,nez)
    real(rp), intent(out) :: dtau_sw_kadd            !< Optical thickness for shortwave radiation above the model top. This is used to calculate the downward shortwave flux at the model top.
    real(rp), intent(in) :: pres(nnode_v,nez)
    real(rp), intent(in) :: qv(nnode_v,nez)
    real(rp), intent(in) :: co2(nnode_v,nez)
 
    integer :: ke_z, p_z
    
    real(rp) :: a, b, c
 
    real(rp) :: a_sw, b_sw, c_sw
    real(rp) :: a_lw, b_lw, c_lw, d_lw
    real(rp) :: a_win, b_win, c_win, d_win
    real(rp) :: tau_sw_
 
    real(rp) :: dsig
 
    real(rp), parameter :: p0 = 1.0e5_rp
    real(rp) :: qv_tmp
    real(rp) :: ln_co2ov360ppm
 
    real(rp) :: dsig_above_mot
    real(rp) :: dtau_sw_kadd_tmp
    integer :: k
    !---------------------------------------------------------
 
    if ( this%optdep_type == optdep_type_v2018_gray ) then
      a = this%OPTDEP_VALLIS_EQ5_PARAMS(1) * this%OPTDEP_VALLIS_EQ5_PARAMS(2)
      b = this%OPTDEP_VALLIS_EQ5_PARAMS(3)
      c = this%OPTDEP_VALLIS_EQ5_PARAMS(4)
      dtau_sw_kadd = 0.0_rp
      
      do ke_z=1, nez
      do p_z=1, nnode_v-1
        qv_tmp = 0.5_rp * ( qv(p_z,ke_z) + qv(p_z+1,ke_z) )
        dsig = max( pres(p_z,ke_z) - pres(p_z+1,ke_z), 0.0_rp ) / p0
        dtau_lw(p_z,ke_z,1) = ( a + b * qv_tmp + c * log( co2(p_z,ke_z) / 360.0_rp ) ) * dsig
        dtau_sw(p_z,ke_z) = 0.0_rp * dsig
      end do
      end do
    
    else if ( this%optdep_type == optdep_type_v2018_lwbnd2swbnd1 ) then
      a_sw = this%OPTDEP_VALLIS_EQ6a_PARAMS(1)
      c_sw = this%OPTDEP_VALLIS_EQ6a_PARAMS(2)
 
      a_lw = this%OPTDEP_VALLIS_EQ7a_PARAMS(1)
      b_lw = this%OPTDEP_VALLIS_EQ7a_PARAMS(2)
      c_lw = this%OPTDEP_VALLIS_EQ7a_PARAMS(3)
      d_lw = this%OPTDEP_VALLIS_EQ7a_PARAMS(4)
 
      a_win = this%OPTDEP_VALLIS_EQ7b_PARAMS(1)
      b_win = this%OPTDEP_VALLIS_EQ7b_PARAMS(2)
      c_win = this%OPTDEP_VALLIS_EQ7b_PARAMS(3)
      d_win = this%OPTDEP_VALLIS_EQ7b_PARAMS(4)
 
      !- Calculate the optical thickness for shortwave radiation above the model top
 
      dsig_above_mot = pres(nnode_v,nez) / p0
      dsig = dsig_above_mot / real(this%OPTDEP_VALLIS_EQ6a_kadd, kind=rp)
      qv_tmp = this%OPTDEP_VALLIS_EQ6a_QV_above_MOT
      ln_co2ov360ppm = log( co2(nnode_v,nez) / 360.0_rp )
 
      tau_sw_ = 0.0_rp
      do k=1, this%OPTDEP_VALLIS_EQ6a_kadd
        b_sw = exp( 1.887e-2_rp / ( tau_sw_ + 9.522e-3_rp ) + 1.603_rp / ( tau_sw_ + 5.194e-1_rp )**2 )
        dtau_sw_kadd_tmp = ( a_sw  + b_sw * qv_tmp +  c_sw * ln_co2ov360ppm ) * dsig
        tau_sw_ = tau_sw_ + dtau_sw_kadd_tmp
      end do
      dtau_sw_kadd = tau_sw_
 
      !-
      do ke_z=nez, 1, -1
        do p_z=nnode_v-1, 1, -1
          qv_tmp = 0.5_rp * ( qv(p_z,ke_z) + qv(p_z+1,ke_z) )
          dsig = max( pres(p_z,ke_z) - pres(p_z+1,ke_z), 0.0_rp ) / p0
          ln_co2ov360ppm = log( co2(p_z,ke_z) / 360.0_rp )
 
          b_sw = exp( 1.887e-2_rp / ( tau_sw_ + 9.522e-3_rp ) + 1.603_rp / ( tau_sw_ + 5.194e-1_rp )**2 )
          dtau_sw(p_z,ke_z)   = ( a_sw  + b_sw * qv_tmp                        +  c_sw * ln_co2ov360ppm ) * dsig
          dtau_lw(p_z,ke_z,1) = ( a_lw  + b_lw * log( c_lw * qv_tmp + 1.0_rp ) +  d_lw * ln_co2ov360ppm ) * dsig
          dtau_lw(p_z,ke_z,2) = ( a_win + qv_tmp * ( b_win + c_win * qv_tmp )  + d_win * ln_co2ov360ppm ) * dsig
 
          tau_sw_ = tau_sw_ + dtau_sw(p_z,ke_z)
        end do
      end do
    end if
 
    return
  end subroutine calc_optical_thick
 
end  module scale_atm_phy_rd_dgm_simple