module rad_Mars_15m

  ! モジュール引用 ; USE statements
  !

  ! 種別型パラメタ
  ! Kind type parameter
  !
  use dc_types, only: DP, &      ! 倍精度実数型. Double precision.
    &                 STRING, &  ! 文字列.       Strings.
    &                 TOKEN      ! キーワード.   Keywords.

  ! 格子点設定
  ! Grid points settings
  !
  use gridset, only: imax, & ! 経度格子点数.
                             ! Number of grid points in longitude
    &                jmax, & ! 緯度格子点数.
                             ! Number of grid points in latitude
    &                kmax    ! 鉛直層数.
                             ! Number of vertical level


  ! 宣言文 ; Declaration statements
  !
  implicit none
  private

  ! 公開手続き
  ! Public procedure
  !
  public :: RadMars15mInit
  public :: RadMars15m


  ! 公開変数
  ! Public variables
  !
  logical, save, public:: rad_Mars_15m_inited = .false.
                              ! 初期設定フラグ.
                              ! Initialization flag.


  ! 非公開変数
  ! Private variables
  !
  real(DP), parameter :: VMRCO2 = 0.95d0
  real(DP), parameter :: AMU    = 1.6605655d-27



  ! Variables for radiation calculation
  !
  ! kg1,kg2,kg3          : maximum number of first, second, third index
  !                      :   for tables of absorption coefficient (cumulative
  !                      : probability function)
  ! kg1(kg1n)            : increament of cumulative probability function
  ! kg2(kg2n)            : pressure (hPa)
  ! kg3(kg3n)            : temperature (K)
  ! dg(kg1n)             : increament of cumulative probability function
  ! lnkg???(...,...,...) : table of absorption coefficient as a function of
  !                      : cumulative probability function, pressure
  !                      : and temperature (m^-1 / (kg m^-2))
  !
  integer , parameter :: kg1n = 16, kg2n = 55, kg3n = 3
  real(DP), save      :: kg1( kg1n ), kg2( kg2n ), kg3( kg3n )
  real(DP), save      :: lnkg( kg1n, kg2n, kg3n )


  ! nl15fn               : maximum number of factors for 15 micron Non-LTE
  !                      : radiative cooling rate calculation
  ! nl15sn               : "reduced" optical depth for 15 micron Non-LTE
  !                      : radiative cooling rate calculation
  ! nl15fa               : parameter for 15 micron Non-LTE
  !                      : radiative cooling rate calculation
  !
  integer , parameter :: nl15fn = 70
  real(DP), save      :: nl15sn( nl15fn )
  real(DP), save      :: nl15fa( nl15fn )


  ! Variables below must have save attribute since these variables are not 
  ! necessarily updated every time steps. 
  !
  ! trans0_res : Transmission from top of the atmosphere to each level.
  ! trans_res  : Transmission between the atmospheric levels
  !
!    real(dp), save :: &
!         trans0_res( km+1, im, jm ), trans_res( km+1, km+1, im, jm )

  real(DP), allocatable, save :: &
    trans_res(:,:,:,:)


  !
  ! 
  !
  real(DP), parameter :: nlte_refp = 1.0d-2


  logical , save :: sw_prep_rv
  data sw_prep_rv /.false./


  integer        , parameter :: lc     = 16


!!$    character(len=lc   ), save      :: ihed( 5 )
!!$    character(len=lc+lc), save      :: titl( 5 )
!!$    data  ihed / 'RAD_P', 'RAD_PH', 'RAD_T', 'RAD_TS', 'RAD_DOD' /
!!$    data  titl / 'RAD_P', 'RAD_PH', 'RAD_T', 'RAD_TS', 'RAD_DOD' /


  real(DP),              save :: rad_time
  real(DP), allocatable, save :: &
    & rad_gp  ( :, :, : ), &
    & rad_gph ( :, :, : ), &
    & rad_gt  ( :, :, : ), &
    & rad_gts ( :, :, : ), &
    & rad_gdod( :, :, : )


  real(DP)             , save :: Rad15mInt

  integer              , save :: nwnl



  !
  ! nwnsl : number of wavenumber sub loop 
  !       : (inner most loop for optimization)
  !
  integer             , save :: nras, nrps
  integer             , save :: nwnsl


!!$  real(DP)    , allocatable, save :: sgmh_f( : ), sgm_f( : )
!!$  real(DP)    , allocatable, save :: gph_f ( :, :, : ), gp_f( :, :, : ), &
!!$    & gth_f( :, :, : )
!!$
!!$  real(DP)    , allocatable, save :: pfh_f ( :, :, : )



  !
  ! *_f( :, ... )             :: variable on fine vertical grids
  !
  ! gvmr_f  ( km, nras+nrps ) :: volume mixing ratio
  ! mmmass_f( km )            :: mean molecular mass
  !
  ! ac( nwnsl, km, nras )   :: absorption coefficient
  !
  ! trans_f( iwnsl, k1, k2 ) :: transmittance between layer interface k1 
  !                          :: and layer midlevel k2
  !
!!$  real(DP)    , allocatable, save :: gvmr_f  ( :, :, :, : )
!!$  real(DP)    , allocatable, save :: mmmass_f( :, :, : )
!!$  real(DP)    , allocatable, save :: ac_f    ( :, :, :, : )
  real(DP)    , allocatable, save :: xyra_Trans (:,:,:,:)
!!$  real(DP)    , allocatable, save :: gdod_f  ( :, :, : )
!!$
!!$  real(DP)    , allocatable, save :: uwflh_f( :, :, : ), dwflh_f( :, :, : )





  real(DP)    , allocatable, save :: &
    & trans_i2i_toa(:,:,:), &          ! f_{1/2}    T_{k+1/2,1/2}
    & trans_i2i_boa(:,:,:), &          ! f_{km+1/2} T_{k+1/2,km+1/2}
    & trans_i2i_s  (:,:,:), &          ! f_{s}      T_{k+1/2,km+1/2}
    & trans_i2m_lli(:,:,:,:), & ! upper layer interface
    & trans_i2m_uli(:,:,:,:) ! lower layer interface


  character(*), parameter:: module_name = 'rad_Mars_15m'
                              ! モジュールの名称.
                              ! Module name
  character(*), parameter:: version = &
    & '$Name:  $' // &
    & '$Id: rad_Mars_15m.f90,v 1.7 2012/11/10 05:00:50 yot Exp $'
                              ! モジュールのバージョン
                              ! Module version


contains

  !**************************************************************************
  ! subroutine rad15m_lowatm_main
  ! calculate radiative heating/cooling rate in CO2 15 micron band
  !**************************************************************************

  subroutine RadMars15m( Time, DelTime, &
    & xyz_Temp, xyr_Press, xyz_Press, xy_SurfTemp, &
    & xyr_DOD067, QeRat, SSA, xy_SurfEmis, &
    & xyr_Rad15mFlux, xyra_DelRad15mFlux &
    & )

    ! メッセージ出力
    ! Message output
    !
    use dc_message, only: MessageNotify

    real(DP)    , intent(in ) :: Time
    real(DP)    , intent(in ) :: DelTime
    real(DP)    , intent(in ) :: xyz_Temp          (0:imax-1, 1:jmax, 1:kmax)
    real(DP)    , intent(in ) :: xyr_Press         (0:imax-1, 1:jmax, 0:kmax)
    real(DP)    , intent(in ) :: xyz_Press         (0:imax-1, 1:jmax, 1:kmax)
    real(DP)    , intent(in ) :: xy_SurfTemp       (0:imax-1, 1:jmax)
    real(DP)    , intent(out) :: xyr_Rad15mFlux    (0:imax-1, 1:jmax, 0:kmax)
    real(DP)    , intent(out) :: xyra_DelRad15mFlux(0:imax-1, 1:jmax, 0:kmax, 0:1)
    real(DP)    , intent(in ) :: xyr_DOD067        (0:imax-1, 1:jmax, 0:kmax)
    real(DP)    , intent(in ) :: QeRat
    real(DP)    , intent(in ) :: SSA
    real(DP)    , intent(in ) :: xy_SurfEmis       (0:imax-1, 1:jmax)


    !
    ! local variables
    !


    ! 実行文 ; Executable statement
    !

    ! 初期化
    ! Initialization
    !
    if ( .not. rad_Mars_15m_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if

    call rad15m_lowatm_newscheme2006( Time, DelTime, &
      & xyz_Temp, xyr_Press, xyz_Press, xy_SurfTemp, &
      & xyr_DOD067, QeRat, SSA, xy_SurfEmis, &
      & xyr_Rad15mFlux, xyra_DelRad15mFlux &
      & )

  end subroutine RadMars15m

  !**************************************************************************
  ! subroutine radiation15m
  ! calculate radiative heating/cooling rate in CO2 15 micron band
  !**************************************************************************

  subroutine rad15m_lowatm_newscheme2006( Time, DelTime, &
    & xyz_Temp, xyr_Press, xyz_Press, xy_SurfTemp, &
    & xyr_DOD067, QeRat, SSA, xy_SurfEmis, &
    & xyr_Rad15mFlux, xyra_DelRad15mFlux &
    & )


    use constants , only : &
      & Grav, &
      & CpDry

    ! メッセージ出力
    ! Message output
    !
    use dc_message, only: MessageNotify

    use ckd_module, only : ckdp


    real(DP)    , intent(in ) :: Time
    real(DP)    , intent(in ) :: DelTime
    real(DP)    , intent(in ) :: xyr_Press(0:imax-1, 1:jmax, 0:kmax)
    real(DP)    , intent(in ) :: xyz_Press(0:imax-1, 1:jmax, 1:kmax)
    real(DP)    , intent(in ) :: xyz_Temp   (0:imax-1, 1:jmax, 1:kmax)
    real(DP)    , intent(in ) :: xy_SurfTemp(0:imax-1, 1:jmax)
    real(DP)    , intent(in ) :: xyr_DOD067(0:imax-1, 1:jmax, 0:kmax)
    real(DP)    , intent(in ) :: QeRat
    real(DP)    , intent(in ) :: SSA
    real(DP)    , intent(in ) :: xy_SurfEmis(0:imax-1, 1:jmax)

    real(DP)    , intent(out) :: xyr_Rad15mFlux (0:imax-1, 1:jmax, 0:kmax)
    real(DP)    , intent(out) :: xyra_DelRad15mFlux(0:imax-1, 1:jmax, 0:kmax, 0:1)

    !
    ! local variables
    !
    real(DP) :: xyr_Temp  (0:imax-1, 1:jmax, 0:kmax)
    real(DP) :: xyz_MMMass(0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyza_VMR(0:imax-1, 1:jmax, 1:kmax, 1:nras+nrps )
    real(DP) :: xyza_AC (0:imax-1, 1:jmax, 1:kmax, 1:nras      )
    real(DP) :: xyr_PF   (0:imax-1, 1:jmax, 0:kmax)
    real(DP) :: xy_SurfPF(0:imax-1, 1:jmax)

    real(DP) :: xy_DPFDT(0:imax-1, 1:jmax)

    real(DP) :: weight_integral
    integer  :: ig, iband

    integer  :: i, j, k, l, m, n
    integer  :: k2

    !
    ! dod      : dust optical depth
    !
    real(DP) :: xyr_DOD(0:imax-1, 1:jmax, 0:kmax)

    !
    ! local variables for pfint
    !

    real(DP) :: MinPress
    real(DP) :: MaxPress

    integer  :: iband_reserve
    real(DP) :: xy_lnPs    (0:imax-1, 1:jmax)
    real(DP) :: xyz_lnPress(0:imax-1, 1:jmax, 1:kmax )
    integer  :: xyz_jj(0:imax-1, 1:jmax, 1:kmax)
    integer  :: xyz_kk(0:imax-1, 1:jmax, 1:kmax)
    integer  :: xy_jj (0:imax-1, 1:jmax)
    integer  :: xy_kk (0:imax-1, 1:jmax)


    ! Surface temperature for calculation of gradient of radiative flux
    real(DP) :: xy_SurfTemp_for_gradcalc(0:imax-1, 1:jmax)
    ! Indices for calculation of gradient of radiative flux
    integer  :: jjs_for_gradcalc(0:imax-1, 1:jmax), kks_for_gradcalc(0:imax-1, 1:jmax)

    real(DP)     :: xyr_PFRat   (0:imax-1, 1:jmax, 0:kmax)
    real(DP)     :: xyz_PFRat   (0:imax-1, 1:jmax, 1:kmax)
    real(DP)     :: xy_SurfPFRat(0:imax-1, 1:jmax)

    logical, save :: FlagCalcTrans

    data FlagCalcTrans / .false. /


    ! 初期化
    ! Initialization
    !
    if ( .not. rad_Mars_15m_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    k = 0
    do j = 1, jmax
      do i = 0, imax-1
!!$        gth(i,j,k) = gt(i,j,k+1)
        xyr_Temp(i,j,k) = &
          &  ( xyz_Temp(i,j,2) - xyz_Temp(i,j,1) ) &
          & / log( xyz_Press(i,j,2) / xyz_Press(i,j,1) ) &
          & * log( xyr_Press(i,j,k) / xyz_Press(i,j,1) ) &
          & + xyz_Temp(i,j,1)
      end do
    end do
    do k = 1, kmax-1
      do j = 1, jmax
        do i = 0, imax-1
          xyr_Temp(i,j,k) = &
            &  ( xyz_Temp(i,j,k+1) - xyz_Temp(i,j,k) ) &
            & / log( xyz_Press(i,j,k+1) / xyz_Press(i,j,k) ) &
            & * log( xyr_Press(i,j,k  ) / xyz_Press(i,j,k) ) &
            & + xyz_Temp(i,j,k)
        end do
      end do
    end do
    k = kmax
    do j = 1, jmax
      do i = 0, imax-1
        xyr_Temp(i,j,k) = xyz_Temp(i,j,k)
      end do
    end do


!!$    do k = 1, km*nvr+1
!!$      do ij = ijs, ije
!!$        gph_f( ij, 1, k ) = sgmh_f( k ) * gph( ij, 1, km+1 )
!!$      end do
!!$    end do
!!$    call calc_lnp( im, jm, km*nvr+1, gph_f , glnph_f , ijs, ije )
!!$
!!$    call increase_vreso_boundary_arr3d( im, jm, km, nvr, gth, gth_f, &
!!$      & "linear", ijs, ije )



    if (  .not. FlagCalcTrans ) then
      if ( Time - dble( int( Time / Rad15mInt ) ) * Rad15mInt < DelTime ) then
        call MessageNotify( 'M', module_name, &
          & 'Transmittance is not saved, but criterion for transmittance calculation is met.' )
      else
        call MessageNotify( 'M', module_name, &
          & 'Transmittance is not saved, and criterion for transmittance calculation ' &
          & // 'is not met. However, transmittance will be calculated.' )
      end if
    end if


    !
    ! Calculation of transmission
    !
    if( ( .not. FlagCalcTrans ) .or. &
      & ( Time - dble( int( Time / Rad15mInt )  ) * Rad15mInt ) < DelTime ) then

      FlagCalcTrans = .true.

!!$      call MessageNotify( 'M', module_name, 'Transmission is calculated.' )

      !
      ! Calculation of "absorption" dust optical depth
      ! This formulation is obtained from Forget et al. [1999].
      !
      do k = 0, kmax
        do j = 1, jmax
          do i = 0, imax-1
            xyr_DOD(i,j,k) = ( 1.0d0 - SSA ) * xyr_DOD067(i,j,k) * QeRat
          end do
        end do
      end do

!!$      call increase_vreso_boundary_arr3d( im, jm, km, nvr, gdod, gdod_f, &
!!$        & "log", ijs, ije )


      !
      ! check pressure
      !
      MinPress = 1.0d100
      MaxPress = 0.0d0
      do j = 1, jmax
        do i = 0, imax-1
          MinPress = min( MinPress, xyz_Press(i,j,kmax) )
          MaxPress = max( MaxPress, xyz_Press(i,j,1   ) )
        end do
      end do
      if( ckdp(1)%lnp(1) > log(MinPress) ) then
        write( 6, * ) 'MARS: pressure is too small.'
        write( 6, * ) MinPress, exp(ckdp(1)%lnp(1))
        stop
      end if
      if( ckdp(1)%lnp(ckdp(1)%nlnp) .lt. log(MaxPress) ) then
        write( 6, * ) 'MARS: pressure is too large.'
        write( 6, * ) MaxPress, exp(ckdp(1)%lnp(ckdp(1)%nlnp))
        stop
      end if


      xyz_MMMass = 43.5d0 * AMU
      do n = 1, nras + nrps
        do k = 1, kmax
          do j = 1, jmax
            do i = 0, imax-1
              xyza_VMR(i,j,k,n) = VMRCO2
            end do
          end do
        end do
      end do


!!$      do n = 1, nras + nrps
!!$        call increase_vreso_b2m_arr3d( im, jm, km, nvr, &
!!$          & gvmrh(:,:,:,n), gvmr_f(:,:,:,n), "log", &
!!$          & ijs, ije )
!!$      end do
!!$      call increase_vreso_b2m_arr3d( im, jm, km, nvr, &
!!$        & mmmassh, mmmass_f, "linear", &
!!$        & ijs, ije )

!!$      call calc_lnp( im, jm, km+1    , gph   , glnph   , ijs, ije )
      call calc_lnp( xyz_Press, xyz_lnPress )
      xy_lnPs(:,:) = log( xyr_Press(:,:,0) )


      !
      ! initialization
      !
      do k = 0, kmax
        do j = 1, jmax
          do i = 0, imax-1
            trans_i2i_toa(i,j,k) = 0.0d0         ! f_{1/2}    T_{k+1/2,1/2}
            trans_i2i_boa(i,j,k) = 0.0d0         ! f_{km+1/2} T_{k+1/2,km+1/2}
            trans_i2i_s  (i,j,k) = 0.0d0         ! f_{s}      T_{k+1/2,km+1/2}
          end do
        end do
      end do
      do k2 = 1, kmax
        do k = 0, kmax
          do j = 1, jmax
            do i = 0, imax-1
              trans_i2m_uli(i,j,k,k2) = 0.0d0
              trans_i2m_lli(i,j,k,k2) = 0.0d0
            end do
          end do
        end do
      end do


      !
      ! loop for wavenumber
      !

      iband_reserve = 0

      do m = 1, nwnl

        call m2ckdpindices( m, ig, iband )


        if( iband .ne. iband_reserve ) then
          call findindices3D( &
            & xyz_Temp, xyz_lnPress, iband, &
            & xyz_jj, xyz_kk &
            & )
          call findindices2D(   &
            & xy_SurfTemp, xy_lnPs,       &
            & iband, xy_jj, xy_kk   &
            & )

          iband_reserve = iband
        end if


        ! IMPORTANT!
        ! This loop for n is confusing.
        ! We have to reconsider about it.
        ! Maybe, the component of ckdp structure has to be reconsidered. 
        ! Now, it cannot include multiple radiatively active species. 
        ! (yot, 2010/09/12)
        !
        do n = 1, nras
          call getlnac_givenindices( &
            & xyz_Temp, xyz_lnPress, xyz_jj, xyz_kk, ig, iband, &
            & xyza_AC(:,:,:,n) &
            & )
        end do
        do n = 1, nras
          xyza_AC(:,:,:,n) = exp( xyza_AC(:,:,:,n) )
        end do

!!$        do n = 1, nras
!!$          call increase_vreso_b2m_arr3d( im, jm, km, nvr, &
!!$            & ach(:,:,:,n), ac_f(:,:,:,n), "log", &
!!$            & ijs, ije )
!!$        end do


        call calc_trans_mp_arr3d(   &
          & nras, nrps, xyr_Press, xyza_VMR, xyz_MMMass, &
          & xyza_AC, xyr_DOD,                     &
          & xyra_Trans                         &
          & )


        call getpfr_givenindices3D( &
          & xyz_Temp, xyz_lnPress, xyz_jj, xyz_kk, ig, iband, &
          & xyz_PFRat &
          & )
        xyr_PFRat(:,:,0) = xyz_PFRat(:,:,1)
        do k = 1, kmax-1
          xyr_PFRat(:,:,k) = ( xyz_PFRat(:,:,k) + xyz_PFRat(:,:,k+1) ) * 0.5_DP
        end do
        xyr_PFRat(:,:,kmax) = xyz_PFRat(:,:,kmax)
        call getpfr_givenindices2D( &
          & xy_SurfTemp, xy_lnPs, xy_jj, xy_kk, ig, iband, &
          & xy_SurfPFRat &
          & )

        do k = 0, kmax
          do j = 1, jmax
            do i = 0, imax-1
              trans_i2i_toa(i,j,k) =       &        ! f_{1/2}    T_{k+1/2,1/2}
                & trans_i2i_toa(i,j,k)     &
                & + xyra_Trans(i,j,k,kmax)      &
                & * xyr_PFRat(i,j,kmax)    &
                & * ckdp(iband)%weight(ig)
              trans_i2i_boa(i,j,k) =       &        ! f_{km+1/2} T_{k+1/2,km+1/2}
                & trans_i2i_boa(i,j,k)     &
                & + xyra_Trans(i,j,k,0)         &
                & * xyr_PFRat(i,j,0)       &
                & * ckdp(iband)%weight(ig)
              trans_i2i_s  (i,j,k) =       &        ! f_{s}      T_{k+1/2,km+1/2}
                & trans_i2i_s  (i,j,k)     &
                & + xyra_Trans(i,j,k,0)         &
                & * xy_SurfPFRat(i,j)      &
                & * ckdp(iband)%weight(ig)
            end do
          end do
        end do

        do k2 = 1, kmax
          do k = 0, kmax
            do j = 1, jmax
              do i = 0, imax-1
                trans_i2m_uli(i,j,k,k2) =                              &
                  & trans_i2m_uli(i,j,k,k2)                            &
                  & + ( xyra_Trans(i,j,k,k2-1) + xyra_Trans(i,j,k,k2) ) * 0.5d0  &
                  & * xyr_PFRat(i,j,k2  )                              &
                  & * ckdp(iband)%weight(ig)
                trans_i2m_lli(i,j,k,k2) =                              &
                  & trans_i2m_lli(i,j,k,k2)                            &
                  & + ( xyra_Trans(i,j,k,k2-1) + xyra_Trans(i,j,k,k2) ) * 0.5d0  &
                  & * xyr_PFRat(i,j,k2-1)                              &
                  & * ckdp(iband)%weight(ig)
              end do
            end do
          end do
        end do


      end do


!!$      call rad15m_rv_put_newscheme2006( time, gt, gph, gp, gts, gdod, ijs, ije )

    else

      if ( trans_i2i_toa(0,1,1) > 1.0d99 ) then
        write( 6, * ) 'transmission function would not be calculated.'
        stop
      end if

    end if


    ! Is this OK?
    iband = 1

    call getpf_arr3d_norat( &
      & xyr_Temp, xy_SurfTemp, iband, &
      & xyr_PF, xy_SurfPF &
      & )


    call calc_rteq_use_meantrans_arr3d( &
      & ( ckdp(iband)%wnbnds(2) - ckdp(iband)%wnbnds(1) ), xy_SurfEmis, &
      & trans_i2i_toa, trans_i2i_boa, trans_i2i_s, &
      & trans_i2m_lli, trans_i2m_uli, &
      & xyr_PF, xy_SurfPF, xyr_Rad15mFlux &
      & )


    do l = 0, 1
      do k = 0, kmax
        do j = 1, jmax
          do i = 0, imax-1
            xyra_DelRad15mFlux(i,j,k,l) = 0.0_DP
          end do
        end do
      end do
    end do


!!$    do k = kmax, 0, -1
!!$      write( 6, * ) gph(0,1,k), gr15mnetflh(0,1,k)
!!$    end do
!!$    stop


!!$      ij = ( ije - ijs + 1 ) / 2
!!$      k  = km + 1
!!$!      write( 6, * ) 'RAD15M : ', gr15mnetflh(ij,1,k), rad_gdr15mnetfldtsh0(ij,1,k-1) * ( gts(ij,1) - rad_gtsbase(ij,1,1) ), emis(ij,1) * 5.67d-8 * gts(ij,1)**4, gts(ij,1)
!!$!      write( 61, * ) gr15mnetflh(ij,1,k-1), rad_gdr15mnetfldtsh0(ij,1,k-1) * ( gts(ij,1) - rad_gtsbase(ij,1,1) ), emis(ij,1) * 5.67d-8 * gts(ij,1)**4, gts(ij,1), gt(ij,1,km), gt(ij,1,km-1), gt(ij,1,km-2), gt(ij,1,km-3), gt(ij,1,km-4), gt(ij,1,km-5)
!!$      write( 6, * ) 'RAD15M : ', gr15mnetflh(ij,1,k), &
!!$        & emis(ij,1) * 5.67d-8 * gts(ij,1)**4, gts(ij,1)
!!$      write( 61, * ) gr15mnetflh(ij,1,k-1), &
!!$        & emis(ij,1) * 5.67d-8 * gts(ij,1)**4, gts(ij,1), &
!!$        & gt(ij,1,km), gt(ij,1,km-1), gt(ij,1,km-2), gt(ij,1,km-3), &
!!$        & gt(ij,1,km-4), gt(ij,1,km-5)
!!$      call flush( 61 )

    !
    ! output variables
    !
!!$    do j = 1, jmax
!!$      do i = 0, imax-1
!!$        goru(i,j) = uwflh_sum(i,j,kmax)
!!$        gord(i,j) = 0.0d0
!!$        gsru(i,j) = uwflh_sum(i,j,0)
!!$        gsrd(i,j) = dwflh_sum(i,j,0)
!!$        gor (i,j) = goru(i,j) - gord(i,j)
!!$        gsr (i,j) = gsru(i,j) - gsrd(i,j)
!!$      end do
!!$    end do



  end subroutine rad15m_lowatm_newscheme2006

  !**************************************************************************

  subroutine calc_trans_mp_arr3d(   &
    & nras, nrps, xyr_Press, xyza_VMR, xyz_MMMass, &
    & xyza_AC, xyr_DOD,                     &
    & xyra_Trans                         &
    & )

    use constants , only : &
      & Grav

    integer , intent(in ) :: nras
    integer , intent(in ) :: nrps
    real(DP), intent(in ) :: xyr_Press(0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(in ) :: xyza_VMR (0:imax-1, 1:jmax, 1:kmax, 1:nras+nrps)
    real(DP), intent(in ) :: xyz_MMMass(0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ) :: xyza_AC(0:imax-1, 1:jmax, 1:kmax, 1:nras)
    real(DP), intent(in ) :: xyr_DOD(0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(out) :: xyra_Trans  (0:imax-1, 1:jmax, 0:kmax, 0:kmax)


    !
    ! local variables
    !
    real(DP)     :: xyz_DelOpDep(0:imax-1, 1:jmax, 1:kmax)
    real(DP)     :: xyz_DelTrans(0:imax-1, 1:jmax, 1:kmax)
    real(DP)     :: xy_Trans(0:imax-1, 1:jmax)
    real(DP), parameter :: DifFac = 1.66_DP


    integer :: k, k2, n
    integer :: ks, ke



    xyra_Trans = 1.0d100


    xyz_DelOpDep = 0.0_DP
    do n = 1, nras
      do k = 1, kmax
        xyz_DelOpDep(:,:,k) = xyz_DelOpDep(:,:,k)                      &
          & + xyza_AC(:,:,k,n) * xyza_VMR(:,:,k,n) / xyz_MMMass(:,:,k) &
          & * ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav
      end do
    end do

    !
    ! add dust optical depth
    !
    do k = 1, kmax
      xyz_DelOpDep(:,:,k) = xyz_DelOpDep(:,:,k) &
        & + xyr_DOD(:,:,k-1) - xyr_DOD(:,:,k)
    end do

    xyz_DelTrans = exp( - xyz_DelOpDep * DifFac )


    !
    ! transmission for "zero thickness" layer ( = 1.0 )
    !
    do ks = 0, kmax
      ke = ks
      xyra_Trans(:,:,ks,ke) = 1.0_DP
    end do

    do ks = 0, kmax
      xy_Trans = 1.0_DP
      do ke = ks+1, kmax
        xy_Trans = xy_Trans * xyz_DelTrans(:,:,ke)
        xyra_Trans(:,:,ks,ke) = xy_Trans
      end do
    end do

    do ks = 0, kmax
      do ke = 0, ks-1
        xyra_Trans(:,:,ks,ke) = xyra_Trans(:,:,ke,ks)
      end do
    end do


  end subroutine calc_trans_mp_arr3d

  !--------------------------------------------------------------------------------------

  subroutine calc_rteq_use_meantrans_arr3d( &
    & dlambda, xy_SurfEmis, &
    & trans_i2i_toa, trans_i2i_boa, trans_i2i_s, &
    & trans_i2m_lli, trans_i2m_uli, &
    & xyr_PF, xy_SurfPF, netflh &
    & )

    ! 物理・数学定数設定
    ! Physical and mathematical constants settings
    !
    use constants0, only: &
      & PI                    ! $ \pi $ .
                              ! 円周率.  Circular constant

    real(DP), intent(in ) :: dlambda
    real(DP), intent(in ) :: xy_SurfEmis(0:imax-1, 1:jmax)
    real(DP), intent(in ) :: trans_i2i_toa(0:imax-1, 1:jmax, 0:kmax)    ! f_{1/2}    T_{k+1/2,1/2}
    real(DP), intent(in ) :: trans_i2i_boa(0:imax-1, 1:jmax, 0:kmax)    ! f_{km+1/2} T_{k+1/2,km+1/2}
    real(DP), intent(in ) :: trans_i2i_s  (0:imax-1, 1:jmax, 0:kmax)    ! f_{s}      T_{k+1/2,km+1/2}
    real(DP), intent(in ) :: trans_i2m_lli(0:imax-1, 1:jmax, 0:kmax, 1:kmax) ! upper layer interface
    real(DP), intent(in ) :: trans_i2m_uli(0:imax-1, 1:jmax, 0:kmax, 1:kmax) ! lower layer interface
    real(DP), intent(in ) :: xyr_PF   (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(in ) :: xy_SurfPF(0:imax-1, 1:jmax)
    real(DP), intent(out) :: netflh(0:imax-1, 1:jmax, 0:kmax)



    !
    ! local variables
    !
    integer :: i, j, k, k2


    do k = 0, kmax
      do j = 1, jmax
        do i = 0, imax-1
          netflh(i,j,k) = 0.0d0
        end do
      end do
    end do

    do k = 0, kmax

      do j = 1, jmax
        do i = 0, imax-1
          netflh(i,j,k) = netflh(i,j,k) &
            & + PI * xy_SurfEmis(i,j) * xy_SurfPF(i,j) * dlambda * trans_i2i_s  (i,j,k) &
            & - PI * xyr_PF(i,j,0   ) * dlambda * trans_i2i_boa(i,j,k) &
            & + PI * xyr_PF(i,j,kmax) * dlambda * trans_i2i_toa(i,j,k)
        end do
      end do

      do k2 = 1, kmax
        do j = 1, jmax
          do i = 0, imax-1
            netflh(i,j,k) = netflh(i,j,k) &
              & - PI * xyr_PF(i,j,k2  ) * dlambda * trans_i2m_uli(i,j,k,k2) &
              & + PI * xyr_PF(i,j,k2-1) * dlambda * trans_i2m_lli(i,j,k,k2)
          end do
        end do
      end do

    end do


  end subroutine calc_rteq_use_meantrans_arr3d

  !**************************************************************************

  subroutine calc_lnp( xyz_Press, xyz_lnPress )

    real(DP), intent(in ) :: xyz_Press  (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(out) :: xyz_lnPress(0:imax-1, 1:jmax, 1:kmax)


    !
    ! local variables
    !


    xyz_lnPress = log( xyz_Press + 1.0d-20 )


  end subroutine calc_lnp

  !**************************************************************************

!!$  subroutine rad15m_2006_rv_read_calctrans( &
!!$    & gt, gph, gp, gts, &
!!$    & gdod, &
!!$    & ijs, ije )
!!$
!!$
!!$    use mars_const, only : vmr_co2, amu
!!$
!!$
!!$    use ckd_module
!!$
!!$
!!$    real(DP), intent(in ) :: gph( im, jm, km+1 ), gp ( im, jm, km )
!!$    real(DP), intent(in ) :: gt( im, jm, km ), gts( im, jm )
!!$    real(DP), intent(in ) :: gdod( im, jm, km+1 )
!!$    integer , intent(in ) :: ijs, ije
!!$
!!$
!!$    !
!!$    ! local variables
!!$    !
!!$    real(DP)                  :: gth( im, jm, km+1 )
!!$    real(DP)                  :: &
!!$      & mmmassh( im, jm, km+1 ), gvmrh( im, jm, km+1, nras+nrps )
!!$    real(DP) :: ach  ( im, jm, km+1, nras )
!!$    real(DP) :: pfh  ( im, jm, km+1 ), pfs  ( im, jm )
!!$    real(DP) :: uwflh( im, jm, km+1 ), dwflh( im, jm, km+1 )
!!$
!!$    real(DP) :: uwflh_sum( im, jm, km+1 ), dwflh_sum( im, jm, km+1 )
!!$
!!$    real(DP) :: pfs_for_gradcalc( im, jm )
!!$    real(DP) :: uwflh_sum_for_gradcalc( im, jm, km+1 ), dwflh_sum_for_gradcalc( im, jm, km+1 )
!!$
!!$    real(DP)                  :: weight_integral
!!$    integer              :: ig, iband
!!$
!!$    integer              :: ij, k, l, m, n
!!$    integer              :: k2
!!$
!!$    real(DP) :: minp, maxp
!!$
!!$    integer :: iband_reserve
!!$    real(DP)     :: glnph  ( im, jm, km+1     )
!!$    real(DP)     :: glnph_f( im, jm, km*nvr+1 )
!!$    real(DP)     :: gts3d1 ( im, jm, 1        )
!!$    integer :: &
!!$      & jjh   ( im, jm, km+1     ), kkh   ( im, jm, km+1     ), &
!!$      & jjh_f ( im, jm, km*nvr+1 ), kkh_f ( im, jm, km*nvr+1 ), &
!!$      & jjs3d1( im, jm, 1        ), kks3d1( im, jm, 1        ), &
!!$      & jjs   ( im, jm )          , kks   ( im, jm )
!!$
!!$
!!$    ! Surface temperature for calculation of gradient of radiative flux
!!$    real(DP) :: gts_for_gradcalc( im, jm )
!!$    ! Indices for calculation of gradient of radiative flux
!!$    integer  :: jjs_for_gradcalc( im, jm ), kks_for_gradcalc( im, jm )
!!$
!!$
!!$    real(DP)     :: glnps3d1( im, jm, 1 )
!!$
!!$    real(DP)     :: pfrh_f( im, jm, km*nvr+1 )
!!$    real(DP)     :: pfr3d1( im, jm, 1        )
!!$
!!$
!!$    k = 1
!!$    do ij = ijs, ije
!!$      gth( ij, 1, k ) = gt( ij, 1, k )
!!$    end do
!!$    do k = 1+1, km+1-1
!!$      do ij = ijs, ije
!!$        gth( ij, 1, k ) &
!!$          = ( gt( ij, 1, k ) - gt( ij, 1, k-1 ) ) &
!!$          / log( gp( ij, 1, k ) / gp( ij, 1, k-1 ) ) &
!!$          * log( gph( ij, 1, k ) / gp( ij, 1, k-1 ) ) &
!!$          + gt( ij, 1, k-1 )
!!$      end do
!!$    end do
!!$    k = km+1
!!$    do ij = ijs, ije
!!$      gth( ij, 1, k ) &
!!$        = ( gt( ij, 1, km ) - gt( ij, 1, km-1 ) ) &
!!$        / log( gp( ij, 1, km ) / gp( ij, 1, km-1 ) ) &
!!$        * log( gph( ij, 1, km+1 ) / gp( ij, 1, km-1 ) ) &
!!$        + gt( ij, 1, km-1 )
!!$    end do
!!$
!!$
!!$    do k = 1, km*nvr+1
!!$      do ij = ijs, ije
!!$        gph_f( ij, 1, k ) = sgmh_f( k ) * gph( ij, 1, km+1 )
!!$      end do
!!$    end do
!!$    call calc_lnp( im, jm, km*nvr+1, gph_f , glnph_f , ijs, ije )
!!$
!!$    call increase_vreso_boundary_arr3d( im, jm, km, nvr, gth, gth_f, &
!!$      & "linear", ijs, ije )
!!$
!!$
!!$
!!$    !
!!$    ! Set interval for radiation/transmission calculation
!!$    !
!!$!      radint = hourm * mins
!!$
!!$
!!$    !
!!$    ! Calculation of transmission
!!$    !
!!$!      if( ( ( time - dble( int( time / rad15mint )  ) * rad15mint ) .lt. dt ) &
!!$!        .or. ( rad_gtsbase( ijs, 1, 1 ) .gt. 1.0d99 ) ) then
!!$!
!!$!
!!$!         write( 6, * ) '########## rad15m in if'
!!$
!!$
!!$    !
!!$    ! Calculation of "absorption" dust optical depth
!!$    ! This formulation is obtained from Forget et al. [1999].
!!$    !
!!$!         do k = 1, km+1
!!$!            do ij = ijs, ije
!!$!               gdod( ij, 1, k ) = ( 1.0d0 - ssa ) * dod067( ij, 1, k ) * qerat
!!$!            end do
!!$!         end do
!!$
!!$    call increase_vreso_boundary_arr3d( im, jm, km, nvr, gdod, gdod_f, &
!!$      & "log", ijs, ije )
!!$
!!$
!!$    !
!!$    ! check pressure
!!$    !
!!$    minp = 1.0d100
!!$    maxp = 0.0d0
!!$    do ij = ijs, ije
!!$      minp = min( minp, gp( ij, 1, 2  ) )
!!$      maxp = max( maxp, gp( ij, 1, km ) )
!!$    end do
!!$    if( ckdp(1)%lnp(1) .gt. log(minp) ) then
!!$      write( 6, * ) 'MARS: pressure is too small.'
!!$      write( 6, * ) minp, exp(ckdp(1)%lnp(1))
!!$      stop
!!$    end if
!!$    if( ckdp(1)%lnp(ckdp(1)%nlnp) .lt. log(maxp) ) then
!!$      write( 6, * ) 'MARS: pressure is too large.'
!!$      write( 6, * ) maxp, exp(ckdp(1)%lnp(ckdp(1)%nlnp))
!!$      stop
!!$    end if
!!$
!!$
!!$    do k = 1, km+1
!!$      do ij = ijs, ije
!!$        mmmassh( ij, 1, k ) = 43.5d0 * amu
!!$      end do
!!$    end do
!!$    do n = 1, nras + nrps
!!$      do k = 1, km+1
!!$        do ij = ijs, ije
!!$          gvmrh( ij, 1, k, n ) = vmr_co2
!!$        end do
!!$      end do
!!$    end do
!!$
!!$
!!$!    do n = 1, nras + nrps
!!$!      call increase_vreso_b2m_arr3d( im, jm, km, nvr, &
!!$!        gvmrh(:,:,:,n), gvmr_f(:,:,:,n), "log", &
!!$!        ijs, ije )
!!$!    end do
!!$!    call increase_vreso_b2m_arr3d( im, jm, km, nvr, &
!!$!      mmmassh, mmmass_f, "linear", &
!!$!      ijs, ije )
!!$
!!$!    call calc_lnp( im, jm, km+1    , gph   , glnph   , ijs, ije )
!!$    call calc_lnph( gph, glnph )
!!$
!!$
!!$    !
!!$    ! initialization
!!$    !
!!$    do k = 0, kmax
!!$      do j = 1, jmax
!!$        do i = 0, imax-1
!!$          trans_i2i_toa(i,j,k) = 0.0d0         ! f_{1/2}    T_{k+1/2,1/2}
!!$          trans_i2i_boa(i,j,k) = 0.0d0         ! f_{km+1/2} T_{k+1/2,km+1/2}
!!$          trans_i2i_s  (i,j,k) = 0.0d0         ! f_{s}      T_{k+1/2,km+1/2}
!!$        end do
!!$      end do
!!$    end do
!!$    do k2 = 1, kmax
!!$      do k = 0, kmax
!!$        do j = 1, jmax
!!$          do i = 0, imax-1
!!$            trans_i2m_uli_f(i,j,k,k2) = 0.0d0
!!$            trans_i2m_lli_f(i,j,k,k2) = 0.0d0
!!$          end do
!!$        end do
!!$      end do
!!$    end do
!!$
!!$
!!$    !
!!$    ! loop for wavenumber
!!$    !
!!$
!!$    iband_reserve = 0
!!$
!!$    do m = 1, nwnl
!!$
!!$      call m2ckdpindices( m, ig, iband )
!!$
!!$
!!$      if( iband .ne. iband_reserve ) then
!!$        call findindices( im, jm, km+1, gth, glnph, iband, &
!!$          jjh, kkh, ijs, ije )
!!$
!!$
!!$        call findindices( im, jm, km*nvr+1, &
!!$          gth_f , glnph_f , &
!!$          iband, jjh_f , kkh_f , ijs, ije )
!!$        do ij = ijs, ije
!!$          gts3d1( ij, 1, 1 ) = gts( ij, 1 )
!!$        end do
!!$        call findindices( im, jm, 1       , &
!!$          gts3d1, glnph_f(:,:,km*nvr+1), &
!!$          iband, jjs3d1, kks3d1, ijs, ije )
!!$        do ij = ijs, ije
!!$          jjs( ij, 1 ) = jjs3d1( ij, 1, 1 )
!!$          kks( ij, 1 ) = kks3d1( ij, 1, 1 )
!!$        end do
!!$
!!$
!!$        iband_reserve = iband
!!$      end if
!!$
!!$
!!$      do n = 1, nras
!!$        call getlnac_givenindices( im, jm, km+1, gth, glnph, jjh, kkh, &
!!$          & ach(:,:,:,n), ig, iband, ijs, ije )
!!$      end do
!!$      do n = 1, nras
!!$        do k = 1, km+1
!!$          do ij = ijs, ije
!!$            ach( ij, 1, k, n ) = exp( ach( ij, 1, k, n ) )
!!$          end do
!!$        end do
!!$      end do
!!$
!!$      do n = 1, nras
!!$        call increase_vreso_b2m_arr3d( im, jm, km, nvr, &
!!$          ach(:,:,:,n), ac_f(:,:,:,n), "log", &
!!$          ijs, ije )
!!$      end do
!!$
!!$
!!$      call calc_trans_mp_arr3d( nras, nrps, im, jm, km*nvr, &
!!$        gph_f, gvmr_f, mmmass_f, &
!!$        ac_f, gdod_f, trans_f, ijs, ije )
!!$
!!$
!!$      do ij = ijs, ije
!!$        gts3d1  ( ij, 1, 1 ) = gts    ( ij, 1 )
!!$        glnps3d1( ij, 1, 1 ) = glnph_f( ij, 1, km*nvr+1 )
!!$        jjs3d1  ( ij, 1, 1 ) = jjs    ( ij, 1 )
!!$        kks3d1  ( ij, 1, 1 ) = kks    ( ij, 1 )
!!$      end do
!!$
!!$      call getpfr_givenindices( im, jm, km*nvr+1, gth_f , glnph_f , jjh_f , kkh_f , pfrh_f, ig, iband, ijs, ije )
!!$      call getpfr_givenindices( im, jm, 1       , gts3d1, glnps3d1, jjs3d1, kks3d1, pfr3d1, ig, iband, ijs, ije )
!!$
!!$
!!$      do k = 1, km*nvr+1
!!$        do ij = ijs, ije
!!$          trans_i2i_toa_f(ij,1,k) = &        ! f_{1/2}    T_{k+1/2,1/2}
!!$            & trans_i2i_toa_f(ij,1,k) &
!!$            & + trans_f(ij,1,k,1   ) &
!!$            & * pfrh_f(ij,1,1  ) &
!!$            & * ckdp(iband)%weight(ig)
!!$          trans_i2i_boa_f(ij,1,k) = &        ! f_{km+1/2} T_{k+1/2,km+1/2}
!!$            & trans_i2i_boa_f(ij,1,k) &
!!$            & + trans_f(ij,1,k,km+1) &
!!$            & * pfrh_f(ij,1,km+1) &
!!$            & * ckdp(iband)%weight(ig)
!!$          trans_i2i_s_f  (ij,1,k) = &        ! f_{s}      T_{k+1/2,km+1/2}
!!$            & trans_i2i_s_f  (ij,1,k) &
!!$            & + trans_f(ij,1,k,km+1) &
!!$            & * pfr3d1(ij,1,1) &
!!$            & * ckdp(iband)%weight(ig)
!!$        end do
!!$      end do
!!$
!!$      do k2 = 1, km*nvr
!!$        do k = 1, km*nvr+1
!!$          do ij = ijs, ije
!!$            trans_i2m_uli_f(ij,1,k,k2) = &
!!$              & trans_i2m_uli_f(ij,1,k,k2) &
!!$              & + ( trans_f(ij,1,k,k2) + trans_f(ij,1,k,k2+1) ) * 0.5d0 &
!!$              & * pfrh_f(ij,1,k2  ) &
!!$              & * ckdp(iband)%weight(ig)
!!$            trans_i2m_lli_f(ij,1,k,k2) = &
!!$              & trans_i2m_lli_f(ij,1,k,k2) &
!!$              & + ( trans_f(ij,1,k,k2) + trans_f(ij,1,k,k2+1) ) * 0.5d0 &
!!$              & * pfrh_f(ij,1,k2+1) &
!!$              & * ckdp(iband)%weight(ig)
!!$          end do
!!$        end do
!!$      end do
!!$
!!$
!!$
!!$    end do
!!$
!!$
!!$
!!$
!!$
!!$!      end if
!!$
!!$
!!$  end subroutine rad15m_2006_rv_read_calctrans

  !**************************************************************************

!!$  subroutine rad15m_readnlte15mfac( fn )
!!$
!!$
!!$    interface
!!$      subroutine findfu( fn, ios, fu, mode )
!!$        use matype
!!$        implicit none
!!$        character(len=*), intent(in )           :: fn
!!$        integer    , intent(out)           :: ios, fu
!!$        character(len=*), intent(in ), optional :: mode
!!$      end subroutine findfu
!!$    end interface
!!$
!!$
!!$    character(len=*), intent(in) :: fn
!!$
!!$
!!$    !
!!$    ! local variables
!!$    !
!!$    character(len=128) :: tmpl
!!$    integer       :: ios, fu
!!$    integer       :: i
!!$
!!$
!!$    call findfu( fn, ios, fu )
!!$    if( ios /= 0 ) then
!!$      write( 6, * ) 'STOP in parse_ctl: ', ios
!!$      stop
!!$    endif
!!$    open( fu, file = fn, status='unknown' )
!!$    read( fu, '(a)' ) tmpl
!!$    do i = 1, nl15fn
!!$      read( fu, * ) nl15sn( i ), nl15fa( i )
!!$    enddo
!!$    close( fu )
!!$
!!$
!!$  end subroutine rad15m_readnlte15mfac

  !**************************************************************************

  !--------------------------------------------------------------------------------------

  subroutine m2ckdpindices( m, ig, iband )

    use ckd_module, only : ckdp, nband

    integer, intent(in ) :: m
    integer, intent(out) :: ig
    integer, intent(out) :: iband


    !
    ! local variables
    !
    integer :: num


    ! The comments below will be removed.


    num = 0
    do iband = 1, nband
      if( num + ckdp( iband ) % ng .ge. m ) exit
      num = num + ckdp( iband ) % ng
    end do
    if( iband > nband ) then
      write( 6, * ) 'Unexpected m'
      write( 6, * ) m
      stop
    end if
    ig = m - num
    if( ig > ckdp( iband ) % ng ) then
      write( 6, * ) 'Unexpected ig'
      write( 6, * ) iband, ig
      stop
    end if


  end subroutine m2ckdpindices

  !--------------------------------------------------------------------------------------

  subroutine getlnac_givenindices( &
    & xyz_Temp, xyz_lnPress, xyz_jj, xyz_kk, ig, iband, &
    & xyz_AC &
    & )

    use ckd_module, only : ckdp

    real(DP), intent(in ) :: xyz_Temp   (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ) :: xyz_lnPress(0:imax-1, 1:jmax, 1:kmax)
    integer , intent(in ) :: xyz_jj  (0:imax-1, 1:jmax, 1:kmax)
    integer , intent(in ) :: xyz_kk  (0:imax-1, 1:jmax, 1:kmax)
    integer , intent(in ) :: ig
    integer , intent(in ) :: iband
    real(DP), intent(out) :: xyz_AC  (0:imax-1, 1:jmax, 1:kmax)


    !
    ! local variables
    !
    real(DP) :: lnac1, lnac2
    integer  :: i, j, k


    do k = 1, kmax
      do j = 1, jmax
        do i = 0, imax-1

          lnac1 &
            & = ( ckdp(iband)%lnac( ig, xyz_jj(i,j,k)  , xyz_kk(i,j,k)+1 )      &
            &   - ckdp(iband)%lnac( ig, xyz_jj(i,j,k)  , xyz_kk(i,j,k)   ) )    &
            &  / ( ckdp(iband)%t( xyz_kk(i,j,k)+1 )                          &
            &    - ckdp(iband)%t( xyz_kk(i,j,k)   ) )                        &
            & * ( xyz_Temp(i,j,k) - ckdp( iband ) % t( xyz_kk(i,j,k) ) )          &
            & + ckdp(iband)%lnac( ig, xyz_jj(i,j,k)  , xyz_kk(i,j,k)   )
          lnac2 &
            & = ( ckdp(iband)%lnac( ig, xyz_jj(i,j,k)+1, xyz_kk(i,j,k)+1 )      &
            &   - ckdp(iband)%lnac( ig, xyz_jj(i,j,k)+1, xyz_kk(i,j,k)   ) )    &
            &  / ( ckdp(iband)%t( xyz_kk(i,j,k)+1 )                          &
            &    - ckdp(iband)%t( xyz_kk(i,j,k)   ) )                        &
            & * ( xyz_Temp(i,j,k) - ckdp( iband ) % t( xyz_kk(i,j,k) ) )          &
            & + ckdp(iband)%lnac( ig, xyz_jj(i,j,k)+1, xyz_kk(i,j,k)   )

          xyz_AC(i,j,k) &
            & = ( lnac2 - lnac1 ) &
            & / ( ckdp( iband ) % lnp( xyz_jj(i,j,k)+1 ) &
            & - ckdp( iband ) % lnp( xyz_jj(i,j,k)   ) ) &
            & * ( xyz_lnPress(i,j,k) - ckdp( iband ) % lnp( xyz_jj(i,j,k) ) ) &
            & + lnac1
        end do
      end do
    end do


  end subroutine getlnac_givenindices

  !--------------------------------------------------------------------------------------
!!$
!!$  subroutine calc_lnph( gph, glnph )
!!$
!!$    real(DP), intent(in ) :: gph  (0:imax-1, 1:jmax, 1:kmax)
!!$    real(DP), intent(out) :: glnph(0:imax-1, 1:jmax, 0:kmax)
!!$
!!$
!!$    !
!!$    ! local variables
!!$    !
!!$    integer :: i, j, k
!!$
!!$
!!$    do k = 0, kmax
!!$      do j = 1, jmax
!!$        do i = 0, imax-1
!!$          glnph(i,j,k) = log( gph(i,j,k) + 1.0d-20 )
!!$        end do
!!$      end do
!!$    end do
!!$
!!$
!!$  end subroutine calc_lnph
!!$
  !--------------------------------------------------------------------------------------

  subroutine findindices( &
    & ks, ke, xyz_Temp, xyz_lnPress, iband, &
    & xyz_jj, xyz_kk &
    & )

    use ckd_module, only : ckdp

    integer , intent(in ) :: ks
    integer , intent(in ) :: ke
    real(DP), intent(in ) :: xyz_Temp   (0:imax-1, 1:jmax, ks:ke)
    real(DP), intent(in ) :: xyz_lnPress(0:imax-1, 1:jmax, ks:ke)
    integer , intent(in ) :: iband
    integer , intent(out) :: xyz_jj(0:imax-1, 1:jmax, ks:ke)
    integer , intent(out) :: xyz_kk(0:imax-1, 1:jmax, ks:ke)


    !
    ! local variables
    !
    integer :: i, j, k, l


    do k = ks, ke
      do j = 1, jmax
        do i = 0, imax-1
          xyz_kk(i,j,k) = 1
        end do
      end do
    end do

    do l = 1+1, ckdp( iband ) % nt - 1
      do k = ks, ke
        do j = 1, jmax
          do i = 0, imax-1
            if( ckdp( iband ) % t( l ) .le. xyz_Temp(i,j,k) ) xyz_kk(i,j,k) = l
          end do
        end do
      end do
    end do

    do k = ks, ke
      do j = 1, jmax
        do i = 0, imax-1
          xyz_jj(i,j,k) = 1
        end do
      end do
    end do
    do l = 1+1, ckdp( iband ) % nlnp - 1
      do k = ks, ke
        do j = 1, jmax
          do i = 0, imax-1
            if( ckdp( iband ) % lnp( l ) <= xyz_lnPress(i,j,k) ) xyz_jj(i,j,k) = l
          end do
        end do
      end do
    end do


  end subroutine findindices

  !--------------------------------------------------------------------------------------

  subroutine findindices3D( &
    & xyz_Temp, xyz_lnPress, iband, &
    & xyz_jj, xyz_kk &
    & )

    real(DP), intent(in ) :: xyz_Temp   (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ) :: xyz_lnPress(0:imax-1, 1:jmax, 1:kmax)
    integer , intent(in ) :: iband
    integer , intent(out) :: xyz_jj(0:imax-1, 1:jmax, 1:kmax)
    integer , intent(out) :: xyz_kk(0:imax-1, 1:jmax, 1:kmax)


    !
    ! local variables
    !


    call findindices(             &
      & 1, kmax, xyz_Temp, xyz_lnPress, iband, &
      & xyz_jj, xyz_kk                    &
      & )


  end subroutine findindices3D

  !--------------------------------------------------------------------------------------

  subroutine findindices2D( &
    & xy_Temp, xy_lnPress, iband, &
    & xy_jj, xy_kk &
    & )

    real(DP), intent(in ) :: xy_Temp   (0:imax-1, 1:jmax)
    real(DP), intent(in ) :: xy_lnPress(0:imax-1, 1:jmax)
    integer , intent(in ) :: iband
    integer , intent(out) :: xy_jj(0:imax-1, 1:jmax)
    integer , intent(out) :: xy_kk(0:imax-1, 1:jmax)


    !
    ! local variables
    !
    real(DP) :: xyz_Temp   (0:imax-1, 1:jmax, 1:1)
    real(DP) :: xyz_lnPress(0:imax-1, 1:jmax, 1:1)
    integer  :: xyz_jj(0:imax-1, 1:jmax, 1:1)
    integer  :: xyz_kk(0:imax-1, 1:jmax, 1:1)


    xyz_Temp   (:,:,1) = xy_Temp
    xyz_lnPress(:,:,1) = xy_lnPress

    call findindices(              &
      & 1, 1, xyz_Temp, xyz_lnPress, iband, &
      & xyz_jj, xyz_kk                 &
      & )

    xy_jj = xyz_jj(:,:,1)
    xy_kk = xyz_kk(:,:,1)


  end subroutine findindices2D

  !--------------------------------------------------------------------------------------

  subroutine getpf_arr3d_norat( &
    & xyz_Temp, xy_SurfTemp, iband, &
    & xyr_PF, xy_SurfPF &
    & )

    use planck_func, only : Integ_PF_GQ_Array3D, Integ_PF_GQ_Array2D

    use ckd_module, only : ckdp

    real(DP), intent(in ) :: xyz_Temp   (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(in ) :: xy_SurfTemp(0:imax-1, 1:jmax)
    integer , intent(in ) :: iband
    real(DP), intent(out) :: xyr_PF   (0:imax-1, 1:jmax, 0:kmax)
    real(DP), intent(out) :: xy_SurfPF(0:imax-1, 1:jmax)


    !
    ! local variables
    !
    integer :: ncp_pfint
    integer :: i, j, k


    ncp_pfint = 5

    call Integ_PF_GQ_Array3D( &
      & ckdp(iband)%wnbnds(1), ckdp(iband)%wnbnds(2), ncp_pfint, &
      & 0, imax-1, 1, jmax, 0, kmax, &
      & xyz_Temp, &
      & xyr_PF &
      & )
    call Integ_PF_GQ_Array2D( &
      & ckdp(iband)%wnbnds(1), ckdp(iband)%wnbnds(2), ncp_pfint, &
      & 0, imax-1, 1, jmax, &
      & xy_SurfTemp, &
      & xy_SurfPF &
      & )

    do k = 0, kmax
      do j = 1, jmax
        do i = 0, imax-1
          xyr_PF(i,j,k) = xyr_PF(i,j,k) / ( ckdp(iband)%wnbnds(2) - ckdp(iband)%wnbnds(1) )
        end do
      end do
    end do
    do j = 1, jmax
      do i = 0, imax-1
        xy_SurfPF(i,j) = xy_SurfPF(i,j) / ( ckdp(iband)%wnbnds(2) - ckdp(iband)%wnbnds(1) )
      end do
    end do


  end subroutine getpf_arr3d_norat

  !--------------------------------------------------------------------------------------

  subroutine getpfr_givenindices( &
    & ks, ke, &
    & xyz_Temp, xyz_lnPress, xyz_jj, xyz_kk, ig, iband, &
    & xyz_PFRat &
    & )

    use ckd_module, only: ckdp

    integer , intent(in ) :: ks
    integer , intent(in ) :: ke
    real(DP), intent(in ) :: xyz_Temp   (0:imax-1, 1:jmax, ks:ke)
    real(DP), intent(in ) :: xyz_lnPress(0:imax-1, 1:jmax, ks:ke)
    integer , intent(in ) :: xyz_jj  (0:imax-1, 1:jmax, ks:ke)
    integer , intent(in ) :: xyz_kk  (0:imax-1, 1:jmax, ks:ke)
    integer , intent(in ) :: ig, iband
    real(DP), intent(out) :: xyz_PFRat(0:imax-1, 1:jmax, ks:ke)


    !
    ! local variables
    !
    real(DP) :: pfr1, pfr2
    integer  :: i, j, k, l


    do k = ks, ke
      do j = 1, jmax
        do i = 0, imax-1

          pfr1 = &
            &   ( ckdp(iband)%pfr( ig, xyz_jj(i,j,k)  , xyz_kk(i,j,k)+1 )     &
            &   - ckdp(iband)%pfr( ig, xyz_jj(i,j,k)  , xyz_kk(i,j,k)   ) )   &
            & / ( ckdp(iband)%t( xyz_kk(i,j,k)+1 )                         &
            &   - ckdp(iband)%t( xyz_kk(i,j,k)   ) )                       &
            & * ( xyz_Temp(i,j,k) - ckdp( iband ) % t( xyz_kk(i,j,k) ) )         &
            & + ckdp(iband)%pfr( ig, xyz_jj(i,j,k)  , xyz_kk(i,j,k)   )
          pfr2 = &
            &   ( ckdp(iband)%pfr( ig, xyz_jj(i,j,k)+1, xyz_kk(i,j,k)+1 )      &
            &   - ckdp(iband)%pfr( ig, xyz_jj(i,j,k)+1, xyz_kk(i,j,k)   ) )    &
            & / ( ckdp(iband)%t( xyz_kk(i,j,k)+1 )                         &
            &   - ckdp(iband)%t( xyz_kk(i,j,k)   ) )                       &
            & * ( xyz_Temp(i,j,k) - ckdp( iband ) % t( xyz_kk(i,j,k) ) )         &
            & + ckdp(iband)%pfr( ig, xyz_jj(i,j,k)+1, xyz_kk(i,j,k)   )


          xyz_PFRat(i,j,k) = &
            &   ( pfr2 - pfr1 ) &
            & / ( ckdp( iband ) % lnp( xyz_jj(i,j,k)+1 ) &
            &   - ckdp( iband ) % lnp( xyz_jj(i,j,k)   ) ) &
            & * ( xyz_lnPress(i,j,k) - ckdp( iband ) % lnp( xyz_jj(i,j,k) ) ) &
            & + pfr1
        end do
      end do
    end do


  end subroutine getpfr_givenindices

  !--------------------------------------------------------------------------------------

  subroutine getpfr_givenindices3D( &
    & xyz_Temp, xyz_lnPress, xyz_jj, xyz_kk, ig, iband, &
    & xyz_PFRat &
    & )

    use ckd_module, only: ckdp

    real(DP), intent(in ) :: xyz_Temp   (0:imax-1, 1:jmax, 1:kmax)
    real(DP), intent(in ) :: xyz_lnPress(0:imax-1, 1:jmax, 1:kmax)
    integer , intent(in ) :: xyz_jj  (0:imax-1, 1:jmax, 1:kmax)
    integer , intent(in ) :: xyz_kk  (0:imax-1, 1:jmax, 1:kmax)
    integer , intent(in ) :: ig, iband
    real(DP), intent(out) :: xyz_PFRat(0:imax-1, 1:jmax, 1:kmax)


    !
    ! local variables
    !

    call getpfr_givenindices( &
      & 1, kmax, &
      & xyz_Temp, xyz_lnPress, xyz_jj, xyz_kk, ig, iband, &
      & xyz_PFRat &
      & )


  end subroutine getpfr_givenindices3D

  !--------------------------------------------------------------------------------------

  subroutine getpfr_givenindices2D( &
    & xy_Temp, xy_lnPress, xy_jj, xy_kk, ig, iband, &
    & xy_PFRat &
    & )

    real(DP), intent(in ) :: xy_Temp   (0:imax-1, 1:jmax)
    real(DP), intent(in ) :: xy_lnPress(0:imax-1, 1:jmax)
    integer , intent(in ) :: xy_jj  (0:imax-1, 1:jmax)
    integer , intent(in ) :: xy_kk  (0:imax-1, 1:jmax)
    integer , intent(in ) :: ig, iband
    real(DP), intent(out) :: xy_PFRat (0:imax-1, 1:jmax)


    !
    ! local variables
    !
    real(DP) :: xyz_Temp   (0:imax-1, 1:jmax, 1:1)
    real(DP) :: xyz_lnPress(0:imax-1, 1:jmax, 1:1)
    integer  :: xyz_jj  (0:imax-1, 1:jmax, 1:1)
    integer  :: xyz_kk  (0:imax-1, 1:jmax, 1:1)
    real(DP) :: xyz_PFRat(0:imax-1, 1:jmax, 1:1)


    xyz_Temp   (:,:,1) = xy_Temp
    xyz_lnPress(:,:,1) = xy_lnPress
    xyz_jj  (:,:,1) = xy_jj
    xyz_kk  (:,:,1) = xy_kk

    call getpfr_givenindices( &
      & 1, 1, &
      & xyz_Temp, xyz_lnPress, xyz_jj, xyz_kk, ig, iband, &
      & xyz_PFRat &
      & )

    xy_PFRat(:,:) = xyz_PFRat(:,:,1)


  end subroutine getpfr_givenindices2D

  !--------------------------------------------------------------------------------------

!!$  function m2ib( m ) result( iband )
!!$
!!$    use matype
!!$    use ckd_module
!!$
!!$    integer(i4b), intent(in ) :: m
!!$    integer(i4b)              :: iband
!!$
!!$
!!$    !
!!$    ! local variables
!!$    !
!!$    integer(i4b) :: ig
!!$    integer(i4b) :: num
!!$
!!$
!!$    num = 0
!!$    do iband = 1, nband
!!$      if( num + ckdp( iband ) % ng .ge. m ) exit
!!$      num = num + ckdp( iband ) % ng
!!$    end do
!!$    if( iband .gt. nband ) then
!!$      write( 6, * ) 'Unexpected m'
!!$      write( 6, * ) m
!!$      stop
!!$    end if
!!$    ig = m - num
!!$    if( ig .gt. ckdp( iband ) % ng ) then
!!$      write( 6, * ) 'Unexpected ig'
!!$      write( 6, * ) iband, ig
!!$      stop
!!$    end if
!!$
!!$
!!$  end function m2ib

  !--------------------------------------------------------------------------------------

!!$  function m2ig( m ) result( ig )
!!$
!!$    use matype
!!$    use ckd_module
!!$
!!$    integer(i4b), intent(in ) :: m
!!$    integer(i4b)              :: ig
!!$
!!$
!!$    !
!!$    ! local variables
!!$    !
!!$    integer(i4b) :: iband
!!$    integer(i4b) :: num
!!$
!!$
!!$    num = 0
!!$    do iband = 1, nband
!!$      if( num + ckdp( iband ) % ng .ge. m ) exit
!!$      num = num + ckdp( iband ) % ng
!!$    end do
!!$    if( iband .gt. nband ) then
!!$      write( 6, * ) 'Unexpected m'
!!$      write( 6, * ) m
!!$      stop
!!$    end if
!!$    ig = m - num
!!$    if( ig .gt. ckdp( iband ) % ng ) then
!!$      write( 6, * ) 'Unexpected ig'
!!$      write( 6, * ) iband, ig
!!$      stop
!!$    end if
!!$
!!$
!!$  end function m2ig

  !--------------------------------------------------------------------------------------

!!$  subroutine getlnac_lblinterface( nwnsl, km, nras, gt, gp, ac, m )
!!$
!!$    use matype
!!$    use ckd_module
!!$
!!$    integer(i4b), intent(in ) :: nwnsl, km, nras
!!$    real(dp)    , intent(in ) :: gt( km ), gp( km )
!!$    real(dp)    , intent(out) :: ac( nwnsl, km, nras )
!!$    integer(i4b), intent(in ) :: m
!!$
!!$
!!$    !
!!$    ! local variables
!!$    !
!!$    real(dp)     :: ac_1d( km )
!!$    integer(i4b) :: ig, iband
!!$    integer(i4b) :: k
!!$
!!$
!!$    if( nwnsl .ne. 1 ) then
!!$      write( 6, * ) 'Unexpected nwnsl.'
!!$      write( 6, * ) nwnsl
!!$      stop
!!$    end if
!!$    if( nras .ne. 1 ) then
!!$      write( 6, * ) 'Unexpected nras.'
!!$      write( 6, * ) nras
!!$      stop
!!$    end if
!!$
!!$    call m2ckdpindices( m, ig, iband )
!!$
!!$    call getlnac_1d( km, gt, gp, ac_1d, ig, iband )
!!$
!!$    do k = 1, km
!!$      ac( :, k, 1 ) = ac_1d( k )
!!$    end do
!!$
!!$
!!$  end subroutine getlnac_lblinterface

  !--------------------------------------------------------------------------------------

!!$  subroutine getlnac_1d( km, gt, gp, ac, ig, iband )
!!$
!!$    use matype
!!$
!!$    integer(i4b), intent(in ) :: km
!!$    real(dp)    , intent(in ) :: gt( km ), gp( km )
!!$    real(dp)    , intent(out) :: ac( km )
!!$    integer(i4b), intent(in ) :: ig, iband
!!$
!!$
!!$    !
!!$    ! local variables
!!$    !
!!$    real(dp)     :: gt3d( 1, 1, km ), gp3d( 1, 1, km )
!!$    real(dp)     :: ac3d( 1, 1, km )
!!$    integer(i4b) :: k
!!$
!!$
!!$    do k = 1, km
!!$      gt3d( 1, 1, k ) = gt( k )
!!$      gp3d( 1, 1, k ) = gp( k )
!!$    end do
!!$
!!$    call getlnac( 1, 1, km, gt3d, gp3d, ac3d, ig, iband, 1, 1 )
!!$
!!$    do k = 1, km
!!$      ac( k ) = ac3d( 1, 1, k )
!!$    end do
!!$
!!$
!!$  end subroutine getlnac_1d

  !--------------------------------------------------------------------------------------
!!$
!!$  subroutine getlnac( im, jm, km, gt, gp, ac, ig, iband, ijs, ije )
!!$
!!$    use matype
!!$    use ckd_module
!!$
!!$    integer(i4b), intent(in ) :: im, jm, km
!!$    real(dp)    , intent(in ) :: gt( im, jm, km ), gp( im, jm, km )
!!$    real(dp)    , intent(out) :: ac( im, jm, km )
!!$    integer(i4b), intent(in ) :: ig, iband
!!$    integer(i4b), intent(in ) :: ijs, ije
!!$
!!$
!!$    !
!!$    ! local variables
!!$    !
!!$    real(dp)     :: glnp( im, jm, km )
!!$    real(dp)     :: lnac1, lnac2
!!$    integer(i4b) :: ij, k, l
!!$    integer(i4b) :: jj( im, jm, km ), kk( im, jm, km )
!!$
!!$
!!$    do k = 1, km
!!$      do ij = ijs, ije
!!$
!!$        glnp( ij, 1, k ) = log( gp( ij, 1, k ) + 1.0d-20 )
!!$
!!$      end do
!!$    end do
!!$
!!$    do k = 1, km
!!$      do ij = ijs, ije
!!$        kk( ij, 1, k ) = 1
!!$      end do
!!$    end do
!!$    do l = 1+1, ckdp( iband ) % nt - 1
!!$      do k = 1, km
!!$        do ij = ijs, ije
!!$          if( ckdp( iband ) % t( l ) .le. gt( ij, 1, k ) ) &
!!$            kk( ij, 1, k ) = l
!!$        end do
!!$      end do
!!$    end do
!!$
!!$    do k = 1, km
!!$      do ij = ijs, ije
!!$        jj( ij, 1, k ) = 1
!!$      end do
!!$    end do
!!$    do l = 1+1, ckdp( iband ) % nlnp - 1
!!$      do k = 1, km
!!$        do ij = ijs, ije
!!$          if( ckdp( iband ) % lnp( l ) .le. glnp( ij, 1, k ) ) &
!!$            jj( ij, 1, k ) = l
!!$        end do
!!$      end do
!!$    end do
!!$
!!$
!!$    do k = 1, km
!!$      do ij = ijs, ije
!!$
!!$        lnac1 &
!!$          = ( ckdp(iband)%lnac( ig, jj(ij,1,k)  , kk(ij,1,k)+1 )      &
!!$          - ckdp(iband)%lnac( ig, jj(ij,1,k)  , kk(ij,1,k)   ) )    &
!!$          / ( ckdp(iband)%t( kk(ij,1,k)+1 )                          &
!!$          - ckdp(iband)%t( kk(ij,1,k)   ) )                        &
!!$          * ( gt(ij,1,k) - ckdp( iband ) % t( kk(ij,1,k) ) )          &
!!$          + ckdp(iband)%lnac( ig, jj(ij,1,k)  , kk(ij,1,k)   )
!!$        lnac2 &
!!$          = ( ckdp(iband)%lnac( ig, jj(ij,1,k)+1, kk(ij,1,k)+1 )      &
!!$          - ckdp(iband)%lnac( ig, jj(ij,1,k)+1, kk(ij,1,k)   ) )    &
!!$          / ( ckdp(iband)%t( kk(ij,1,k)+1 )                          &
!!$          - ckdp(iband)%t( kk(ij,1,k)   ) )                        &
!!$          * ( gt(ij,1,k) - ckdp( iband ) % t( kk(ij,1,k) ) )          &
!!$          + ckdp(iband)%lnac( ig, jj(ij,1,k)+1, kk(ij,1,k)   )
!!$
!!$        ac( ij, 1, k ) &
!!$          = ( lnac2 - lnac1 ) &
!!$          / ( ckdp( iband ) % lnp( jj(ij,1,k)+1 ) &
!!$          - ckdp( iband ) % lnp( jj(ij,1,k)   ) ) &
!!$          * ( glnp(ij,1,k) - ckdp( iband ) % lnp( jj(ij,1,k) ) ) &
!!$          + lnac1
!!$      end do
!!$    end do
!!$
!!$
!!$  end subroutine getlnac

  !--------------------------------------------------------------------------------------

!!$  subroutine getpf_arr3d( im, jm, km, gt, gts, gp, gps, &
!!$    pfarr, pfsarr, ig, iband, ijs, ije )
!!$
!!$    use matype
!!$    use ckd_module
!!$    use pf_module
!!$
!!$    integer(i4b), intent(in ) :: im, jm, km
!!$    real(dp)    , intent(in ) :: gt( im, jm, km ), gts( im, jm ), gp( im, jm, km ), gps( im, jm )
!!$    real(dp)    , intent(out) :: pfarr( im, jm, km ), pfsarr( im, jm )
!!$    integer(i4b), intent(in ) :: ig, iband
!!$    integer(i4b), intent(in ) :: ijs, ije
!!$
!!$
!!$    !
!!$    ! local variables
!!$    !
!!$    real(dp)     :: pfr   ( im, jm, km )
!!$    real(dp)     :: gts3d1( im, jm, 1  ), gps3d1( im, jm, 1 ), pfr3d1( im, jm, 1 )
!!$    integer(i4b) :: ncp_pfint
!!$    integer(i4b) :: ij, k
!!$
!!$
!!$    ncp_pfint = 5
!!$
!!$    call getpfr( im, jm, km, gt    , gp    , pfr   , ig, iband, ijs, ije )
!!$    do ij = ijs, ije
!!$      gts3d1( ij, 1, 1 ) = gts( ij, 1 )
!!$      gps3d1( ij, 1, 1 ) = gps( ij, 1 )
!!$    end do
!!$    call getpfr( im, jm, 1 , gts3d1, gps3d1, pfr3d1, ig, iband, ijs, ije )
!!$
!!$
!!$    call pfint_gq_array( ckdp(iband)%wnbnds(1), ckdp(iband)%wnbnds(2), &
!!$      ncp_pfint, im, jm, km, gt, pfarr, &
!!$      ijs, ije )
!!$    call pfint_gq_array2d( ckdp(iband)%wnbnds(1), ckdp(iband)%wnbnds(2), &
!!$      ncp_pfint, im, jm, gts, pfsarr, &
!!$      ijs, ije )
!!$
!!$    do k = 1, km
!!$      do ij = ijs, ije
!!$        pfarr( ij, 1, k ) = pfarr( ij, 1, k ) &
!!$          / ( ckdp(iband)%wnbnds(2) - ckdp(iband)%wnbnds(1) )       &
!!$          * pfr( ij, 1, k )
!!$      end do
!!$    end do
!!$    do ij = ijs, ije
!!$      pfsarr( ij, 1 ) = pfsarr( ij, 1 ) &
!!$        / ( ckdp(iband)%wnbnds(2) - ckdp(iband)%wnbnds(1) )       &
!!$        * pfr3d1( ij, 1, 1 )
!!$    end do
!!$
!!$
!!$  end subroutine getpf_arr3d

  !--------------------------------------------------------------------------------------

!!$  subroutine getpf_arr3d_givenindices( im, jm, km, gt, gts, glnp, glnps, &
!!$    jj, kk, jjs, kks, &
!!$    pfarr, pfsarr, ig, iband, ijs, ije )
!!$
!!$    use matype
!!$    use ckd_module
!!$    use pf_module
!!$
!!$    integer(i4b), intent(in ) :: im, jm, km
!!$    real(dp)    , intent(in ) :: gt( im, jm, km ), gts( im, jm ), glnp( im, jm, km ), glnps( im, jm )
!!$    integer(i4b), intent(in ) :: jj ( im, jm, km ), kk ( im, jm, km )
!!$    integer(i4b), intent(in ) :: jjs( im, jm )    , kks( im, jm )
!!$    real(dp)    , intent(out) :: pfarr( im, jm, km ), pfsarr( im, jm )
!!$    integer(i4b), intent(in ) :: ig, iband
!!$    integer(i4b), intent(in ) :: ijs, ije
!!$
!!$
!!$    !
!!$    ! local variables
!!$    !
!!$    real(dp)     :: pfr   ( im, jm, km )
!!$    real(dp)     :: gts3d1( im, jm, 1  ), glnps3d1( im, jm, 1 ), pfr3d1( im, jm, 1 )
!!$    integer(i4b) :: jjs3d1( im, jm, 1  ), kks3d1  ( im, jm, 1 )
!!$    integer(i4b) :: ncp_pfint
!!$    integer(i4b) :: ij, k
!!$
!!$
!!$    ncp_pfint = 5
!!$
!!$
!!$    do ij = ijs, ije
!!$      gts3d1  ( ij, 1, 1 ) = gts  ( ij, 1 )
!!$      glnps3d1( ij, 1, 1 ) = glnps( ij, 1 )
!!$      jjs3d1  ( ij, 1, 1 ) = jjs  ( ij, 1 )
!!$      kks3d1  ( ij, 1, 1 ) = kks  ( ij, 1 )
!!$    end do
!!$
!!$
!!$    call getpfr_givenindices( im, jm, km, gt    , glnp    , jj    , kk    , pfr   , ig, iband, ijs, ije )
!!$    call getpfr_givenindices( im, jm, 1 , gts3d1, glnps3d1, jjs3d1, kks3d1, pfr3d1, ig, iband, ijs, ije )
!!$
!!$
!!$    call pfint_gq_array( ckdp(iband)%wnbnds(1), ckdp(iband)%wnbnds(2), &
!!$      ncp_pfint, im, jm, km, gt, pfarr, &
!!$      ijs, ije )
!!$    call pfint_gq_array2d( ckdp(iband)%wnbnds(1), ckdp(iband)%wnbnds(2), &
!!$      ncp_pfint, im, jm, gts, pfsarr, &
!!$      ijs, ije )
!!$
!!$    do k = 1, km
!!$      do ij = ijs, ije
!!$        pfarr( ij, 1, k ) = pfarr( ij, 1, k ) &
!!$          / ( ckdp(iband)%wnbnds(2) - ckdp(iband)%wnbnds(1) )       &
!!$          * pfr( ij, 1, k )
!!$      end do
!!$    end do
!!$    do ij = ijs, ije
!!$      pfsarr( ij, 1 ) = pfsarr( ij, 1 ) &
!!$        / ( ckdp(iband)%wnbnds(2) - ckdp(iband)%wnbnds(1) )       &
!!$        * pfr3d1( ij, 1, 1 )
!!$    end do
!!$
!!$
!!$  end subroutine getpf_arr3d_givenindices

  !--------------------------------------------------------------------------------------

!!$  subroutine getpf_arr3d_givenindices_norat( im, jm, km, gt, gts, glnp, glnps, &
!!$    pfarr, pfsarr, iband, ijs, ije )
!!$
!!$    use matype
!!$    use ckd_module
!!$    use pf_module
!!$
!!$    integer(i4b), intent(in ) :: im, jm, km
!!$    real(dp)    , intent(in ) :: gt( im, jm, km ), gts( im, jm ), glnp( im, jm, km ), glnps( im, jm )
!!$    real(dp)    , intent(out) :: pfarr( im, jm, km ), pfsarr( im, jm )
!!$    integer(i4b), intent(in ) :: iband
!!$    integer(i4b), intent(in ) :: ijs, ije
!!$
!!$
!!$    !
!!$    ! local variables
!!$    !
!!$    integer(i4b) :: ncp_pfint
!!$    integer(i4b) :: ij, k
!!$
!!$
!!$    ncp_pfint = 5
!!$
!!$
!!$    call pfint_gq_array( ckdp(iband)%wnbnds(1), ckdp(iband)%wnbnds(2), &
!!$      ncp_pfint, im, jm, km, gt, pfarr, &
!!$      ijs, ije )
!!$    call pfint_gq_array2d( ckdp(iband)%wnbnds(1), ckdp(iband)%wnbnds(2), &
!!$      ncp_pfint, im, jm, gts, pfsarr, &
!!$      ijs, ije )
!!$
!!$    do k = 1, km
!!$      do ij = ijs, ije
!!$        pfarr( ij, 1, k ) = pfarr( ij, 1, k ) &
!!$          / ( ckdp(iband)%wnbnds(2) - ckdp(iband)%wnbnds(1) )
!!$      end do
!!$    end do
!!$    do ij = ijs, ije
!!$      pfsarr( ij, 1 ) = pfsarr( ij, 1 ) &
!!$        / ( ckdp(iband)%wnbnds(2) - ckdp(iband)%wnbnds(1) )
!!$    end do
!!$
!!$
!!$  end subroutine getpf_arr3d_givenindices_norat

  !--------------------------------------------------------------------------------------

!!$  subroutine getpf_arr2d_givenindices( im, jm, gts, glnps, &
!!$    & jjs, kks, &
!!$    & pfsarr, ig, iband, ijs, ije )
!!$
!!$    use matype
!!$    use ckd_module
!!$    use pf_module
!!$
!!$    integer(i4b), intent(in ) :: im, jm
!!$    real(dp)    , intent(in ) :: gts( im, jm ), glnps( im, jm )
!!$    integer(i4b), intent(in ) :: jjs( im, jm )    , kks( im, jm )
!!$    real(dp)    , intent(out) :: pfsarr( im, jm )
!!$    integer(i4b), intent(in ) :: ig, iband
!!$    integer(i4b), intent(in ) :: ijs, ije
!!$
!!$
!!$    !
!!$    ! local variables
!!$    !
!!$    real(dp)     :: gts3d1( im, jm, 1  ), glnps3d1( im, jm, 1 ), pfr3d1( im, jm, 1 )
!!$    integer(i4b) :: jjs3d1( im, jm, 1  ), kks3d1  ( im, jm, 1 )
!!$    integer(i4b) :: ncp_pfint
!!$    integer(i4b) :: ij, k
!!$
!!$
!!$    ncp_pfint = 5
!!$
!!$
!!$    do ij = ijs, ije
!!$      gts3d1  ( ij, 1, 1 ) = gts  ( ij, 1 )
!!$      glnps3d1( ij, 1, 1 ) = glnps( ij, 1 )
!!$      jjs3d1  ( ij, 1, 1 ) = jjs  ( ij, 1 )
!!$      kks3d1  ( ij, 1, 1 ) = kks  ( ij, 1 )
!!$    end do
!!$
!!$
!!$    call getpfr_givenindices( im, jm, 1 , gts3d1, glnps3d1, jjs3d1, kks3d1, pfr3d1, ig, iband, ijs, ije )
!!$
!!$
!!$    call pfint_gq_array2d( ckdp(iband)%wnbnds(1), ckdp(iband)%wnbnds(2), &
!!$      ncp_pfint, im, jm, gts, pfsarr, &
!!$      ijs, ije )
!!$
!!$    do ij = ijs, ije
!!$      pfsarr( ij, 1 ) = pfsarr( ij, 1 ) &
!!$        / ( ckdp(iband)%wnbnds(2) - ckdp(iband)%wnbnds(1) )       &
!!$        * pfr3d1( ij, 1, 1 )
!!$    end do
!!$
!!$
!!$  end subroutine getpf_arr2d_givenindices

  !--------------------------------------------------------------------------------------

!!$  subroutine getpf_lblinterface( nwnsl, km, gt, gts, gp, gps, &
!!$    pfarr, pfsarr, iwnsls, iwnsle, m )
!!$
!!$    use matype
!!$    use ckd_module
!!$    use pf_module
!!$
!!$
!!$    integer(i4b), intent(in ) :: nwnsl, km
!!$    real(dp)    , intent(in ) :: gt( km ), gts, gp( km ), gps
!!$    real(dp)    , intent(out) :: pfarr( nwnsl, km ), pfsarr( nwnsl )
!!$    integer(i4b), intent(in ) :: iwnsls, iwnsle, m
!!$
!!$
!!$    !
!!$    ! local variables
!!$    !
!!$    real(dp)     :: pfr( km )
!!$    real(dp)     :: gt1( 1 ), gp1( 1 ), pfr1( 1 )
!!$    integer(i4b) :: ncp_pfint
!!$    integer(i4b) :: ig, iband
!!$    integer(i4b) :: k, iwnsl
!!$
!!$
!!$    call m2ckdpindices( m, ig, iband )
!!$
!!$    call getpfr_1d( km, gt , gp , pfr , ig, iband )
!!$    gt1( 1 ) = gts
!!$    gp1( 1 ) = gps
!!$    call getpfr_1d( 1 , gt1, gp1, pfr1, ig, iband )
!!$
!!$    ncp_pfint = 5
!!$
!!$    do k = 1, km
!!$      do iwnsl = iwnsls, iwnsle
!!$        pfarr( iwnsl, k ) &
!!$          = pfint_gq( ckdp(iband)%wnbnds(1), ckdp(iband)%wnbnds(2), &
!!$          ncp_pfint, gt( k ) )                          &
!!$          / ( ckdp(iband)%wnbnds(2) - ckdp(iband)%wnbnds(1) )       &
!!$          * pfr( k )
!!$      end do
!!$    end do
!!$    do iwnsl = iwnsls, iwnsle
!!$      pfsarr( iwnsl ) &
!!$        = pfint_gq( ckdp(iband)%wnbnds(1), ckdp(iband)%wnbnds(2), &
!!$        ncp_pfint, gts )                          &
!!$        / ( ckdp(iband)%wnbnds(2) - ckdp(iband)%wnbnds(1) )       &
!!$        * pfr1( 1 )
!!$    end do
!!$
!!$
!!$  end subroutine getpf_lblinterface

  !--------------------------------------------------------------------------------------

!!$  subroutine getpfr_1d( km, gt, gp, pfr, ig, iband )
!!$
!!$    use matype
!!$
!!$    integer(i4b), intent(in ) :: km
!!$    real(dp)    , intent(in ) :: gt( km ), gp( km )
!!$    real(dp)    , intent(out) :: pfr( km )
!!$    integer(i4b), intent(in ) :: ig, iband
!!$
!!$
!!$    !
!!$    ! local variables
!!$    !
!!$    real(dp)     :: gt3d ( 1, 1, km ), gp3d( 1, 1, km )
!!$    real(dp)     :: pfr3d( 1, 1, km )
!!$    integer(i4b) :: k
!!$
!!$
!!$    do k = 1, km
!!$      gt3d( 1, 1, k ) = gt( k )
!!$      gp3d( 1, 1, k ) = gp( k )
!!$    end do
!!$
!!$    call getpfr( 1, 1, km, gt3d, gp3d, pfr3d, ig, iband, 1, 1 )
!!$
!!$    do k = 1, km
!!$      pfr( k ) = pfr3d( 1, 1, k )
!!$    end do
!!$
!!$
!!$  end subroutine getpfr_1d

  !--------------------------------------------------------------------------------------

!!$  subroutine getpfr( im, jm, km, gt, gp, pfr, ig, iband, ijs, ije )
!!$
!!$    use matype
!!$    use ckd_module
!!$
!!$    integer(i4b), intent(in ) :: im, jm, km
!!$    real(dp)    , intent(in ) :: gt( im, jm, km ), gp( im, jm, km )
!!$    real(dp)    , intent(out) :: pfr( im, jm, km )
!!$    integer(i4b), intent(in ) :: ig, iband
!!$    integer(i4b), intent(in ) :: ijs, ije
!!$
!!$
!!$    !
!!$    ! local variables
!!$    !
!!$    real(dp)     :: glnp( im, jm, km )
!!$    real(dp)     :: pfr1, pfr2
!!$    integer(i4b) :: ij, k, l
!!$    integer(i4b) :: jj( im, jm, km ), kk( im, jm, km )
!!$
!!$
!!$    do k = 1, km
!!$      do ij = ijs, ije
!!$
!!$        glnp( ij, 1, k ) = log( gp( ij, 1, k ) + 1.0d-20 )
!!$
!!$      end do
!!$    end do
!!$
!!$
!!$    do k = 1, km
!!$      do ij = ijs, ije
!!$        kk( ij, 1, k ) = 1
!!$      end do
!!$    end do
!!$    do l = 1+1, ckdp( iband ) % nt - 1
!!$      do k = 1, km
!!$        do ij = ijs, ije
!!$          if( ckdp( iband ) % t( l ) .le. gt( ij, 1, k ) ) &
!!$            kk( ij, 1, k ) = l
!!$        end do
!!$      end do
!!$    end do
!!$
!!$    do k = 1, km
!!$      do ij = ijs, ije
!!$        jj( ij, 1, k ) = 1
!!$      end do
!!$    end do
!!$    do l = 1+1, ckdp( iband ) % nlnp - 1
!!$      do k = 1, km
!!$        do ij = ijs, ije
!!$          if( ckdp( iband ) % lnp( l ) .le. glnp( ij, 1, k ) ) &
!!$            jj( ij, 1, k ) = l
!!$        end do
!!$      end do
!!$    end do
!!$
!!$
!!$    do k = 1, km
!!$      do ij = ijs, ije
!!$
!!$        pfr1 &
!!$          = ( ckdp(iband)%pfr( ig, jj(ij,1,k)  , kk(ij,1,k)+1 )      &
!!$          - ckdp(iband)%pfr( ig, jj(ij,1,k)  , kk(ij,1,k)   ) )    &
!!$          / ( ckdp(iband)%t( kk(ij,1,k)+1 )                          &
!!$          - ckdp(iband)%t( kk(ij,1,k)   ) )                        &
!!$          * ( gt(ij,1,k) - ckdp( iband ) % t( kk(ij,1,k) ) )          &
!!$          + ckdp(iband)%pfr( ig, jj(ij,1,k)  , kk(ij,1,k)   )
!!$        pfr2 &
!!$          = ( ckdp(iband)%pfr( ig, jj(ij,1,k)+1, kk(ij,1,k)+1 )      &
!!$          - ckdp(iband)%pfr( ig, jj(ij,1,k)+1, kk(ij,1,k)   ) )    &
!!$          / ( ckdp(iband)%t( kk(ij,1,k)+1 )                          &
!!$          - ckdp(iband)%t( kk(ij,1,k)   ) )                        &
!!$          * ( gt(ij,1,k) - ckdp( iband ) % t( kk(ij,1,k) ) )          &
!!$          + ckdp(iband)%pfr( ig, jj(ij,1,k)+1, kk(ij,1,k)   )
!!$
!!$
!!$        pfr( ij, 1, k ) &
!!$          = ( pfr2 - pfr1 ) &
!!$          / ( ckdp( iband ) % lnp( jj(ij,1,k)+1 ) &
!!$          - ckdp( iband ) % lnp( jj(ij,1,k)   ) ) &
!!$          * ( glnp(ij,1,k) - ckdp( iband ) % lnp( jj(ij,1,k) ) ) &
!!$          + pfr1
!!$      end do
!!$    end do
!!$
!!$
!!$  end subroutine getpfr

  !--------------------------------------------------------------------------------------

!!$  subroutine getweight_lblinterface( nwnsl, weight, m )
!!$
!!$    use matype
!!$    use ckd_module
!!$    use pf_module
!!$
!!$    integer(i4b), intent(in ) :: nwnsl
!!$    real(dp)    , intent(out) :: weight
!!$    integer(i4b), intent(in ) :: m
!!$
!!$
!!$    !
!!$    ! local variables
!!$    !
!!$    integer(i4b) :: ig, iband
!!$
!!$
!!$    if( nwnsl .ne. 1 ) then
!!$      write( 6, * ) 'Now, nwnsl must be 1.'
!!$      write( 6, * ) nwnsl
!!$      stop
!!$    end if
!!$
!!$
!!$    call m2ckdpindices( m, ig, iband )
!!$
!!$    weight = ( ckdp(iband)%wnbnds(2) - ckdp(iband)%wnbnds(1) ) &
!!$      * ckdp(iband)%weight(ig)
!!$
!!$
!!$  end subroutine getweight_lblinterface

  !--------------------------------------------------------------------------------------

!!$  subroutine getweight_gcm( ig, iband, weight )
!!$
!!$    use matype
!!$    use ckd_module
!!$
!!$    integer(i4b), intent(in ) :: ig, iband
!!$    real(dp)    , intent(out) :: weight
!!$
!!$
!!$    !
!!$    ! local variables
!!$    !
!!$
!!$
!!$    weight = ( ckdp(iband)%wnbnds(2) - ckdp(iband)%wnbnds(1) ) &
!!$      * ckdp(iband)%weight(ig)
!!$
!!$
!!$  end subroutine getweight_gcm


  !**************************************************************************

  subroutine RadMars15mInit

    ! モジュール引用 ; USE statements
    !

    ! NAMELIST ファイル入力に関するユーティリティ
    ! Utilities for NAMELIST file input
    !
    use namelist_util, only: namelist_filename, NmlutilMsg, NmlutilAryValid

    ! ファイル入出力補助
    ! File I/O support
    !
    use dc_iounit, only: FileOpen

    use ckd_module, only : ckd_input, ckdp, nband


    !
    ! local variables
    !
    integer:: unit_nml        ! NAMELIST ファイルオープン用装置番号.
                              ! Unit number for NAMELIST file open
    integer:: iostat_nml      ! NAMELIST 読み込み時の IOSTAT.
                              ! IOSTAT of NAMELIST read

    character(STRING) :: rad15mkg_fn
    character(STRING) :: rad15mnf_fn

    integer           :: m

    namelist /rad_Mars_15m_nml/ &
      & rad15mkg_fn, &
!!$      & rad15mnf_fn, &
      & Rad15mInt


    ! 実行文 ; Executable statement
    !

    if ( rad_Mars_15m_inited ) return

    ! デフォルト値の設定
    ! Default values settings
    !

    rad15mkg_fn = "./kg15m"
!!$    rad15mnf_fn = "./nlte15mfactor"
    Rad15mInt   = 925.0_DP

    ! NAMELIST の読み込み
    ! NAMELIST is input
    !
    if ( trim(namelist_filename) /= '' ) then
      call FileOpen( unit_nml, &          ! (out)
        & namelist_filename, mode = 'r' ) ! (in)

      rewind( unit_nml )
      read( unit_nml,                &  ! (in)
        & nml = rad_Mars_15m_nml,    &  ! (out)
        & iostat = iostat_nml )         ! (out)
      close( unit_nml )

      call NmlutilMsg( iostat_nml, module_name ) ! (in)
    end if


!!$    allocate( rad_gp   ( im, jm,   km ) )
!!$    allocate( rad_gph  ( im, jm, 0:km ) )
!!$    allocate( rad_gt   ( im, jm,   km ) )
!!$    allocate( rad_gts  ( im, jm,   1  ) )
!!$    allocate( rad_gdod ( im, jm, 0:km ) )


    nras = 1
    nrps = 0

!!$    allocate( sgmh_f          ( km*nvr+1 ), &
!!$      &       sgm_f           ( km*nvr   ) )
!!$    allocate( gph_f    ( im, jm, km*nvr+1 ), &
!!$      &       gp_f     ( im, jm, km*nvr   ), &
!!$      &       gth_f    ( im, jm, km*nvr+1 ) )
!!$
!!$    allocate( gvmr_f   ( im, jm, km*nvr  , nras + nrps ) )
!!$    allocate( mmmass_f ( im, jm, km*nvr   ) )
!!$    allocate( ac_f     ( im, jm, km*nvr  , nras        ) )
!!$
!!$    allocate( gdod_f   ( im, jm, km*nvr+1 ) )
!!$
!!$
    allocate( xyra_Trans  (0:imax-1, 1:jmax, 0:kmax, 0:kmax) )
!!$    allocate( pfh_f    ( im, jm, km*nvr+1 ) )
!!$
!!$    allocate( uwflh_f  ( im, jm, km*nvr+1 ), &
!!$      &       dwflh_f  ( im, jm, km*nvr+1 ) )


    allocate( &
      & trans_i2i_toa(0:imax-1, 1:jmax, 0:kmax), &
      & trans_i2i_boa(0:imax-1, 1:jmax, 0:kmax), &
      & trans_i2i_s  (0:imax-1, 1:jmax, 0:kmax), &
      & trans_i2m_lli(0:imax-1, 1:jmax, 0:kmax, 1:kmax), &
      & trans_i2m_uli(0:imax-1, 1:jmax, 0:kmax, 1:kmax)  &
      & )

    trans_i2i_toa(:,:,:)   = 1.0d100
    trans_i2i_boa(:,:,:)   = 1.0d100
    trans_i2i_s  (:,:,:)   = 1.0d100
    trans_i2m_lli(:,:,:,:) = 1.0d100
    trans_i2m_uli(:,:,:,:) = 1.0d100


    !
    ! check
    !
    if( nras .ne. 1 ) then
      write( 6, * ) 'nras is not 1.'
      write( 6, * ) nras
      stop
    end if

    call ckd_input( rad15mkg_fn )

    ! check
    if( nband /= 1 ) then
      write( 6, * ) ' nband is not 1.'
      write( 6, * ) nband
      stop
    end if

    nwnl = 0
    do m = 1, nband
      nwnl = nwnl + ckdp( m ) % ng
    end do



!!$    call increase_vreso_boundary( km, nvr, sgmh, sgmh_f, "log" )
!!$    do k = 1, km * nvr
!!$      sgm_f( k ) = sqrt( sgmh_f( k ) * sgmh_f( k+1 ) )
!!$    end do


!!$    call rad15m_readnlte15mfac( rad15mnf_fn )


    !
    ! This routine must be called after rad15m_readkgtbl.
    !
!!$      call rad15m_rv_read( time )
!!$    call rad15m_rv_read_newscheme2006( time )


    rad_Mars_15m_inited = .true.


  end subroutine RadMars15mInit

  !**************************************************************************

!!$  subroutine rad15m_rv_put_newscheme2006( time, gt, gph, gp, gts, gdod, ijs, ije )
!!$
!!$    real(DP), intent(in ) :: time
!!$    real(DP), intent(in ) :: gt ( im, jm, km   )
!!$    real(DP), intent(in ) :: gph( im, jm, km+1 ), gp( im, jm, km )
!!$    real(DP), intent(in ) :: gts( im, jm )
!!$    real(DP), intent(in ) :: gdod( im, jm, km+1 )
!!$    integer , intent(in ) :: ijs, ije
!!$
!!$
!!$    !
!!$    ! local variables
!!$    !
!!$    integer  :: ij, k
!!$    real(DP) :: q15m( im, jm, km )
!!$    real(DP) :: gdf15m( im, jm )
!!$
!!$
!!$    q15m(:,:,:) = 1.0d100
!!$    gdf15m(:,:) = 1.0d100
!!$
!!$
!!$    rad_time = time
!!$
!!$
!!$    do k = 1, km
!!$      do ij = ijs, ije
!!$        rad_gt( ij, 1, k ) = gt( ij, 1, k )
!!$      end do
!!$    end do
!!$    do k = 1, km+1
!!$      do ij = ijs, ije
!!$        rad_gph( ij, 1, k-1 ) = gph( ij, 1, k )
!!$      end do
!!$    end do
!!$    do k = 1, km
!!$      do ij = ijs, ije
!!$        rad_gp( ij, 1, k ) = gp( ij, 1, k )
!!$      end do
!!$    end do
!!$    do ij = ijs, ije
!!$      rad_gts( ij, 1, 1 ) = gts( ij, 1 )
!!$    end do
!!$    do k = 1, km+1
!!$      do ij = ijs, ije
!!$        rad_gdod( ij, 1, k-1 ) = gdod( ij, 1, k )
!!$      end do
!!$    end do
!!$
!!$
!!$    sw_prep_rv = .true.
!!$
!!$
!!$  end subroutine rad15m_rv_put_newscheme2006

  !**************************************************************************

!!$  subroutine rad15m_rv_read_newscheme2006( timei )
!!$
!!$    use mamicro   , only : ntask, ijstart, ijend
!!$
!!$
!!$    real(DP), intent(in ) :: timei
!!$
!!$
!!$    !
!!$    ! local variables
!!$    !
!!$    integer, parameter :: mch = -1
!!$    logical     , parameter :: lgf = .false.
!!$
!!$    real(DP)                :: time
!!$    integer            :: ierr, ierr1, ierr2, ierr3, ierr4, ierr5
!!$    integer            :: iter
!!$
!!$
!!$    call ainini( rad_gp  , time, iter, ierr1, ihed(1), timei, &
!!$      &          mch     , lgf , im  , jm   , km   )
!!$    call ainini( rad_gph , time, iter, ierr2, ihed(2), timei, &
!!$      &          mch     , lgf , im  , jm   , km+1 )
!!$    call ainini( rad_gt  , time, iter, ierr3, ihed(3), timei, &
!!$      &          mch     , lgf , im  , jm   , km   )
!!$    call ainini( rad_gts , time, iter, ierr4, ihed(4), timei, &
!!$      &          mch     , lgf , im  , jm   , 1    )
!!$    call ainini( rad_gdod, time, iter, ierr5, ihed(5), timei, &
!!$      &          mch     , lgf , im  , jm   , km+1 )
!!$
!!$
!!$    ierr = max( ierr1, ierr2, ierr3, ierr4, ierr5 )
!!$
!!$    if( ierr .gt. 0 ) then
!!$      write( 6, * ) 'MARS: Radiative restart variables cannot be found.'
!!$    end if
!!$
!!$
!!$    rad_time = time
!!$
!!$
!!$    if( ierr .lt. 1 ) sw_prep_rv = .true.
!!$
!!$
!!$    if( ierr .lt. 1 ) then
!!$
!!$      write( 6, * ) 'MARS: Transmission is calculated.'
!!$
!!$      call rad15m_2006_rv_read_calctransW
!!$
!!$    end if
!!$
!!$
!!$  end subroutine rad15m_rv_read_newscheme2006

  !**************************************************************************

!!$  subroutine rad15m_2006_rv_read_calctransW
!!$
!!$    use mamicro   , only : ntask, ijstart, ijend
!!$
!!$
!!$    !
!!$    ! local variables
!!$    !
!!$    integer            :: m
!!$    integer            :: ij, k
!!$    integer            :: ijs, ije
!!$
!!$    real(DP)                :: &
!!$      gph( im, jm, km+1 ), gp( im, jm, km ), &
!!$      gt ( im, jm, km   ), gts( im, jm ), &
!!$      gdod( im, jm, km+1 )
!!$
!!$
!!$!cdir pardo for
!!$    do m = 1,ntask
!!$      ijs = ijstart( m )
!!$      ije = ijend( m )
!!$
!!$
!!$      do k = 1, km
!!$        do ij = ijs, ije
!!$          gp  ( ij, 1, k ) = rad_gp  ( ij, 1, k )
!!$          gt  ( ij, 1, k ) = rad_gt  ( ij, 1, k )
!!$        end do
!!$      end do
!!$      do ij = ijs, ije
!!$        gts( ij, 1 ) = rad_gts( ij, 1, 1 )
!!$      end do
!!$      do k = 0, km
!!$        do ij = ijs, ije
!!$          gph ( ij, 1, k+1 ) = rad_gph ( ij, 1, k )
!!$          gdod( ij, 1, k+1 ) = rad_gdod( ij, 1, k )
!!$        end do
!!$      end do
!!$
!!$
!!$      call rad15m_2006_rv_read_calctrans( &
!!$        & gt, gph, gp, gts, &
!!$        & gdod, &
!!$        & ijs, ije )
!!$
!!$    end do
!!$
!!$
!!$  end subroutine rad15m_2006_rv_read_calctransW

  !**************************************************************************

!!$  subroutine rad15m_rv_wrt_newscheme2006( ifl )
!!$
!!$    use maaxisr, only : signame, sighname
!!$
!!$
!!$    integer, intent(in ) :: ifl
!!$
!!$
!!$    !
!!$    ! local variables
!!$    !
!!$    integer     , parameter :: izero  = 0
!!$    integer(i8b)     , parameter :: notime = 0
!!$    character(len=lc), parameter :: nounit = '   '
!!$    character(len= 3), parameter :: dfmt   = 'UR8'
!!$    real(DP)         , parameter :: miss   = -1.0d30
!!$
!!$    integer                 :: iter, idate( 6 )
!!$
!!$
!!$    if( sw_prep_rv ) then
!!$      call ansymd( idate, rad_time )
!!$
!!$
!!$      call aonjrf( rad_gp  , ihed(1), titl(1), 'Pa'   , signame , &
!!$        &          km      , iter   , idate  , izero  , nounit  , &
!!$        &          notime  , idate  , notime , dfmt   , miss    , ifl )
!!$      call aonjrf( rad_gph , ihed(2), titl(2), 'Pa'   , sighname, &
!!$        &          km+1    , iter   , idate  , izero  , nounit  , &
!!$        &          notime  , idate  , notime , dfmt   , miss    , ifl )
!!$      call aonjrf( rad_gt  , ihed(3), titl(3), 'K'    , signame , &
!!$        &          km      , iter   , idate  , izero  , nounit  , &
!!$        &          notime  , idate  , notime , dfmt   , miss    , ifl )
!!$      call aonjrf( rad_gts , ihed(4), titl(4), 'K'    , sighname, &
!!$        &          1       , iter   , idate  , izero  , nounit  , &
!!$        &          notime  , idate  , notime , dfmt   , miss    , ifl )
!!$      call aonjrf( rad_gdod, ihed(5), titl(5), '1'    , sighname, &
!!$        &          km+1    , iter   , idate  , izero  , nounit  , &
!!$        &          notime  , idate  , notime , dfmt   , miss    , ifl )
!!$    end if
!!$
!!$
!!$
!!$  end subroutine rad15m_rv_wrt_newscheme2006

  !**************************************************************************




end module rad_Mars_15m