!2012/09/14
!2012/11/07(最終更新日)

!= 1次元放射対流平衡放射伝達モデル(Nakajima et al., 1992)
!
!= Radiative transfer with One-Dimensional Radiative-Convective Equilibrium Model
!
! Authors::   Masanori Onishi
! Version::   $Id: nakajima_RadiatConvEqModel.f90,v 1.0 2012-11-07 05:19:53 onishi Exp $ 
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!

module subprogs
  !
  != 1次元放射対流平衡放射伝達モデル(Nakajima et al., 1992)
  !
  != Radiative transfer with One-Dimensional Radiative-Convective Equilibrium Model
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  !
  !== Procedures List
  !
  ! func_MoistAd      :: 湿潤断熱減率の関数
  ! func_OptDep       :: 光学的厚さの関数
  ! func_Flux         :: Flux計算のための関数
  ! calc_AdiabLapse   :: 湿潤断熱減率で温度、フラックスを計算
  ! calc_RaiseTemp    :: 成層圏の温度を上昇
  ! calc_Tropopause   :: 対流圏界面の計算
  ! ------------            :: ------------
  ! func_MoistAd      :: function of moist adiabatic lapse rate
  ! func_OptDep       :: function of optical depth
  ! func_Flux         :: function for calculating flux
  ! calc_AdiabLapse   :: Calculate temperature and flux with moist adiabat
  ! calc_RaiseTemp    :: Raise temperature in stratosphere
  ! calc_Tropopause   :: Calculate to find toropopause
  !
  implicit none
  !
contains

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

  function func_MoistAd(                                       &
    & GasRUniv, LatentHeat, CpDry, CpWet, pv0, z_Press, z_Temp &
    & ) result(a)
    !
    ! 湿潤断熱減率の関数. 
    !
    ! function of moist adiabatic lapse rate
    !

    real(8), intent(in ) :: GasRUniv, LatentHeat, CpDry, CpWet, pv0, z_Press, z_Temp
    real(8) :: a, xv, xn

    ! 実行文 ; Executable statement
    !
    xv = ( pv0 / z_Press ) * exp( -LatentHeat / ( GasRUniv * z_Temp ) )
    xn = 1.0_8 - xv

    a = ( (( GasRUniv * z_Temp ) / ( z_Press * CpDry )) +        &
      &   ( xv / xn ) * ( LatentHeat / ( z_Press * CpDry ) ) ) / &
      & ( xn + ( xv * ( CpWet / CpDry ) ) +                      &
      &   ( xv / xn ) * ( LatentHeat**2 / ( GasRUniv * z_Temp*z_Temp * CpDry ) ) )

  end function func_MoistAd

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

  function func_OptDep(                                  & 
    & GasRUniv, LatentHeat, Grav, MolWtDry, MolWtWet,    &
    & AbsCoefDryCom, AbsCoefWetCom, pv0, z_Press, z_Temp &
    ) result(a)
    !
    ! 光学的厚さの関数.
    !
    ! function of optical depth
    !

    real(8), intent(in ) :: GasRUniv, LatentHeat, Grav, MolWtDry, MolWtWet, &
      & AbsCoefDryCom, AbsCoefWetCom, pv0, z_Press, z_Temp
    real(8) :: a, xv, xn

    ! 実行文 ; Executable statement
    !

    ! 湿潤大気
    ! moist atomosphere
    !
    xv = ( pv0 / z_Press ) * exp( -LatentHeat / ( GasRUniv * z_Temp ) )
    xn = 1.0_8 - xv

    a = ( AbsCoefWetCom * xv * MolWtWet + AbsCoefDryCom * xn * MolWtDry ) / &
      & ( ( xv * MolWtWet + xn * MolWtDry ) * Grav )

  end function func_OptDep

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

  function func_Flux( paiBkk, taukk, tauk ) result(a)
    !
    ! Flux計算のための関数.
    !
    ! function for calculating flux
    !

    real(8), intent(in ) :: paiBkk, taukk, tauk
    real(8) :: a

    ! 実行文 ; Executable statement
    !
    a = 1.5_8 * paiBkk * exp( -1.5_8 * ( taukk - tauk ) )

  end function func_Flux

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

  subroutine calc_AdiabLapse(                                        &
    & GasRUniv, StB, LatentHeat, Grav, MolWtDry, MolWtWet          , & ! (in )
    & AbsCoefDryCom, AbsCoefWetCom, p_v_0, CpWet, CpDry            , & ! (in )
    & eps, kmax, z_Press                                           , & ! (in )
    & knum, z_Temp, z_OptDep, z_RadLUwFlux, z_RadLDwFlux, z_FluxConv & ! (out)
    & )
    !
    ! 湿潤断熱減率で温度、フラックスを計算
    ! Calculate temperature and flux with moist adiabat
    !
    
    ! 宣言文 ; Declaration statements
    ! 
    real(8), intent(in ) :: GasRUniv
                              ! Universal gas constant
    real(8), intent(in ) :: StB
                              ! Stefan-Boltzmann constant
    real(8), intent(in ) :: Grav
                              ! Gravitational acceleration
    real(8), intent(in ) :: MolWtDry
                              ! Mean molecular weight of dry air
    real(8), intent(in ) :: CpDry
                              ! Specific heat of air at constant pressure
    real(8), intent(in ) :: MolWtWet
                              ! Mean molecular weight of condensible elements
    real(8), intent(in ) :: CpWet
                              ! Specific heat of condensible elements at constant pressure
    real(8), intent(in ) :: LatentHeat
                              ! Latent heat of condensation
    real(8), intent(in ) :: p_v_0
                              ! Constant for the water vapor saturation curve
    real(8), intent(in ) :: AbsCoefWetCom
                              ! Absorption coefficient of condensable component
    real(8), intent(in ) :: AbsCoefDryCom
                              ! Absorption coefiicient of noncondensable component
    real(8), intent(in ) :: eps
                              ! Resolution of calculation (temperature, optical depth, flux)
    integer, intent(in ) :: kmax
                              ! Grid settings: Number of vertical level
    real(8), intent(in ) :: z_Press      (1:kmax  )
                              ! pressure
    real(8), intent(out) :: z_Temp       (1:kmax  )
                              ! temperature
    real(8), intent(out) :: z_OptDep     (1:kmax  )
                              ! optical depth
    real(8), intent(out) :: z_FluxConv   (1:kmax-1)
                              ! flux convergence
    real(8), intent(out) :: z_RadLUwFlux (1:kmax  )
                              ! upoard longwave flux
    real(8), intent(out) :: z_RadLDwFlux (1:kmax  )
                              ! downward longwave flux
    integer, intent(out) :: knum         (1:kmax  )
                              ! work variable

    ! 作業変数
    ! Work variables
    !
    integer :: k, kk, kkk, kkkk, pDN
    real(8) :: k1, k2, k3, k4, ph, tauh, TempRK, preTempRK, DelOptDepDK, &
                 preDelOptDepDK, OptDepDK, preRadLFlux, err, errtau
    real(8) :: DelOptDep(1:kmax)

    ! 実行文 ; Executable statement
    !

    ! 温度、光学的厚さの鉛直構造を計算
    ! calculate vertical temperature and optical depth
    !
    ! temperature: Runge Kutta, optical depth: 台形法
    !
    do k = kmax, 2, -1

      pDN = int( z_Press(k) - z_Press(k-1) ) + 2 !初めに適当な刻み幅を与える
      ph = ( z_Press(k-1) - z_Press(k) ) / pDN
      TempRK = z_Temp(k)  !T(k)の値からT(k-1)を求める
      DelOptDepDK = 0.0_8 !Tau(k-1)-Tau(k)の間のDelOptDepDKをDelOptDep(k)とする

      DelOptDepDK = func_OptDep(                                      &
        & GasRUniv, LatentHeat, Grav, MolWtDry, MolWtWet,             &
        & AbsCoefDryCom, AbsCoefWetCom, p_v_0 , z_Press(k), z_Temp(k) &
        & ) * 0.5_8 * (-ph)

      do kk = 1, pDN

        k1 = func_MoistAd(                               &
          & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
          & z_Press(k) + ph * (kk-1)                   , &
          & TempRK                                       )
        k2 = func_MoistAd(                               &
          & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
          & z_Press(k) + ph * (kk-1) +  0.5_8 * ph     , &
          & TempRK                   + (0.5_8 * ph) * k1 )
        k3 = func_MoistAd(                               &
          & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
          & z_Press(k) + ph * (kk-1) +  0.5_8 * ph     , &
          & TempRK                   + (0.5_8 * ph) * k2 )
        k4 = func_MoistAd(                               &
          & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
          & z_Press(k) + ph * (kk-1) + ph              , &
          & TempRK                   + ph * k3           )

        TempRK = TempRK + &
          & ph * ( k1 + 2.0_8 * k2 + 2.0_8 * k3 + k4 ) / 6.0_8

        if (kk .ne. pDN) then
          DelOptDepDK = DelOptDepDK + func_OptDep(            &
            & GasRUniv, LatentHeat, Grav, MolWtDry, MolWtWet, &
            & AbsCoefDryCom, AbsCoefWetCom, p_v_0           , &
            & z_Press(k) + ph * (kk-1), TempRK                &
            & ) * (-ph)
        else
          DelOptDepDK = DelOptDepDK + func_OptDep(            &
            & GasRUniv, LatentHeat, Grav, MolWtDry, MolWtWet, &
            & AbsCoefDryCom, AbsCoefWetCom, p_v_0           , &
            & z_Press(k) + ph * (kk-1), TempRK                &
            & ) * 0.5_8 * (-ph)
        end if

      end do !kk

      preDelOptDepDK = DelOptDepDK
      preTempRK = TempRK

      do kkk = 1, 100 !精度が出るまでループ(上限は与える)

        TempRK = z_Temp(k) !積分値の初期化
        DelOptDepDK = 0.0_8
        pDN = pDN * 2
        ph = ( z_Press(k-1) - z_Press(k) ) / pDN

        DelOptDepDK = func_OptDep(                                     &
          & GasRUniv, LatentHeat, Grav, MolWtDry, MolWtWet           , &
          & AbsCoefDryCom, AbsCoefWetCom, p_v_0, z_Press(k), z_Temp(k) &
          & ) * 0.5_8 * (-ph)

        do kk = 1, pDN

          k1 = func_MoistAd(                               &
            & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
            & z_Press(k) + ph * (kk-1)                   , &
            & TempRK                                       )
          k2 = func_MoistAd(                               &
            & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
            & z_Press(k) + ph * (kk-1) +  0.5_8 * ph     , &
            & TempRK                   + (0.5_8 * ph) * k1 )
          k3 = func_MoistAd(                               &
            & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
            & z_Press(k) + ph * (kk-1) +  0.5_8 * ph     , &
            & TempRK                   + (0.5_8 * ph) * k2 )
          k4 = func_MoistAd(                               &
            & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
            & z_Press(k) + ph * (kk-1) + ph              , &
            & TempRK                   + ph * k3           )

          TempRK = TempRK + &
            & ph * ( k1 + 2.0_8 * k2 + 2.0_8 * k3 + k4 ) / 6.0_8

          if (kk .ne. pDN) then
            DelOptDepDK = DelOptDepDK + func_OptDep(            &
              & GasRUniv, LatentHeat, Grav, MolWtDry, MolWtWet, &
              & AbsCoefDryCom, AbsCoefWetCom, p_v_0           , &
              & z_Press(k) + ph * (kk-1), TempRK                &
              & ) * (-ph)
          else
            DelOptDepDK = DelOptDepDK + func_OptDep(            &
              & GasRUniv, LatentHeat, Grav, MolWtDry, MolWtWet, &
              & AbsCoefDryCom, AbsCoefWetCom, p_v_0           , &
              & z_Press(k) + ph * (kk-1), TempRK                &
              & ) * 0.5_8 * (-ph)
          end if

        end do !kk

        !!$ print *, "k=", k, "kkk=", kkk, "TempRK=", TempRK, "DelOptDepDK=", DelOptDepDK

        ! 収束判定
        ! convergence test
        !
        err = abs( ( TempRK - preTempRK ) / TempRK )
        errtau = abs( ( DelOptDepDK - preDelOptDepDK ) / DelOptDepDK )

        If (  ( err    < eps ) .or.                                             &
          &   (log(abs(TempRK-z_Temp(k))) < minexponent(abs(TempRK-z_Temp(k)))) &
          & .and.                                                               &
          &   ( errtau < eps ) .or.                                             &
          &   ( log(abs(DelOptDepDK)    ) < minexponent(abs(DelOptDepDK     ))) &
           ) then

          If ( log(abs(TempRK-z_Temp(k))) < minexponent(abs(TempRK-z_Temp(k)))) then
            TempRK = z_Temp(k)
            !!$ print *, "TempRK < 計算限界"
          end if
          If ( log(abs(DelOptDepDK     )) < minexponent(abs(DelOptDepDK     ))) then
            DelOptDepDK = 0.0_8
            !!$ print *, "DelOptDepDK < 計算限界"
          end if

          knum(k) = kkk
          !!$ print *, "err=", err, "errtau=", errtau, "kkk=", kkk
          exit
        end if

        preTempRK = TempRK
        preDelOptDepDK = DelOptDepDK
      end do !kkk

      z_Temp(k-1) = TempRK
      DelOptDep(k) = DelOptDepDK
    end do !k

    ! 全光学的厚さの計算
    ! calclate optical depth
    !
    z_OptDep(:) = 0.0_8
    do k = 1, kmax
      do kk = 1, k
        z_OptDep(k) = z_OptDep(k) + DelOptDep(kk)
      end do ! kk in k
    end do ! k

    ! Fluxの計算
    ! calculate flux
    !
    ! in: Temp, OptDep
    ! out: z_RadLUwFlux, z_RadLDwFlux
    !

    ! 上向きフラックスの計算
    ! calculate upward flux
    !
    z_RadLUwFlux(kmax) = StB * z_Temp(kmax)**4

    do k = kmax, 2, -1

      !==k~k-1層のFluxの計算
      pDN = ( abs(int(z_Press(k) - z_Press(k-1))) + 2) * 2 * (knum(k) + 1 ) !T, Tauの収束した刻み
      ph = ( z_Press(k-1) - z_Press(k) ) / pDN
      tauh = ( z_OptDep(k-1) - z_OptDep(k) ) / pDN
      TempRK = z_Temp(k)  !T(k)の値からT(k-1)に向かってsub gridの計算
      OptDepDK = z_OptDep(k) !Tau(k)の値からTau(k-1)に向かってsub gridの計算

      z_RadLUwFlux(k-1) = 0.5_8 * (-tauh) * (                             &
        & func_Flux(StB * z_Temp(k)**4  , z_OptDep(k)  , z_OptDep(k-1)) + &
        & func_Flux(StB * z_Temp(k-1)**4, z_OptDep(k-1), z_OptDep(k-1))   )

      OptDepDK = OptDepDK + func_OptDep(                             &
        & GasRUniv, LatentHeat, Grav, MolWtDry, MolWtWet,            &
        & AbsCoefDryCom, AbsCoefWetCom, p_v_0, z_Press(k), z_Temp(k) &
        & ) * 0.5_8 * ph

      do kk = 1, pDN-1 ! in k

        k1 = func_MoistAd(                               &
          & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
          & z_Press(k) + ph * (kk-1)                   , &
          & TempRK                                       )
        k2 = func_MoistAd(                               &
          & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
          & z_Press(k) + ph * (kk-1) +  0.5_8 * ph     , &
          & TempRK                   + (0.5_8 * ph) * k1 )
        k3 = func_MoistAd(                               &
          & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
          & z_Press(k) + ph * (kk-1) +  0.5_8 * ph     , &
          & TempRK                   + (0.5_8 * ph) * k2 )
        k4 = func_MoistAd(                               &
          & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
          & z_Press(k) + ph * (kk-1) + ph              , &
          & TempRK                   + ph * k3           )

        TempRK = TempRK + &
          & ph * ( k1 + 2.0_8 * k2 + 2.0_8 * k3 + k4 ) / 6.0_8
 
        if ( kk .ne. pDN ) then
          OptDepDK = OptDepDK + func_OptDep(                  &
            & GasRUniv, LatentHeat, Grav, MolWtDry, MolWtWet, &
            & AbsCoefDryCom, AbsCoefWetCom, p_v_0           , &
            & z_Press(k) + ph * (kk-1), TempRK                &
            & ) * ph
        else
          OptDepDK = OptDepDK + func_OptDep(                  &
            & GasRUniv, LatentHeat, Grav, MolWtDry, MolWtWet, &
            & AbsCoefDryCom, AbsCoefWetCom, p_v_0           , &
            & z_Press(k) + ph * (kk-1), TempRK                &
            & ) * 0.5_8 * ph
        end if
          z_RadLUwFlux(k-1) = z_RadLUwFlux(k-1) + (-tauh) * &
            & func_Flux( StB * TempRK**4  , OptDepDK  , z_OptDep(k-1) )
      end do !kk in k

      ! Fluxの収束判定
      ! convergence test for flux
      !
      preRadLFlux = z_RadLUwFlux(k-1)

      do kkk = 1, 100 !精度が出るまでループ(上限は与える
        pDN = pDN * 2 !2倍ずつ与える
        ph = ( z_Press(k-1) - z_Press(k) ) / pDN
        tauh = ( z_OptDep(k-1) - z_OptDep(k) ) / pDN
        TempRK = z_Temp(k)  !T(k)の値からT(k-1)に向かってsub gridの計算
        OptDepDK = z_OptDep(k) !Tau(k)の値からTau(k-1)に向かってsub gridの計算

        z_RadLUwFlux(k-1) = 0.5_8 * (-tauh) * (                             &
          & func_Flux(StB * z_Temp(k)**4  , z_OptDep(k)  , z_OptDep(k-1)) + &
          & func_Flux(StB * z_Temp(k-1)**4, z_OptDep(k-1), z_OptDep(k-1))   )

        OptDepDK = OptDepDK + func_OptDep(                             &
          & GasRUniv, LatentHeat, Grav, MolWtDry, MolWtWet,            &
          & AbsCoefDryCom, AbsCoefWetCom, p_v_0, z_Press(k), z_Temp(k) &
          & ) * 0.5_8 * ph

        do kk = 1, pDN-1 !in kkk in k

          k1 = func_MoistAd(                               &
            & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
            & z_Press(k) + ph * (kk-1)                   , &
            & TempRK                                       )
          k2 = func_MoistAd(                               &
            & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
            & z_Press(k) + ph * (kk-1) +  0.5_8 * ph     , &
            & TempRK                   + (0.5_8 * ph) * k1 )
          k3 = func_MoistAd(                               &
            & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
            & z_Press(k) + ph * (kk-1) +  0.5_8 * ph     , &
            & TempRK                   + (0.5_8 * ph) * k2 )
          k4 = func_MoistAd(                               &
            & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
            & z_Press(k) + ph * (kk-1) + ph              , &
            & TempRK             + ph * k3                 )

          TempRK = TempRK + &
            & ph * ( k1 + 2.0_8 * k2 + 2.0_8 * k3 + k4 ) / 6.0_8
  
          if ( kk .ne. pDN ) then
            OptDepDK = OptDepDK + func_OptDep(                  &
              & GasRUniv, LatentHeat, Grav, MolWtDry, MolWtWet, &
              & AbsCoefDryCom, AbsCoefWetCom, p_v_0           , &
              & z_Press(k) + ph * (kk-1), TempRK                &
              & ) * ph
          else
            OptDepDK = OptDepDK + func_OptDep(                  &
              & GasRUniv, LatentHeat, Grav, MolWtDry, MolWtWet, &
              & AbsCoefDryCom, AbsCoefWetCom, p_v_0           , &
              & z_Press(k) + ph * (kk-1), TempRK                &
              & ) * 0.5_8 * ph
          end if
            z_RadLUwFlux(k-1) = z_RadLUwFlux(k-1) + (-tauh) * &
              & func_Flux( StB * TempRK**4  , OptDepDK  , z_OptDep(k-1) )
        end do !kk in kkk, k

        ! 収束判定
        ! convergence test
        !
        err = abs( (z_RadLUwFlux(k-1) - preRadLFlux) / z_RadLUwFlux(k-1) )
        If ( ( err < eps ) .or.                                                    &
          &  ( log(abs(z_RadLUwFlux(k-1))) < minexponent(abs(z_RadLUwFlux(k-1))) ) &
          & ) then

          If (log(abs(z_RadLUwFlux(k-1))) < minexponent(abs(z_RadLUwFlux(k-1))) ) then
            z_RadLUwFlux(k-1) = 0.0_8
            !!$ print *, "FluxUp < 計算限界"
          end if
          !!$ print *, "err=", err, "kkk=", kkk
          exit
        end if
        preRadLFlux = z_RadLUwFlux(k-1)
      end do !kkk in k !====FluxUp収束完了====

      ! 1つ前の層までのFluxの寄与
      ! contribution of flux from pre layer
      !
      z_RadLUwFlux(k-1) = z_RadLUwFlux(k-1) + &
        & z_RadLUwFlux(k) * exp(-1.5_8 * (z_OptDep(k) - z_OptDep(k-1)))

    end do !k

    ! 下向きフラックスの計算
    ! calculate downward flux
    !
    z_RadLDwFlux(1) = 0.0_8

    do k = 2, kmax

      !==k~k+1層のFluxの計算
      pDN = ( abs(int(z_Press(k) - z_Press(k-1))) + 2) * 2 * (knum(k) + 1 ) !T, Tauの収束した刻み
      ph = ( z_Press(k-1) - z_Press(k) ) / pDN
      tauh = ( z_OptDep(k-1) - z_OptDep(k) ) / pDN
      TempRK = z_Temp(k)  !T(k)の値からT(k-1)に向かってsub gridの計算
      OptDepDK = z_OptDep(k) !Tau(k)の値からTau(k-1)に向かってsub gridの計算

      z_RadLDwFlux(k) = 0.5_8 * (-tauh) * (                           &
        & func_Flux(StB * z_Temp(k)**4  , z_OptDep(k), z_OptDep(k)) + &
        & func_Flux(StB * z_Temp(k-1)**4, z_OptDep(k), z_OptDep(k-1)) )

      OptDepDK = OptDepDK + func_OptDep(                             &
        & GasRUniv, LatentHeat, Grav, MolWtDry, MolWtWet           , &
        & AbsCoefDryCom, AbsCoefWetCom, p_v_0, z_Press(k), z_Temp(k) &
        & ) * 0.5_8 * ph

      do kk = 1, pDN-1 !in k

        k1 = func_MoistAd(                               &
          & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
          & z_Press(k) + ph * (kk-1)                   , &
          & TempRK                                       )
        k2 = func_MoistAd(                               &
          & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
          & z_Press(k) + ph * (kk-1) +  0.5_8 * ph     , &
          & TempRK                   + (0.5_8 * ph) * k1 )
        k3 = func_MoistAd(                               &
          & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
          & z_Press(k) + ph * (kk-1) +  0.5_8 * ph     , &
          & TempRK                   + (0.5_8 * ph) * k2 )
        k4 = func_MoistAd(                               &
          & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
          & z_Press(k) + ph * (kk-1) + ph              , &
          & TempRK                   + ph * k3           )

        TempRK = TempRK + &
          & ph * ( k1 + 2.0_8 * k2 + 2.0_8 * k3 + k4 ) / 6.0_8

        if ( kk .ne. pDN ) then
          OptDepDK = OptDepDK + func_OptDep(                  &
            & GasRUniv, LatentHeat, Grav, MolWtDry, MolWtWet, &
            & AbsCoefDryCom, AbsCoefWetCom, p_v_0           , &
            & z_Press(k) + ph * (kk-1), TempRK                &
            & ) * ph
        else
          OptDepDK = OptDepDK + func_OptDep(                  &
            & GasRUniv, LatentHeat, Grav, MolWtDry, MolWtWet, &
            & AbsCoefDryCom, AbsCoefWetCom, p_v_0           , &
            & z_Press(k) + ph * (kk-1), TempRK                &
            & ) * 0.5_8 * ph
        end if
          z_RadLDwFlux(k) = z_RadLDwFlux(k) + (-tauh) * &
            & func_Flux( StB * TempRK**4, z_OptDep(k), OptDepDK )

      end do !kk in k


      ! Fluxの収束判定
      ! convergence test for flux
      !
      preRadLFlux = z_RadLDwFlux(k)

      do kkk = 1, 100 ! in k !精度が出るまでループ(上限は与える)
        pDN = pDN * 2 !2倍ずつ与える
        ph = ( z_Press(k-1) - z_Press(k) ) / pDN
        tauh = ( z_OptDep(k-1) - z_OptDep(k) ) / pDN
        TempRK = z_Temp(k)  !T(k)の値からT(k-1)に向かってsub gridの計算
        OptDepDK = z_OptDep(k) !Tau(k)の値からTau(k-1)に向かってsub gridの計算

        z_RadLDwFlux(k) = 0.5_8 * (-tauh) * (                             &
          & func_Flux(StB * z_Temp(k)**4  , z_OptDep(k)  , z_OptDep(k)) + &
          & func_Flux(StB * z_Temp(k-1)**4, z_OptDep(k), z_OptDep(k-1))   )

        OptDepDK = OptDepDK + func_OptDep(                             &
          & GasRUniv, LatentHeat, Grav, MolWtDry, MolWtWet           , &
          & AbsCoefDryCom, AbsCoefWetCom, p_v_0, z_Press(k), z_Temp(k) &
          & ) * 0.5_8 * ph

        do kk = 1, pDN-1 !in kkk, k

          k1 = func_MoistAd(                               &
            & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
            & z_Press(k) + ph * (kk-1)                   , &
            & TempRK                                       )
          k2 = func_MoistAd(                               &
            & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
            & z_Press(k) + ph * (kk-1) +  0.5_8 * ph     , &
            & TempRK                   + (0.5_8 * ph) * k1 )
          k3 = func_MoistAd(                               &
            & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
            & z_Press(k) + ph * (kk-1) +  0.5_8 * ph     , &
            & TempRK                   + (0.5_8 * ph) * k2 )
          k4 = func_MoistAd(                               &
            & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
            & z_Press(k) + ph * (kk-1) + ph              , &
            & TempRK                   + ph * k3           )

          TempRK = TempRK + &
            & ph * ( k1 + 2.0_8 * k2 + 2.0_8 * k3 + k4 ) / 6.0_8

          if ( kk .ne. pDN ) then
            OptDepDK = OptDepDK + func_OptDep(                  & 
              & GasRUniv, LatentHeat, Grav, MolWtDry, MolWtWet, &
              & AbsCoefDryCom, AbsCoefWetCom, p_v_0           , &
              & z_Press(k) + ph * (kk-1), TempRK                &
              & ) * ph
          else
            OptDepDK = OptDepDK + func_OptDep(                  &
              & GasRUniv, LatentHeat, Grav, MolWtDry, MolWtWet, &
              & AbsCoefDryCom, AbsCoefWetCom, p_v_0           , &
              & z_Press(k) + ph * (kk-1), TempRK                &
              & ) * 0.5_8 * ph
          end if
          z_RadLDwFlux(k) = z_RadLDwFlux(k) + (-tauh) * &
            & func_Flux( StB * TempRK**4, z_OptDep(k), OptDepDK )

        end do !kk in kkk, k

        ! 収束判定
        ! convergence test
        !
        err = abs( (z_RadLDwFlux(k) - preRadLFlux) / z_RadLDwFlux(k) )
        If ( ( err < eps ) .or.                                                &
          &  ( log(abs(z_RadLDwFlux(k))) < minexponent(abs(z_RadLDwFlux(k))) ) &
          & ) then

          If (log(abs(z_RadLDwFlux(k))) < minexponent(abs(z_RadLDwFlux(k))) ) then
            z_RadLDwFlux(k) = 0.0_8
            !!$ print *, "FluxDown < 計算限界"
          end if
          !!$ print *, "err=", err, "kkk=", kkk
          exit
        end if
        preRadLFlux = z_RadLDwFlux(k)
      end do !kkk in k
 
      ! 1つ前の層までのFluxの寄与
      ! contribution of flux from pre layer
      !
      z_RadLDwFlux(k) = z_RadLDwFlux(k) + &
        & z_RadLDwFlux(k-1) * exp(-1.5_8 * (z_OptDep(k) - z_OptDep(k-1)))

    end do !k

    ! Flux Convergenceの計算
    ! calculation of flux convergence
    !
    ! z_FluxConv(k) is the flux convergence in the layer between grid(k) and grid(k+1)
    !
    do k = kmax-1, 1, -1
      z_FluxConv(k) = z_RadLUwFlux(k+1) - z_RadLUwFlux(k) - &
        &           ( z_RadLDwFlux(k+1) - z_RadLDwFlux(k) )
    end do

  end subroutine calc_AdiabLapse

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

  subroutine calc_RaiseTemp(                                   &
    & GasRUniv, StB, LatentHeat, Grav, MolWtDry, MolWtWet,     & ! (in ) 
    & AbsCoefDryCom, AbsCoefWetCom, p_v_0, CpWet, CpDry,       & ! (in )
    & eps, DelTime, kmax, ktp, knum, z_Press,                  & ! (in )
    & z_Temp, z_OptDep, z_RadLUwFlux, z_RadLDwFlux, z_FluxConv & ! (out)
    & )
    !
    ! 成層圏の温度を上昇
    ! Raise temperature in stratosphere
    !

    ! 宣言文 ; Declaration statements
    ! 
    real(8), intent(in ) :: GasRUniv
                              ! Universal gas constant
    real(8), intent(in ) :: StB
                              ! Stefan-Boltzmann constant
    real(8), intent(in ) :: Grav
                              ! Gravitational acceleration
    real(8), intent(in ) :: MolWtDry
                              ! Mean molecular weight of dry air
    real(8), intent(in ) :: CpDry
                              ! Specific heat of air at constant pressure
    real(8), intent(in ) :: MolWtWet
                              ! Mean molecular weight of condensible elements
    real(8), intent(in ) :: CpWet
                              ! Specific heat of condensible elements at constant pressure
    real(8), intent(in ) :: LatentHeat
                              ! Latent heat of condensation
    real(8), intent(in ) :: p_v_0
                              ! Constant for the water vapor saturation curve
    real(8), intent(in ) :: AbsCoefWetCom
                              ! Absorption coefficient of condensable component
    real(8), intent(in ) :: AbsCoefDryCom
                              ! Absorption coefiicient of noncondensable component
    real(8), intent(in ) :: eps
                              ! Resolution of calculation (temperature, optical depth, flux)
    real(8), intent(in ) :: DelTime
                              ! Time step
    integer, intent(in ) :: kmax
                              ! Grid settings: Number of vertical level
    integer, intent(in ) :: ktp
                              ! Grid number of toropopause
    real(8), intent(in ) :: z_Press      (1:kmax  )
                              ! pressure
    real(8), intent(out) :: z_Temp       (1:kmax  )
                              ! temperature
    real(8), intent(out) :: z_OptDep     (1:kmax  )
                              ! optical depth
    real(8), intent(out) :: z_FluxConv   (1:kmax-1)
                              ! flux convergence
    real(8), intent(out) :: z_RadLUwFlux (1:kmax  )
                              ! upoard longwave flux
    real(8), intent(out) :: z_RadLDwFlux (1:kmax  )
                              ! downward longwave flux
    integer, intent(out) :: knum         (1:kmax  )
                              ! work variable

    ! 作業変数
    ! Work variables
    !
    integer :: k, kk, kkk, kkkk, pDN
    real(8) :: k1, k2, k3, k4, ph, tauh, TempRK, OptDepDK, preRadLFlux, &
               err, x_v_max
    real(8) :: DelOptDep(1:kmax), DelTemp(1:kmax), x_n(1:kmax), x_v(1:kmax)

    ! 実行文 ; Executable statement
    !

    ! 加熱
    ! raising temperature by flux convergence
    !
    ! in: Temp, z_FluxConv, ktp, DelTime
    ! out: Temp, x

    ! 混合比 x の決定
    ! calculate mixing rate x
    !
    do k = 1, ktp
      x_v(k) = (p_v_0/z_Press(ktp)) * exp(-LatentHeat/(GasRUniv * z_Temp(ktp))) !まずは圏界面と同じ混合比
      x_v_max  = (p_v_0/z_Press(k)) * exp(-LatentHeat/(GasRUniv * z_Temp(k)  )) !飽和水蒸気がmax
      if(x_v(k) > x_v_max) then !飽和水蒸気量までしか含むことができない
         x_v(k) = x_v_max
      end if
      x_n(k) = 1.0_8 - x_v(k)
    end do ! k

    ! z_FluxConvで温度を上昇
    ! raise temperature by z_FluxConv
    !
    do k = 1, ktp-1
      DelTemp(k) = (                                              &
        & ( (0.5_8 * (x_n(k)  * MolWtDry + x_v(k)   * MolWtWet) + &
        &    0.5_8 * (x_n(k+1)* MolWtDry + x_v(k+1) * MolWtWet)   &
        &   ) * Grav                                              &
        & ) * z_FluxConv(k) * DelTime                             &
        &           ) / (                                         &
        & ( 0.5_8 * (x_n(k)   * CpDry + x_v(k)   * CpWet)       + &
        &   0.5_8 * (x_n(k+1) * CpDry + x_v(k+1) * CpWet)         &
        & ) * (z_Press(k+1) - z_Press(k))                         &
        &                )
    end do ! k

    z_Temp(1) = z_Temp(1) + DelTemp(1) !最も上の層だけ、下から2倍の加熱(上から加熱されない)
    do k = 2, ktp-1 !温度を上げるグリッドの上下のFluxConvergenceが正負なら温度変えない
      if(.not.((z_FluxConv(k-1) < 0.0_8) .and. (z_FluxConv(k) >  0.0_8))) then
        z_Temp(k) = z_Temp(k) + 0.5_8 * (DelTemp(k) + DelTemp(k-1))
      end if
    end do ! k

    ! 光学的厚さの決定
    ! calculate optical depth
    !
    ! in: Temp, ktp, x
    ! out: DelOptDep, OptDep
    !

    ! 混合比xの再決定(温度上昇後、光学的厚さを求めるため)
    !
    do k = 1, ktp
      x_v(k)  = (p_v_0/z_Press(ktp)) * exp(-LatentHeat/(GasRUniv * z_Temp(ktp)))
      x_v_max = (p_v_0/z_Press(k  )) * exp(-LatentHeat/(GasRUniv * z_Temp(k)  ))
      !!$ print *, "k=", k, "x_v=", x_v(k), "x_v_max=", x_v_max
      if(x_v(k) > x_v_max) then
        x_v(k) = x_v_max
      end if
      x_n(k) = 1.0_8 - x_v(k)
      !!$ print *, "x_v=", x_v(k)
    end do ! k

    ! 光学的厚さの計算
    ! calculate optical depth
    !
    do k = 2, ktp
      DelOptDep(k) = (                                            &
        & ( z_Press(k)-z_Press(k-1) ) / (                         &
        & ( 0.5_8 * (MolWtWet * x_v(k-1) + MolWtDry * x_n(k-1)) + &
        &   0.5_8 * (MolWtWet * x_v(k)   + MolWtDry * x_n(k)  ) ) &
        & * Grav )    ) * (                                       &
        & 0.5_8 * ( AbsCoefWetCom * MolWtWet * x_v(k-1)         + &
        &           AbsCoefDryCom * MolWtDry * x_n(k-1) )       + &
        & 0.5_8 * ( AbsCoefWetCom * MolWtWet * x_v(k)           + &
        &           AbsCoefDryCom * MolWtDry * x_n(k)   )         &
        &                 )
      !!$ print *, "DelOptDep(", k, ")=", DelOptDep(k)
    end do ! k

    do k = ktp+1, kmax
      DelOptDep(k) = z_OptDep(k) - z_OptDep(k-1)
      !!$ print *, "DelOptDep(", k, ")=", DelOptDep(k)
    end do ! k

    z_OptDep(1) = 0.0_8
    do k = 2, kmax
      z_OptDep(k) = z_OptDep(k-1) + DelOptDep(k)
      !!$ print *, Tau(k)
    end do ! k

    ! Fluxの計算
    ! calculate flux
    !
    ! in: Temp, OptDep, ktp
    ! out: z_RadLUwFlux, z_RadLDwFlux
    !
    ! T, tauを地表からRungeKuttaで再計算しながらFluxを求めていく
    !

    ! FluxUpの計算
    ! calculate upward flux
    !
    z_RadLUwFlux(kmax) = StB * z_Temp(kmax)**4

    do k = kmax, ktp+1, -1 !対流圏の計算

      !==k結-1層のFluxの計算
      pDN = ( abs(int(z_Press(k) - z_Press(k-1))) + 2) * 2 * (knum(k) + 1 ) !T, Tauの収束した刻み
      ph = ( z_Press(k-1) - z_Press(k) ) / pDN
      tauh = ( z_OptDep(k-1) - z_OptDep(k) ) / pDN
      TempRK = z_Temp(k)  !T(k)の値からT(k-1)に向かってsub gridの計算
      OptDepDK = z_OptDep(k) !Tau(k)の値からTau(k-1)に向かってsub gridの計算

      z_RadLUwFlux(k-1) = 0.5_8 * (-tauh) * (                             &
        & func_Flux(StB * z_Temp(k)**4  , z_OptDep(k)  , z_OptDep(k-1)) + &
        & func_Flux(StB * z_Temp(k-1)**4, z_OptDep(k-1), z_OptDep(k-1))   )

      OptDepDK = OptDepDK + func_OptDep(                             &
        & GasRUniv, LatentHeat, Grav, MolWtDry, MolWtWet           , &
        & AbsCoefDryCom, AbsCoefWetCom, p_v_0, z_Press(k), z_Temp(k) &
        & ) * 0.5_8 * ph

      do kk = 1, pDN-1 ! in k

        k1 = func_MoistAd(                               &
          & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
          & z_Press(k) + ph * (kk-1)                   , &
          & TempRK                                       )
        k2 = func_MoistAd(                               &
          & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
          & z_Press(k) + ph * (kk-1) +  0.5_8 * ph     , &
          & TempRK                   + (0.5_8 * ph) * k1 )
        k3 = func_MoistAd(                               &
          & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
          & z_Press(k) + ph * (kk-1) +  0.5_8 * ph     , &
          & TempRK                   + (0.5_8 * ph) * k2 )
        k4 = func_MoistAd(                               &
          & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
          & z_Press(k) + ph * (kk-1) + ph              , &
          & TempRK                   + ph * k3           )

        TempRK = TempRK + &
          & ph * ( k1 + 2.0_8 * k2 + 2.0_8 * k3 + k4 ) / 6.0_8

        if ( kk .ne. pDN ) then
          OptDepDK = OptDepDK + func_OptDep(                  &
            & GasRUniv, LatentHeat, Grav, MolWtDry, MolWtWet, &
            & AbsCoefDryCom, AbsCoefWetCom, p_v_0           , &
            & z_Press(k) + ph * (kk-1), TempRK ) * ph
        else
          OptDepDK = OptDepDK + func_OptDep(                  &
            & GasRUniv, LatentHeat, Grav, MolWtDry, MolWtWet, &
            & AbsCoefDryCom, AbsCoefWetCom, p_v_0           , &
            & z_Press(k) + ph * (kk-1), TempRK ) * 0.5_8 * ph
        end if
          z_RadLUwFlux(k-1) = z_RadLUwFlux(k-1) + (-tauh) * &
            & func_Flux( StB * TempRK**4  , OptDepDK  , z_OptDep(k-1) )
      end do !kk in k

      ! FluxUp収束するまで計算
      !
      preRadLFlux = z_RadLUwFlux(k-1)

      do kkk = 1, 100 ! in k !精度が出るまでループ(上限は与える
        pDN = pDN * 2 !2倍ずつ与える
        ph = ( z_Press(k-1) - z_Press(k) ) / pDN
        tauh = ( z_OptDep(k-1) - z_OptDep(k) ) / pDN
        TempRK = z_Temp(k)  !T(k)の値からT(k-1)に向かってsub gridの計算
        OptDepDK = z_OptDep(k) !Tau(k)の値からTau(k-1)に向かってsub gridの計算

        z_RadLUwFlux(k-1) = 0.5_8 * (-tauh) * (                             &
          & func_Flux(StB * z_Temp(k)**4  , z_OptDep(k)  , z_OptDep(k-1)) + &
          & func_Flux(StB * z_Temp(k-1)**4, z_OptDep(k-1), z_OptDep(k-1))   )

        OptDepDK = OptDepDK + func_OptDep(                             &
          & GasRUniv, LatentHeat, Grav, MolWtDry, MolWtWet           , &
          & AbsCoefDryCom, AbsCoefWetCom, p_v_0, z_Press(k), z_Temp(k) &
          & ) * 0.5_8 * ph

        do kk = 1, pDN-1 !in kkk, k

          k1 = func_MoistAd(                               &
            & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
            & z_Press(k) + ph * (kk-1)                   , &
            & TempRK                                       )
          k2 = func_MoistAd(                               &
            & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
            & z_Press(k) + ph * (kk-1) +  0.5_8 * ph     , &
            & TempRK                   + (0.5_8 * ph) * k1 )
          k3 = func_MoistAd(                               &
            & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
            & z_Press(k) + ph * (kk-1) +  0.5_8 * ph     , &
            & TempRK                   + (0.5_8 * ph) * k2 )
          k4 = func_MoistAd(                               &
            & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
            & z_Press(k) + ph * (kk-1) + ph              , &
            & TempRK             + ph * k3                 )

          TempRK = TempRK + &
            & ph * ( k1 + 2.0_8 * k2 + 2.0_8 * k3 + k4 ) / 6.0_8

          if ( kk .ne. pDN ) then
            OptDepDK = OptDepDK + func_OptDep(                  &
              & GasRUniv, LatentHeat, Grav, MolWtDry, MolWtWet, &
              & AbsCoefDryCom, AbsCoefWetCom, p_v_0           , &
              & z_Press(k) + ph * (kk-1), TempRK ) * ph
          else
            OptDepDK = OptDepDK + func_OptDep(                  &
              & GasRUniv, LatentHeat, Grav, MolWtDry, MolWtWet, &
              & AbsCoefDryCom, AbsCoefWetCom, p_v_0           , &
              & z_Press(k) + ph * (kk-1), TempRK ) * 0.5_8 * ph
          end if
            z_RadLUwFlux(k-1) = z_RadLUwFlux(k-1) + (-tauh) *   &
              & func_Flux( StB * TempRK**4, OptDepDK, z_OptDep(k-1) )
         end do !kk in kkk, k

         ! 収束判定
         ! convergence test
         !
         err = abs( (z_RadLUwFlux(k-1) - preRadLFlux) / z_RadLUwFlux(k-1) )
         If ( ( err < eps ) .or.                                                    &
           &  ( log(abs(z_RadLUwFlux(k-1))) < minexponent(abs(z_RadLUwFlux(k-1))) ) &
           & ) then

           If ( log(abs(z_RadLUwFlux(k-1))) < minexponent(abs(z_RadLUwFlux(k-1))) ) then
             z_RadLUwFlux(k-1) = 0.0_8
             !!$ print *, "FluxUp < 計算限界"
           end if
           !!$ print *, "err=", err, "kkk=", kkk
           exit
         end if
         preRadLFlux = z_RadLUwFlux(k-1)
      end do !kkk in k, kkkk !====FluxUp収束完了====

      ! 1つ前の層までのFluxの寄与
      ! contribution of flux from pre layer
      !
      z_RadLUwFlux(k-1) = z_RadLUwFlux(k-1) + &
        & z_RadLUwFlux(k) * exp(-1.5_8 * (z_OptDep(k) - z_OptDep(k-1)))

    end do !k

    do k = ktp, 2, -1 !k<=ktp成層圏!!!!!!FluxUp!!!!!

      pDN = ( abs(int(z_Press(k) - z_Press(k-1))) + 2 ) !!こちらは今回初めて
      ph = ( z_Press(k-1) - z_Press(k) ) / pDN
      tauh = ( z_OptDep(k-1) - z_OptDep(k) ) / pDN

      !台形法
      !台形量端の計算!!FluxUp
      z_RadLUwFlux(k-1) = 0.5_8 * (-tauh) * (                         &
        & func_Flux(StB * z_Temp(k)**4  , z_OptDep(k)  , z_OptDep(k-1)) + &
        & func_Flux(StB * z_Temp(k-1)**4, z_OptDep(k-1), z_OptDep(k-1))   )

      !台形間の計算!!FluxUp
      do kk = 1, pDN-1 !in k
         TempRK  = z_Temp(k)                       + &
           & (z_Temp(k-1) - z_Temp(k)) * (ph * kk) / &
           & (z_Press(k-1) - z_Press(k))             ! 直線

         OptDepDK  = z_OptDep(k)                       + &
           & (z_OptDep(k-1) - z_OptDep(k)) * (ph * kk) / &
           & (z_Press(k-1) - z_Press(k))                 ! 直線

         z_RadLUwFlux(k-1) = z_RadLUwFlux(k-1) + &
           & (-tauh) * func_Flux( StB * TempRK**4, OptDepDK, z_OptDep(k-1) )

      end do !kk in k !!FluxUp

      ! FluxUp収束するまで計算
      !
      preRadLFlux = z_RadLUwFlux(k-1)

      do kkk = 1, 100 !in k !精度が出るまでループ(上限は与える
        pDN = pDN * 2 !2倍ずつ与える
        ph = ( z_Press(k-1) - z_Press(k) ) / pDN
        tauh = ( z_OptDep(k-1) - z_OptDep(k) ) / pDN

        !台形法
        !台形量端の計算!!FluxUp
        z_RadLUwFlux(k-1) = 0.5_8 * (-tauh) * (                             &
          & func_Flux(StB * z_Temp(k)**4  , z_OptDep(k)  , z_OptDep(k-1)) + &
          & func_Flux(StB * z_Temp(k-1)**4, z_OptDep(k-1), z_OptDep(k-1))   )

        !台形間の計算!!FluxUp
        do kk = 1, pDN-1 !in kkk, k
          TempRK  = z_Temp(k)                       + &
            & (z_Temp(k-1) - z_Temp(k)) * (ph * kk) / &
            & (z_Press(k-1) - z_Press(k))             ! 直線

          OptDepDK  = z_OptDep(k)                       + &
            & (z_OptDep(k-1) - z_OptDep(k)) * (ph * kk) / &
            & (z_Press(k-1) - z_Press(k))                 ! 直線

          z_RadLUwFlux(k-1) = z_RadLUwFlux(k-1) + &
            & (-tauh) * func_Flux( StB * TempRK**4, OptDepDK, z_OptDep(k-1) )

        end do !kk in kkk, k

        ! 収束判定
        ! convergence test
        !
        err = abs((z_RadLUwFlux(k-1) - preRadLFlux)/z_RadLUwFlux(k-1))
        If ( ( err < eps ) .or.                                                    &
          &  ( log(abs(z_RadLUwFlux(k-1))) < minexponent(abs(z_RadLUwFlux(k-1))) ) &
          & ) then

          If (log(abs(z_RadLUwFlux(k-1))) < minexponent(abs(z_RadLUwFlux(k-1))) ) then
            z_RadLUwFlux(k-1) = 0.0_8
            !!$ print *, "FluxUp < 計算限界"
          end if
          !!$ print *, "err=", err, "kkk=", kkk
          exit
        end if
        preRadLFlux = z_RadLUwFlux(k-1)
      end do !kkk in k !====FluxUp収束完了====

      ! 1つ前の層までのFluxの寄与
      ! contribution of flux from pre layer
      !
      z_RadLUwFlux(k-1) = z_RadLUwFlux(k-1) + &
        & z_RadLUwFlux(k) * exp(-1.5_8 * (z_OptDep(k) - z_OptDep(k-1)))

    end do !k !!FluxUp計算終了!!===========FluxUp

    ! ここからFluxDownの計算
    ! calculate downward flux hereafter
    !
    z_RadLDwFlux(1) = 0.0_8

    do k = 2, ktp !k<=ktp成層圏!!FluxDown

      pDN = ( abs(int(z_Press(k) - z_Press(k-1))) + 2 ) !!こちらは今回初めて
      ph = ( z_Press(k-1) - z_Press(k) ) / pDN
      tauh = ( z_OptDep(k-1) - z_OptDep(k) ) / pDN

      !台形法
      !台形量端の計算 !!FluxDown
      z_RadLDwFlux(k) = 0.5_8 * (-tauh) * ( &
        & func_Flux(StB * z_Temp(k)**4  , z_OptDep(k)  , z_OptDep(k)) + &
        & func_Flux(StB * z_Temp(k-1)**4, z_OptDep(k), z_OptDep(k-1))   )

      !台形間の計算 !!FluxDown
      do kk = 1, pDN-1 ! in k
        TempRK  = z_Temp(k)                       + &
          & (z_Temp(k-1) - z_Temp(k)) * (ph * kk) / &
          & (z_Press(k-1) - z_Press(k))             ! 直線

        OptDepDK  = z_OptDep(k)                       + &
          & (z_OptDep(k-1) - z_OptDep(k)) * (ph * kk) / &
          & (z_Press(k-1) - z_Press(k))                 ! 直線

        z_RadLDwFlux(k) = z_RadLDwFlux(k) + &
          & (-tauh) * func_Flux( StB * TempRK**4, z_OptDep(k), OptDepDK )
      end do !kk in k !!FluxDown

      ! FluxDown収束するまで計算
      preRadLFlux = z_RadLDwFlux(k)

      do kkk = 1, 100 !in k !精度が出るまでループ(上限は与える
        pDN = pDN * 2 !2倍ずつ与える
        ph = ( z_Press(k-1) - z_Press(k) ) / pDN
        tauh = ( z_OptDep(k-1) - z_OptDep(k) ) / pDN

        !台形法
        !台形量端の計算 !!FluxDown
        z_RadLDwFlux(k) = 0.5_8 * (-tauh) * (                             &
          & func_Flux(StB * z_Temp(k)**4  , z_OptDep(k)  , z_OptDep(k)) + &
          & func_Flux(StB * z_Temp(k-1)**4, z_OptDep(k), z_OptDep(k-1))   )

        !台形間の計算 !!FluxDown
        do kk = 1, pDN-1 ! in kkk, k
          TempRK  = z_Temp(k)                       + &
            & (z_Temp(k-1) - z_Temp(k)) * (ph * kk) / &
            & (z_Press(k-1) - z_Press(k))             ! 直線

          OptDepDK  = z_OptDep(k)                       + &
            & (z_OptDep(k-1) - z_OptDep(k)) * (ph * kk) / &
            & (z_Press(k-1) - z_Press(k))                 ! 直線

          z_RadLDwFlux(k) = z_RadLDwFlux(k) + &
            & (-tauh) * func_Flux( StB * TempRK**4, z_OptDep(k), OptDepDK )
        end do !kk in kkk, k !!FluxDown

        ! 収束判定
        ! convergence test
        !
        err = abs( (z_RadLDwFlux(k) - preRadLFlux) / z_RadLDwFlux(k) )
        If ( ( err < eps ) .or. &
          &  ( log(abs(z_RadLDwFlux(k))) < minexponent(abs(z_RadLDwFlux(k))) ) &
          & ) then

          If (log(abs(z_RadLDwFlux(k))) < minexponent(abs(z_RadLDwFlux(k))) ) then
            z_RadLDwFlux(k) = 0.0_8
            !!$ print *, "FluxDown < 計算限界"
          end if
          !!$ print *, "err=", err, "kkk=", kkk
          exit
        end if
        preRadLFlux = z_RadLDwFlux(k)
      end do !kkk in k

      ! 1つ前の層までのFluxの寄与
      ! contribution of flux from pre layer
      !
      z_RadLDwFlux(k) = z_RadLDwFlux(k) + &
        & z_RadLDwFlux(k-1) * exp(-1.5_8 * (z_OptDep(k) - z_OptDep(k-1)))

    end do !k ! FluxDown成層圏計算終了

    ! FluxDownの計算 対流圏
    !
    do k = ktp+1, kmax

      !==k結+1層のFluxの計算
      pDN = ( abs(int(z_Press(k) - z_Press(k-1))) + 2) * 2 * (knum(k) + 1 ) !T, Tauの収束した刻み
      ph = ( z_Press(k-1) - z_Press(k) ) / pDN
      tauh = ( z_OptDep(k-1) - z_OptDep(k) ) / pDN
      TempRK = z_Temp(k)  !T(k)の値からT(k-1)に向かってsub gridの計算
      OptDepDK = z_OptDep(k) !Tau(k)の値からTau(k-1)に向かってsub gridの計算

      z_RadLDwFlux(k) = 0.5_8 * (-tauh) * (                             &
        & func_Flux(StB * z_Temp(k)**4  , z_OptDep(k)  , z_OptDep(k)) + &
        & func_Flux(StB * z_Temp(k-1)**4, z_OptDep(k), z_OptDep(k-1))   )

      OptDepDK= OptDepDK + &
             func_OptDep(GasRUniv, LatentHeat, Grav                 , &
             MolWtDry, MolWtWet, AbsCoefDryCom, AbsCoefWetCom, p_v_0, &
             z_Press(k), z_Temp(k)  ) * 0.5_8 * ph

      do kk = 1, pDN-1 !in k

        k1 = func_MoistAd(                               &
          & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
          & z_Press(k) + ph * (kk-1)                   , &
          & TempRK                                       )
        k2 = func_MoistAd(                               &
          & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
          & z_Press(k) + ph * (kk-1) +  0.5_8 * ph     , &
          & TempRK                   + (0.5_8 * ph) * k1 )
        k3 = func_MoistAd(                               &
          & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
          & z_Press(k) + ph * (kk-1) +  0.5_8 * ph     , &
          & TempRK                   + (0.5_8 * ph) * k2 )
        k4 = func_MoistAd(                               &
          & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
          & z_Press(k) + ph * (kk-1) + ph              , &
          & TempRK                   + ph * k3           )

        TempRK = TempRK + &
          & ph * ( k1 + 2.0_8 * k2 + 2.0_8 * k3 + k4 ) / 6.0_8

        if ( kk .ne. pDN ) then
          OptDepDK = OptDepDK + func_OptDep(                  &
            & GasRUniv, LatentHeat, Grav, MolWtDry, MolWtWet, &
            & AbsCoefDryCom, AbsCoefWetCom, p_v_0           , &
            & z_Press(k) + ph * (kk-1), TempRK ) * ph
        else
          OptDepDK = OptDepDK + func_OptDep(                  &
            & GasRUniv, LatentHeat, Grav, MolWtDry, MolWtWet, &
            & AbsCoefDryCom, AbsCoefWetCom, p_v_0           , &
            & z_Press(k) + ph * (kk-1), TempRK ) * 0.5_8 * ph
        end if
          z_RadLDwFlux(k) = z_RadLDwFlux(k) + (-tauh) * &
            & func_Flux( StB * TempRK**4, z_OptDep(k), OptDepDK )

      end do !kk in k

      ! FluxDown収束するまで計算
      !
      preRadLFlux = z_RadLDwFlux(k)

      do kkk = 1, 100 ! in k !精度が出るまでループ(上限は与える)
        pDN = pDN * 2 !2倍ずつ与える
        ph = ( z_Press(k-1) - z_Press(k) ) / pDN
        tauh = ( z_OptDep(k-1) - z_OptDep(k) ) / pDN
        TempRK = z_Temp(k)  !T(k)の値からT(k-1)に向かってsub gridの計算
        OptDepDK = z_OptDep(k) !Tau(k)の値からTau(k-1)に向かってsub gridの計算

        z_RadLDwFlux(k) = 0.5_8 * (-tauh) * (                             &
          & func_Flux(StB * z_Temp(k)**4  , z_OptDep(k)  , z_OptDep(k)) + &
          & func_Flux(StB * z_Temp(k-1)**4, z_OptDep(k), z_OptDep(k-1))   )

        OptDepDK = OptDepDK + func_OptDep(                             &
          & GasRUniv, LatentHeat, Grav, MolWtDry, MolWtWet           , &
          & AbsCoefDryCom, AbsCoefWetCom, p_v_0, z_Press(k), z_Temp(k) &
          & ) * 0.5_8 * ph

        do kk = 1, pDN-1 !in kkk, k

          k1 = func_MoistAd(                               &
            & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
            & z_Press(k) + ph * (kk-1)                   , &
            & TempRK                                       )
          k2 = func_MoistAd(                               &
            & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
            & z_Press(k) + ph * (kk-1) +  0.5_8 * ph     , &
            & TempRK                   + (0.5_8 * ph) * k1 )
          k3 = func_MoistAd(                               &
            & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
            & z_Press(k) + ph * (kk-1) +  0.5_8 * ph     , &
            & TempRK                   + (0.5_8 * ph) * k2 )
          k4 = func_MoistAd(                               &
            & GasRUniv, LatentHeat, CpDry, CpWet, p_v_0  , &
            & z_Press(k) + ph * (kk-1) + ph              , &
            & TempRK                   + ph * k3           )

          TempRK = TempRK + &
            & ph * ( k1 + 2.0_8 * k2 + 2.0_8 * k3 + k4 ) / 6.0_8

          if ( kk .ne. pDN ) then
            OptDepDK = OptDepDK + func_OptDep(                  &
              & GasRUniv, LatentHeat, Grav, MolWtDry, MolWtWet, &
              & AbsCoefDryCom, AbsCoefWetCom, p_v_0           , &
              & z_Press(k) + ph * (kk-1), TempRK ) * ph
          else
            OptDepDK = OptDepDK + func_OptDep(                  &
              & GasRUniv, LatentHeat, Grav, MolWtDry, MolWtWet, &
              & AbsCoefDryCom, AbsCoefWetCom, p_v_0           , &
              & z_Press(k) + ph * (kk-1), TempRK ) * 0.5_8 * ph
          end if
            z_RadLDwFlux(k) = z_RadLDwFlux(k) + (-tauh) * &
              & func_Flux( StB * TempRK**4, z_OptDep(k), OptDepDK )

        end do !kk in kkk, k

        ! 収束判定
        ! convergence test
        !
        err = abs( (z_RadLDwFlux(k) - preRadLFlux) / z_RadLDwFlux(k) )
        If ( ( err < eps ) .or.                                                &
          &  ( log(abs(z_RadLDwFlux(k))) < minexponent(abs(z_RadLDwFlux(k))) ) &
          & ) then

          If ( log(abs(z_RadLDwFlux(k))) < minexponent(abs(z_RadLDwFlux(k))) ) then
            z_RadLDwFlux(k) = 0.0_8
            !!$ print *, "FluxDown < 計算限界"
          end if
          !!$ print *, "err=", err, "kkk=", kkk
          exit
        end if
        preRadLFlux = z_RadLDwFlux(k)
      end do !kkk in k

      ! 1つ前の層までのFluxの寄与
      ! contribution of flux from pre layer
      !
      z_RadLDwFlux(k) = z_RadLDwFlux(k) + &
        & z_RadLDwFlux(k-1) * exp(-1.5_8 * (z_OptDep(k) - z_OptDep(k-1)))

    end do !k !==FluxDown対流圏計算終了

    ! FluxConvergenceの決定
    ! calculate flux convergence
    !
    ! in: FluxUp, FluxDown
    ! out:z_FluxConv
    !
    ! z_FluxConv(k) is the flux convergence in the layer between grid(k) and grid(k+1)
    !

    do k = kmax-1, 1, -1
      z_FluxConv(k) = z_RadLUwFlux(k+1) - z_RadLUwFlux(k) - &
        &           ( z_RadLDwFlux(k+1) - z_RadLDwFlux(k) )
    end do ! k

  end subroutine calc_RaiseTemp

  !----------------------------------------------------------------------------
  
  subroutine calc_Tropopause(                                            &
    & GasRUniv, StB, LatentHeat, Grav, MolWtDry, MolWtWet,               & ! (in )
    & AbsCoefDryCom, AbsCoefWetCom, p_v_0, eps, kmax,                    & ! (in )
    & z_Press, z_Temp, z_OptDep, z_RadLUwFlux, z_RadLDwFlux, z_FluxConv, & ! (in )
    & val_tropo                                                          & ! (out)
    & )
    !
    ! 対流圏界面の計算
    ! Calculate to find toropopause
    !

    ! 宣言文 ; Declaration statements
    ! 
    real(8), intent(in ) :: GasRUniv
                              ! Universal gas constant
    real(8), intent(in ) :: StB
                              ! Stefan-Boltzmann constant
    real(8), intent(in ) :: Grav
                              ! Gravitational acceleration
    real(8), intent(in ) :: MolWtDry
                              ! Mean molecular weight of dry air
    real(8), intent(in ) :: MolWtWet
                              ! Mean molecular weight of condensible elements
    real(8), intent(in ) :: LatentHeat
                              ! Latent heat of condensation
    real(8), intent(in ) :: p_v_0
                              ! Constant for the water vapor saturation curve
    real(8), intent(in ) :: AbsCoefWetCom
                              ! Absorption coefficient of condensable component
    real(8), intent(in ) :: AbsCoefDryCom
                              ! Absorption coefiicient of noncondensable component
    real(8), intent(in ) :: eps
                              ! Resolution of calculation (temperature, optical depth, flux)
    integer, intent(in ) :: kmax
                              ! Grid settings: Number of vertical level
    real(8), intent(in ) :: z_Press      (1:kmax  )
                              ! pressure
    real(8), intent(in ) :: z_Temp       (1:kmax  )
                              ! temperature
    real(8), intent(in ) :: z_OptDep     (1:kmax  )
                              ! optical depth
    real(8), intent(in ) :: z_FluxConv   (1:kmax-1)
                              ! flux convergence
    real(8), intent(in ) :: z_RadLUwFlux (1:kmax  )
                              ! upoard longwave flux
    real(8), intent(in ) :: z_RadLDwFlux (1:kmax  )
                              ! downward longwave flux
    real(8), intent(out) :: val_tropo    (1:7     )
                              ! values of toropopause

    ! 作業変数
    ! Work variables
    !
    integer :: k, kk, kkk, kkkk, pDN, ktp, ktpRN, numktp
    real(8) :: ph, tauh, TempRK, OptDepDK, preRadLFlux,               &
      & err, x_v_max, FluxConv_plus, FluxConv_minus
    real(8) :: subPress(1:5), subTemp(1:5), subx_v(1:5), subx_n(1:5), &
      & subDelOptDep(1:5), subOptDep(1:5),                            &
      & subRadLUwFlux(1:5), subRadLDwFlux(1:5), subFluxConv(1:4),     &
      & Press_ktp(1:3), Temp_ktp(1:3), OptDep_ktp(1:3),               &
      & RadLUwFlux_ktp(1:3), RadLDwFlux_ktp(1:3)

    ! 実行文 ; Executable statement
    !

    ! FluxConvergenceが負から正になるグリッドを探す. k = ktp
    ! find the grid that sign of flux convergence changes from negative to positive
    !

    ! 分割する層の決定 k= kpt-1  結pt+1の間に圏界面がある
    ! toropopause exists between the grid(kpt-1) and grid(ktp+1)
    !
    do k = kmax-1, 2, -1
      if ( z_FluxConv(k)<0.0 .and. z_FluxConv(k-1)>0.0 ) then
        ktp = k
        exit
      end if
    end do !k

    do k = 1, 3
      Press_ktp(k) = z_Press(k + ktp -2)
      Temp_ktp(k) = z_Temp(k + ktp -2)
      OptDep_ktp(k) = z_OptDep(k + ktp -2)
      RadLUwFlux_ktp(k) = z_RadLUwFlux(k + ktp -2)
      RadLDwFlux_ktp(k) = z_RadLDwFlux(k + ktp -2)
    end do !k

    ! 圏界面候補ktpが決まったところでループ
    do kkkk = 1, 100 ! 圏界面が決まるまでループ

      ! k = ktpを挟む2層を2つに分割

      ! 圧力を新しく設定
      subPress(1) = Press_ktp(1)
      subPress(2) = 0.5_8 * (Press_ktp(1) + Press_ktp(2))
      subPress(3) = Press_ktp(2)
      subPress(4) = 0.5_8 * (Press_ktp(2) + Press_ktp(3))
      subPress(5) = Press_ktp(3)

      ! 分割した層の温度は直線近似で与える
      subTemp(1) = Temp_ktp(1)
      subTemp(2) = 0.5_8 * (Temp_ktp(1) + Temp_ktp(2))
      subTemp(3) = Temp_ktp(2)
      subTemp(4) = 0.5_8 * (Temp_ktp(2) + Temp_ktp(3))
      subTemp(5) = Temp_ktp(3)

      ! 分割した層の光学的厚さを与える
      ! calculate optical depth in new layer
      ! in: Temp, ktp, x
      ! out: DelOptDep, OptDep

      ! 混合比xの再決定(温度上昇後、光学的厚さを求めるため)
      do k = 1, 5 ! in kkkk
        subx_v(k) = (p_v_0/Press_ktp(2)) * exp(-LatentHeat/(GasRUniv * Temp_ktp(2)))
        x_v_max = (p_v_0/subPress(k)) * exp(-LatentHeat/(GasRUniv * subTemp(k)))
        if(subx_v(k) > x_v_max) then
          subx_v(k) = x_v_max
        end if
        subx_n(k) = 1.0_8 - subx_v(k)
      end do ! k in kkkk

      ! 光学的厚さの計算
      ! calculate optical depth
      !
      do k = 2, 5 ! in kkkk
        subDelOptDep(k) = (                                               &
          & ( subPress(k)-subPress(k-1) ) / (                             &
          & ( 0.5_8 * (MolWtWet * subx_v(k-1) + MolWtDry * subx_n(k-1)) + &
          &   0.5_8 * (MolWtWet * subx_v(k)   + MolWtDry * subx_n(k)  ) ) &
          & * Grav )      ) * (                                           &
          & 0.5_8 * ( AbsCoefWetCom * MolWtWet * subx_v(k-1)            + &
          &           AbsCoefDryCom * MolWtDry * subx_n(k-1) )          + &
          & 0.5_8 * ( AbsCoefWetCom * MolWtWet * subx_v(k)              + &
          &           AbsCoefDryCom * MolWtDry * subx_n(k)   )            &
          &                    )
      end do ! k in kkkk

      subOptDep(1) = OptDep_ktp(1)
      subOptDep(2) = OptDep_ktp(1) + subDelOptDep(2)
      subOptDep(3) = OptDep_ktp(2)
      subOptDep(4) = OptDep_ktp(2) + subDelOptDep(4)
      subOptDep(5) = OptDep_ktp(3)

      ! 新しくできたグリッドのFluxを計算
      ! calculate the flux of new grid
      !
      subRadLUwFlux(5) = RadLUwFlux_ktp(3)

      do k = 5, 2, -1 ! in kkkk !!!!!FluxUp!!!!!

        pDN = ( abs(int(subPress(k) - subPress(k-1))) + 2 ) !!こちらは今回初めて
        ph = ( subPress(k-1) - subPress(k) ) / pDN
        tauh = ( subOptDep(k-1) - subOptDep(k) ) / pDN

        ! 台形法
        ! 台形量端の計算!!FluxUp
        subRadLUwFlux(k-1) = 0.5_8 * (-tauh) * (                                 &
          & func_Flux( StB * subTemp(k)**4  , subOptDep(k)  , subOptDep(k-1) ) + &
          & func_Flux( StB * subTemp(k-1)**4, subOptDep(k-1), subOptDep(k-1) )   )

        ! 台形間の計算!!FluxUp
        do kk = 1, pDN-1 !in k, kkkk
          TempRK  = subTemp(k) +                          &
            & ( subTemp(k-1) - subTemp(k) ) * (ph * kk) / &
            & ( subPress(k-1) - subPress(k) )             ! 直線

          OptDepDK  = subOptDep(k) +                          &
            & ( subOptDep(k-1) - subOptDep(k) ) * (ph * kk) / &
            & ( subPress(k-1) - subPress(k) )                 ! 直線
  
          subRadLUwFlux(k-1) = subRadLUwFlux(k-1) + &
            & (-tauh) * func_Flux( StB * TempRK**4, OptDepDK, subOptDep(k-1) )

        end do !kk in k, kkkk  !!FluxUp

        ! FluxUp収束するまで計算
        ! 
        preRadLFlux = subRadLUwFlux(k-1)

        do kkk = 1, 100 !in k, kkkk !精度が出るまでループ(上限は与える
          pDN = pDN * 2 !2倍ずつ与える
          ph = ( subPress(k-1) - subPress(k) ) / pDN
          tauh = ( subOptDep(k-1) - subOptDep(k) ) / pDN

          ! 台形法
          ! 台形量端の計算!!FluxUp
          subRadLUwFlux(k-1) = 0.5_8 * (-tauh) * (                               &
            & func_Flux(StB * subTemp(k)**4  , subOptDep(k)  , subOptDep(k-1)) + &
            & func_Flux(StB * subTemp(k-1)**4, subOptDep(k-1), subOptDep(k-1))   )

          ! 台形間の計算!!FluxUp
          do kk = 1, pDN-1 !in kkk, k, kkkk
            TempRK  = subTemp(k) +                          &
              & ( subTemp(k-1) - subTemp(k) ) * (ph * kk) / &
              & ( subPress(k-1) - subPress(k) )             ! 直線

            OptDepDK  = subOptDep(k) + &
              & ( subOptDep(k-1) - subOptDep(k) ) * (ph * kk) / &
              & ( subPress(k-1) - subPress(k) )                 ! 直線
  
           subRadLUwFlux(k-1) = subRadLUwFlux(k-1) + &
              & (-tauh) * func_Flux( StB * TempRK**4, OptDepDK, subOptDep(k-1) )

           end do !kk in kkk, k, kkkk

           ! 収束判定
           ! convergence test
           !
           err = abs( (subRadLUwFlux(k-1) - preRadLFlux ) / subRadLUwFlux(k-1) )
           If ( ( err < eps ) .or.                                                      &
             &  ( log(abs(subRadLUwFlux(k-1))) < minexponent(abs(subRadLUwFlux(k-1))) ) &
             & ) then

             If (log(abs(subRadLUwFlux(k-1))) < minexponent(abs(subRadLUwFlux(k-1))) ) then
               subRadLUwFlux(k-1) = 0.0_8
               !!$ print *, "FluxUp < 計算限界"
             end if
             !!$ print *, "err=", err, "kkk=", kkk
             exit
           end if
           preRadLFlux = subRadLUwFlux(k-1)
        end do !kkk in k, kkkk !====FluxUp収束完了====

        ! 1つ前の層までのFluxの寄与
        ! contribution of flux from pre layer
        !
        subRadLUwFlux(k-1) = subRadLUwFlux(k-1) + &
          & subRadLUwFlux(k) * exp(-1.5_8 * (subOptDep(k) - subOptDep(k-1)))

      end do !k, kkkk !!FluxUp計算終了!!===========FluxUp

      ! FluxDownの計算
      ! caluculate downward flux
      !
      subRadLDwFlux(1) = RadLDwFlux_ktp(1)

      do k = 2, 5 ! in kkkk  !!!FluxDown

        pDN = ( abs(int(subPress(k) - subPress(k-1))) + 2 ) !!こちらは今回初めて
        ph = ( subPress(k-1) - subPress(k) ) / pDN
        tauh = ( subOptDep(k-1) - subOptDep(k) ) / pDN

        ! 台形法
        ! 台形量端の計算 !!FluxDown
        subRadLDwFlux(k) = 0.5_8 * (-tauh) * (                               &
          & func_Flux(StB * subTemp(k)**4  , subOptDep(k)  , subOptDep(k)) + &
          & func_Flux(StB * subTemp(k-1)**4, subOptDep(k), subOptDep(k-1))   )

        ! 台形間の計算 !!FluxDown
        do kk = 1, pDN-1 ! in k, kkkk
          TempRK  = subTemp(k) +                          &
            & ( subTemp(k-1) - subTemp(k) ) * (ph * kk) / & 
            & ( subPress(k-1) - subPress(k) )             ! 直線

          OptDepDK  = subOptDep(k) +                          &
            & ( subOptDep(k-1) - subOptDep(k) ) * (ph * kk) / &
            & ( subPress(k-1) - subPress(k) )                 ! 直線

          subRadLDwFlux(k) = subRadLDwFlux(k) + &
            & (-tauh) * func_Flux(StB * TempRK**4, subOptDep(k), OptDepDK)
        end do !kk in k, kkkk !!FluxDown

        ! FluxDown収束するまで計算
        !
        preRadLFlux = subRadLDwFlux(k)

        do kkk = 1, 100 !in k, kkkk !精度が出るまでループ(上限は与える
          pDN = pDN * 2 !2倍ずつ与える
          ph = ( subPress(k-1) - subPress(k) ) / pDN
          tauh = ( subOptDep(k-1) - subOptDep(k) ) / pDN

          ! 台形法
          ! 台形量端の計算 !!FluxDown
          subRadLDwFlux(k) = 0.5_8 * (-tauh) * (                               &
            & func_Flux(StB * subTemp(k)**4  , subOptDep(k)  , subOptDep(k)) + &
            & func_Flux(StB * subTemp(k-1)**4, subOptDep(k), subOptDep(k-1))   )

          ! 台形間の計算 !!FluxDown
          do kk = 1, pDN-1 ! in kkk, k, kkkk
            TempRK  = subTemp(k) +                          &
              & ( subTemp(k-1) - subTemp(k) ) * (ph * kk) / &
              & ( subPress(k-1) - subPress(k) )             ! 直線

            OptDepDK  = subOptDep(k) +                          &
              & ( subOptDep(k-1) - subOptDep(k) ) * (ph * kk) / &
              & ( subPress(k-1) - subPress(k) )                 ! 直線
 
            subRadLDwFlux(k) = subRadLDwFlux(k) + &
              & (-tauh) * func_Flux( StB * TempRK**4, subOptDep(k), OptDepDK )
          end do !kk in kkk, k, kkkk !!FluxDown

          ! 収束判定
          ! convergence test
          !
          err = abs( (subRadLDwFlux(k) - preRadLFlux) / subRadLDwFlux(k) )
          If ( ( err < eps ) .or.                                                  &
            &  ( log(abs(subRadLDwFlux(k))) < minexponent(abs(subRadLDwFlux(k))) ) &
            & ) then

            If ( log(abs(subRadLDwFlux(k))) < minexponent(abs(subRadLDwFlux(k))) ) then
              subRadLDwFlux(k) = 0.0_8
              !!$ print *, "FluxDown < 計算限界"
            end if
            !!$ print *, "err=", err, "kkk=", kkk
            exit
          end if
          preRadLFlux = subRadLDwFlux(k)
        end do !kkk in k, kkkk

        ! 1つ前の層までのFluxの寄与
        ! contribution of flux from pre layer
        !
        subRadLDwFlux(k) = subRadLDwFlux(k) + &
          & subRadLDwFlux(k-1) * exp(-1.5_8 * (subOptDep(k) - subOptDep(k-1)))

      end do !k in kkkk ! FluxDown成層圏計算終了

      ! FluxConvergenceを計算
      !
      ! FluxConvergenceの計算と境界面の判定
      !
      subFluxConv(:) = 0.0_8
      do k = 1, 4 ! in kkkk
        subFluxConv(k) = subRadLUwFlux(k+1) - subRadLUwFlux(k) - &
          &            ( subRadLDwFlux(k+1) - subRadLDwFlux(k) )
      end do ! k in kkkk

      do k = 4, 2, -1 ! in kkkk
        if( subFluxConv(k) < 0.0 .and. subFluxConv(k-1) > 0.0 ) then
          numktp = k
          exit
        end if
      end do ! k in kkkk

      do k = 1, 3 ! in kkkk
        Press_ktp(k) = subPress(k + numktp - 2)
        Temp_ktp(k) = subTemp(k + numktp - 2)
        OptDep_ktp(k) = subOptDep(k + numktp - 2)
        RadLUwFlux_ktp(k) = subRadLUwFlux(k + numktp - 2)
        RadLDwFlux_ktp(k) = subRadLDwFlux(k + numktp - 2)
      end do ! k in kkkk

      select case(numktp)
        case(2)
          FluxConv_minus = FluxConv_minus + subFluxConv(3) + subFluxConv(4)
        case(3)
          FluxConv_plus = FluxConv_plus + subFluxConv(1)
          FluxConv_minus = FluxConv_minus + subFluxConv(4)
        case(4)
          FluxConv_plus = FluxConv_plus + subFluxConv(1) + subFluxConv(2)
      end select

      if( abs( subFluxConv(numktp  ) ) < 0.000001_8 .and. &
        &      subFluxConv(numktp-1)   < 0.0000001_8      ) then
        FluxConv_plus  = FluxConv_plus  + subFluxConv(numktp-1)
        FluxConv_minus = FluxConv_minus + subFluxConv(numktp)
        ktpRN = kkkk
        exit
      end if

    end do !kkkk !圏界面ktpが決まるまでループ ! 圏界面決定

    if( z_Press(ktp) > Press_ktp(2) ) then
      FluxConv_plus = FluxConv_plus - z_FluxConv(ktp-1)
    else
      FluxConv_minus = FluxConv_minus - z_FluxConv(ktp)
    end if

    val_tropo(1) = Press_ktp(2)
    val_tropo(2) = Temp_ktp(2)
    val_tropo(3) = OptDep_ktp(2)
    val_tropo(4) = RadLUwFlux_ktp(2)
    val_tropo(5) = RadLDwFlux_ktp(2)
    val_tropo(6) = FluxConv_plus
    val_tropo(7) = FluxConv_minus

  end subroutine calc_Tropopause

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

end module subprogs



program NakajimaTest
  !
  != Nakajima et al., 1992の検証プログラム
  !
  != Program for Nakajima et al,. 1992
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! Nakajima et al.1992の計算
  !
  ! 1. 湿潤断熱減率に従って上空まで温度、Fluxを計算
  !     湿潤断熱減率の計算: RungeKutta
  !     光学的厚さの計算: 台形
  ! 2. 対流圏界面の計算(FluxConvergenceの正負判定)
  ! 3. 対流圏界面を挿入
  ! 4. 圏界面水蒸気混合比で飽和
  ! 5. 放射平衡計算
  !    5.1 加熱
  !        FluxConvergenceがグリッドごとに正負(負正はOK)を繰り返すときには温度を変えない
  !    5.2 光学的厚さの決定
  !    5.3 Flux, FluxConvergenceの計算
  !        Total FluxConvergenceの微分が小さくなったらstep.6へ
  !        Total FluxConvergenceの値が小さくなったら計算終了
  ! 6. 圏界面移動
  ! 7. 5-6をループ
  !
  !== NAMELIST
  !
  ! NAMELIST#conditions_nml
  !

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

  ! 1次元放射対流平衡放射伝達モデル(Nakajima et al., 1992)
  ! Radiative transfer with One-Dimensional Radiative-Convective Equilibrium Model
  !
  use subprogs

  ! ヒストリデータ出力
  ! History data output
  !
  use gtool_history

  ! 宣言文 ; Declaration statements
  !
  implicit none

  ! 物理定数, 計算条件
  ! Physical constants, condition of calculation
  !
  real(8) :: GasRUniv = 8.314_8
                              ! $ R^{*} $ [J K-1 mol-1].
                              ! 普遍気体定数.
                              ! Universal gas constant
  real(8) :: StB = 5.67e-8_8
                              ! $ \sigma_{SB} $ [W m-2 K-4]. 
                              ! ステファンボルツマン定数. 
                              ! Stefan-Boltzmann constant
  real(8) :: Grav = 9.8_8
                              ! $ g $ [m s-2]. 
                              ! 重力加速度. 
                              ! Gravitational acceleration
  real(8) :: MolWtDry = 18.0e-3_8
                              ! $ M $ [kg mol-1]. 
                              ! 乾燥大気の平均分子量. 
                              ! Mean molecular weight of dry air
  real(8) :: CpDry
                              ! $ C_p $ [J mol-1 K-1]. 
                              ! 乾燥大気の定圧比熱. 
                              ! Specific heat of air at constant pressure
  real(8) :: MolWtWet = 18.0e-3_8
                              ! $ M_v $ [kg mol-1]. 
                              ! 凝結成分の平均分子量. 
                              ! Mean molecular weight of condensible elements
  real(8) :: CpWet
                              ! $ C_v $ [J mol-1 K-1] . 
                              ! 凝結成分の定圧比熱. 
                              ! Specific heat of condensible elements at constant pressure
  real(8) :: LatentHeat = 43655.0_8
                              ! $ L $ [J mol-1] . 
                              ! 凝結の潜熱. 
                              ! Latent heat of condensation
  real(8) :: p_v_0 = 1.4e+11_8
                              ! $ p*_0 $ [Pa] .
                              ! 飽和水蒸気曲線の係数.
                              ! Constant for the water vapor saturation curve
  real(8) :: AbsCoefWetCom = 0.01_8
                              ! $ k_v $ [m2 kg-1] .
                              ! 凝結成分の吸収係数.
                              ! Absorption coefficient of condensable component
  real(8) :: AbsCoefDryCom = 0.0_8
                              ! $ k_n $ [m2 kg-1] .
                              ! 非凝結成分の吸収係数.
                              ! Absorption coefiicient of noncondensable component
  real(8) :: SurfTemp           
                              ! $ T_s $ [K] .
                              ! 地表面温度.
                              ! Surface temperature
  real(8) :: Humidity = 1.0_8
                              ! $ h $ [1] .
                              ! 相対湿度.
                              ! Relative humidity
  real(8) :: p_n_s = 1.0e+5_8
                              ! $ pn_0 $ [Pa] .
                              ! 非凝結大気の地表での圧力
                              ! Amount of noncondensable component at the bottom of atmosphere
  real(8) :: sigma_min
                              ! $ sigma_min $ [1] .
                              ! 大気上端での圧力の大きさ(地表に対する割合)
                              ! Ratio of pressure at the top to the bottom
  integer :: kmax
                              ! 格子点設定: 鉛直総数
                              ! Grid settings: Number of vertical level
  integer :: changeP = 2
                              ! 格子点設定: gridの並べ方(圧力等間隔: 1, ln(圧力)等間隔: 2)
                              ! Grid settings: Arrangement of grid
                              ! (equal interval of 1. pressure or 2. ln(pressure))
  real(8) :: DelTime
                              ! $ dt $ [s] .
                              ! タイムステップ.
                              ! Time step
  real(8) :: eps
                              ! 計算精度(温度、光学的厚さ、フラックス).
                              ! Resolution of calculation (temperature, optical depth, flux) 
  real(8) :: eps_FluxConv
                              ! 計算精度(放射平衡の収束)
                              ! Resolution of caluculation (radiative equilibrium)

  ! 変数
  ! variables
  !
  real(8), allocatable :: z_Press(:)
                              ! $ p $ [Pa].
                              ! 圧力
                              ! pressure
  real(8), allocatable :: z_Temp(:)
                              ! $ T $ [K].
                              ! 温度
                              ! temperature
  real(8), allocatable :: z_OptDep(:)
                              ! $ tau $ [1].
                              ! 光学的厚さ
                              ! optical depth
  real(8), allocatable :: z_RadLUwFlux(:)
                              ! $ F_up $ [W m-2].
                              ! 上向き長波フラックス
                              ! upward longwave flux
  real(8), allocatable :: z_RadLDwFlux(:)
                              ! $ F_down $ [W m-2].
                              ! 下向き長波フラックス
                              ! downward longwave flux
  real(8), allocatable :: z_FluxConv(:)
                              ! flux convergence [W m-2].
  real(8), allocatable :: z_Press_tp(:)
                              ! $ p $ [Pa].
                              ! 圧力(対流圏界面を含む)
                              ! pressure(+ tropopause)
  real(8), allocatable :: z_Temp_tp(:)
                              ! $ T $ [K].
                              ! 温度(対流圏界面を含む)
                              ! temperature(+ tropopause)
  real(8), allocatable :: z_OptDep_tp(:)
                              ! $ tau $ [1].
                              ! 光学的厚さ(対流圏界面を含む)
                              ! optical depth(+ tropopause)
  real(8), allocatable :: z_RadLUwFlux_tp(:)
                              ! $ F_up $ [W m-2].
                              ! 上向き長波フラックス(対流圏界面を含む)
                              ! upward longwave flux(+ tropopause)
  real(8), allocatable :: z_RadLDwFlux_tp(:)
                              ! $ F_down $ [W m-2].
                              ! 下向き長波フラックス(対流圏界面を含む)
                              ! downward longwave flux(+ tropopause)
  real(8), allocatable :: z_FluxConv_tp(:)
                              ! flux convergence(+ tropopause) [W m-2].
  real(8) :: val_tropo(1:7)   ! 対流圏界面の物理量
                              ! values of toropopause
  real(8) :: preOptDep_tp     ! 
                              ! previous value of tropopause optical depth
  real(8) :: preSumFluxConv   !
                              ! previous value of summation of flux convergence
  real(8) :: SumFluxConv      !
                              ! summation of flux congergence [W m-2].
  real(8) :: DelSumFluxConv   !
                              ! differentiation of sum of flux convergence

  ! 作業変数
  ! work variables
  !
  real(8) :: com_g            ! 圧力グリッドを計算するための係数
                              ! Coeficient for pressure grid
  integer :: ktp              ! 対流圏界面のグリッド
                              ! Grid number of toropopause
  integer, allocatable :: knum(:)
                              ! 積分計算のための作業変数
                              ! Work variables 
  integer, allocatable :: kn(:)
                              ! netCDF書き出しのためのk次元
                              ! k-dimension for netCDF output
  integer :: k, kk, kkk       ! 鉛直方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in vertical direction
  integer :: count            ! 収束までの計算回数をカウントする作業変数
                              ! Work variable for radiative eq. caluclation loop
  integer :: fi = 10          ! namelistファイルの装置変数
                              ! device number for namelist file
  integer :: fo = 11          ! 出力ファイルの装置変数
                              ! device number for output file
  logical :: OutText = .true.         
                              ! textフォーマットの書き出し設定
                              ! text format output(yes: .true., no: .false.)
  logical :: OutNetCDF_initial = .true.
                              ! netCDFフォーマットの書き出し設定(moist ad. data)
                              ! netCDF format output(yes: .true., no: .false.)
  logical :: OutNetCDF_run  = .true.
                              ! netCDFフォーマットの書き出し設定(running results)
                              ! netCDF format output(yes: .true., no: .false.)
  logical :: OutNetCDF_final = .true.
                              ! netCDFフォーマットの書き出し設定(final results)
                              ! netCDF format output(yes: .true., no: .false.)

  ! NAMELIST 変数群
  ! NAMELIST group name
  !
  namelist /conditions_nml/ &
    & MolWtDry,             &
    & CpDry,                &
    & MolWtWet,             &
    & CpWet,                &
    & AbsCoefWetCom,        &
    & AbsCoefDryCom,        &
    & SurfTemp,             &
    & Humidity,             &
    & p_n_s,                &
    & sigma_min,            &
    & kmax,                 &
    & changeP,              &
    & DelTime,              &
    & eps,                  &
    & eps_FluxConv

  ! 実行文 ; Executable statement
  !

  ! NAMELIST変数の読み込み
  ! Reading NAMELIST variables
  !
  open (fi, file='nakajima_RadiatConvEqModel.nml')
  read (fi, nml=conditions_nml)
  close(fi)

  ! 変数の割付
  ! Allocation of variables
  !
  allocate( z_Press        (1:kmax  ) )
  allocate( z_Temp         (1:kmax  ) )
  allocate( z_OptDep       (1:kmax  ) )
  allocate( z_RadLUwFlux   (1:kmax  ) )
  allocate( z_RadLDwFlux   (1:kmax  ) )
  allocate( z_FluxConv     (1:kmax-1) )

  allocate( z_Press_tp     (1:kmax+1) )
  allocate( z_Temp_tp      (1:kmax+1) )
  allocate( z_OptDep_tp    (1:kmax+1) )
  allocate( z_RadLUwFlux_tp(1:kmax+1) )
  allocate( z_RadLDwFlux_tp(1:kmax+1) )
  allocate( z_FluxConv_tp  (1:kmax  ) )

  allocate( knum           (1:kmax  ) )
  allocate( kn             (1:kmax+1) )

  ! textデータによる書き出しの設定
  ! Setting of data output in text format
  !
  if( OutText ) then
    open(fo, file = 'output_RCEM.d')
    write(fo, nml=conditions_nml)
  end if

  ! netCDF書き出しの設定(湿潤断熱減率計算結果)
  ! Setting of data output in netCDF(results of moist ad. )
  !
  if( OutNetCDF_initial ) then
    call HistoryCreate(                              &
      & file='RCEMnetCDF_initial.nc'               , &
      & title='RCEMnetCDF_initial'                 , &
      & source='Nakajima et al, 1992'              , &
      & institution='GFD_Dennou Club davis project', &
      & dims=(/'p'/), dimsizes=(/kmax/)            , &
      & longnames=(/'p-coordinate'/)               , &
      & units=(/'p'/)                                )

    call HistoryAddVariable(               &
      & varname='temp', dims=(/'p'/)     , &
      & longname='temperature'           , &
      & units='K', xtype='double'          )
    call HistoryAddVariable(               &
      & varname='optdep', dims=(/'p'/)   , &
      & longname='optical depth'         , &
      & units='1', xtype='double'          )
    call HistoryAddVariable(               &
      & varname='fluxup', dims=(/'p'/)   , &
      & longname='longwave upward flux'  , &
      & units='W m-2', xtype='double'      )
    call HistoryAddVariable(               &
      & varname='fluxdown', dims=(/'p'/) , &
      & longname='longwave downward flux', &
      & units='W m-2', xtype='double'      )
    call HistoryAddVariable(               &
      & varname='fluxconv', dims=(/'p'/) , &
      & longname='flux convergence'      , &
      & units='W m-2', xtype='double'      )
  end if

  ! 地表面温度の設定
  ! Setting of surface temperature
  !
  z_Temp(:) = 0.0_8
  z_Temp(kmax) = SurfTemp

  ! 圧力の鉛直分布の設定
  ! Setting of pressure vertical profile
  !
  p_v_0 = p_v_0 * Humidity
  z_Press(kmax) = p_n_s + p_v_0 * exp( -LatentHeat/(GasRUniv*SurfTemp))

  if(changeP == 2) then    ! ln(pressure): 等間隔
    com_g = (-log(sigma_min))/(kmax-1)
    do k = 1, kmax-1
      z_Press(k) = z_Press(kmax) * sigma_min * exp(com_g * (k - 1))
    end do
  end if

  if(changeP == 1) then    ! pressure : 等間隔
    z_Press(1) = z_Press(kmax) * sigma_min
    do k = 2, kmax
      z_Press(k) = z_Press(1) + &
        & ((z_Press(kmax) - z_Press(1)) * (k-1))/(kmax - 1)
    end do
  end if

  ! 湿潤断熱減率の計算
  ! Calculation of moist adiabatic lapse late
  !
  call calc_AdiabLapse(                                              &
    & GasRUniv, StB, LatentHeat, Grav, MolWtDry, MolWtWet          , & ! (in )
    & AbsCoefDryCom, AbsCoefWetCom, p_v_0, CpWet, CpDry            , & ! (in )
    & eps, kmax, z_Press                                           , & ! (in )
    & knum, z_Temp, z_OptDep, z_RadLUwFlux, z_RadLDwFlux, z_FluxConv & ! (out)
    & )
  ! netCDF書き出し(湿潤断熱減率計算結果)
  ! data output in netCDF(results of moist ad. )
  !
  if( OutNetCDF_initial ) then !!Gtool output initial result
    call HistoryPut('p'       ,z_Press     )
    call HistoryPut('temp'    ,z_Temp      )
    call HistoryPut('optdep'  ,z_OptDep    )
    call HistoryPut('fluxup'  ,z_RadLUwFlux)
    call HistoryPut('fluxdown',z_RadLDwFlux)
    call HistoryPut('fluxconv',z_FluxConv  )
    call HistoryClose
  end if

  ! netCDF書き出しの設定(放射平衡計算途中結果)
  ! Setting of data output in netCDF(running history of radiative eq. calc. )
  !
  if( OutNetCDF_run ) then !!Gtool output running
    do k = 1, kmax+1
      kn(k) = k
    end do

    call HistoryCreate(                              &
      & file='RCEMnetCDF_run.nc'                   , &
      & title='RCEMnetCDF_run'                     , &
      & source='Nakajima et al, 1992'              , &
      & institution='GFD_Dennou Club davis project', &
      & dims=(/'k','t'/), dimsizes=(/kmax+1,0/)    , &
      & longnames=(/'k-coordinate','time        '/), &
      & units=(/'1','s'/)                          , &
      & origin=real(0.0), interval=real(DelTime)     )

    call HistoryPut('k',kn)

    call HistoryAddVariable(                  &
      & varname='press', dims=(/'k','t'/)   , &
      & longname='pressure'                 , &
      & units='N/m2', xtype='double'          )
    call HistoryAddVariable(                  &
      & varname='temp', dims=(/'k','t'/)    , &
      & longname='temperature'              , &
      & units='K', xtype='double'             )
    call HistoryAddVariable(                  &
      & varname='optdep', dims=(/'k','t'/)  , &
      & longname='optical depth'            , &
      & units='1', xtype='double'             )
    call HistoryAddVariable(                  &
      & varname='fluxup', dims=(/'k','t'/)  , &
      & longname='longwave upward flux'     , &
      & units='W m-2', xtype='double'         )
    call HistoryAddVariable(                  &
      & varname='fluxdown', dims=(/'k','t'/), &
      & longname='longwave downward flux'   , &
      & units='W m-2', xtype='double'         )
    call HistoryAddVariable(                  &
      & varname='fluxconv', dims=(/'k','t'/), &
      & longname='flux convergence'         , & 
      & units='W m-2', xtype='double'         )
  end if

  ! 圏界面の計算
  ! calculate tropopause
  !
  !in: z_Press, z_Temp, z_OptDep, z_RadLUwFlux, z_RadLDwFlux, z_FluxConv
  !out: ktp, val_tropo

  do k = kmax-2, 3, -1 !始めをkmax-2にしています。1番下の層を外しています。
    if ( z_FluxConv(k)<0 .and. z_FluxConv(k-1)>0 ) then
      ktp = k
      exit
    end if
  end do ! k

  call calc_Tropopause(                                                  &
    & GasRUniv, StB, LatentHeat, Grav, MolWtDry, MolWtWet,               & ! (in )
    & AbsCoefDryCom, AbsCoefWetCom, p_v_0, eps, kmax,                    & ! (in )
    & z_Press, z_Temp, z_OptDep, z_RadLUwFlux, z_RadLDwFlux, z_FluxConv, & ! (in )
    & val_tropo                                                          & ! (out)
    & )

  ! 圏界面挿入
  ! insert grid of tropopause
  !
  do k = 1, ktp-1
    z_Press_tp      (k) = z_Press      (k)
    z_Temp_tp       (k) = z_Temp       (k)
    z_OptDep_tp     (k) = z_OptDep     (k)
    z_RadLUwFlux_tp (k) = z_RadLUwFlux (k)
    z_RadLDwFlux_tp (k) = z_RadLDwFlux (k)
    z_FluxConv_tp   (k) = z_FluxConv   (k)
  end do

  if(val_tropo(1) > z_Press(ktp)) then ! 始めの圏界面より実際は下にある場合
     z_Press_tp      (ktp  ) = z_Press      (ktp)
     z_Press_tp      (ktp+1) = val_tropo    (1  )
     z_Temp_tp       (ktp  ) = z_Temp       (ktp)
     z_Temp_tp       (ktp+1) = val_tropo    (2  )
     z_OptDep_tp     (ktp  ) = z_OptDep     (ktp)
     z_OptDep_tp     (ktp+1) = val_tropo    (3  )
     z_RadLUwFlux_tp (ktp  ) = z_RadLUwFlux (ktp)
     z_RadLUwFlux_tp (ktp+1) = val_tropo    (4  )
     z_RadLDwFlux_tp (ktp  ) = z_RadLDwFlux (ktp)
     z_RadLDwFlux_tp (ktp+1) = val_tropo    (5  )
     z_FluxConv_tp   (ktp  ) = val_tropo    (6  )
     z_FluxConv_tp   (ktp+1) = val_tropo    (7  )
   else                                 ! 始めの圏界面より実際は上にある場合
     z_Press_tp      (ktp  ) = val_tropo    (1  )
     z_Press_tp      (ktp+1) = z_Press      (ktp)
     z_Temp_tp       (ktp  ) = val_tropo    (2  )
     z_Temp_tp       (ktp+1) = z_Temp       (ktp)
     z_OptDep_tp     (ktp  ) = val_tropo    (3  )
     z_OptDep_tp     (ktp+1) = z_OptDep     (ktp)
     z_RadLUwFlux_tp (ktp  ) = val_tropo    (4  )
     z_RadLUwFlux_tp (ktp+1) = z_RadLUwFlux (ktp)
     z_RadLDwFlux_tp (ktp  ) = val_tropo    (5  )
     z_RadLDwFlux_tp (ktp+1) = z_RadLDwFlux (ktp)
     z_FluxConv_tp   (ktp-1) = val_tropo    (6  )
     z_FluxConv_tp   (ktp  ) = val_tropo    (7  )
     z_FluxConv_tp   (ktp+1) = z_FluxConv   (ktp)
  end if

  do k = ktp+2, kmax+1
    z_Press_tp      (k) = z_Press      (k-1)
    z_Temp_tp       (k) = z_Temp       (k-1)
    z_OptDep_tp     (k) = z_OptDep     (k-1)
    z_RadLUwFlux_tp (k) = z_RadLUwFlux (k-1)
    z_RadLDwFlux_tp (k) = z_RadLDwFlux (k-1)
  end do

  do k = ktp+2, kmax
    z_FluxConv_tp(k) = z_FluxConv(k-1)
  end do

  ! 成層圏を圏界面の水蒸気混合比で飽和
  ! saturate stratosphere
  !
  do k = kmax-1, 3, -1
    if ( z_FluxConv_tp(k)<0 .and. z_FluxConv_tp(k-1)>0 ) then
      ktp = k
      exit
    end if
  end do

  preOptDep_tp = z_OptDep_tp(ktp)
  do k = 1, ktp !ktp-1から
    z_Temp_tp(k) = z_Temp_tp(ktp) / (                    &
      & 1.0_8 - z_Temp_tp(ktp) * (GasRUniv/LatentHeat) * &
      & log(z_Press_tp(k)/z_Press_tp(ktp))               &
      &                              )
    z_OptDep_tp(k) = ( z_Press_tp(k) / z_Press_tp(ktp) ) * (                              &
      & AbsCoefWetCom * MolWtWet * p_v_0 * exp(-LatentHeat/(GasRUniv * z_Temp_tp(k)))   + &
      & AbsCoefDryCom * MolWtDry * (                                                      &
      & z_Press_tp(ktp) - p_v_0 * exp(-LatentHeat/(GasRUniv * z_Temp_tp(k)))              &
      &                             )                      ) / ( (                        &
      & MolWtWet * p_v_0 * exp(-LatentHeat/(GasRUniv * z_Temp_tp(k)))                   + &
      & MolWtDry * (z_Press_tp(ktp) - p_v_0 * exp(-LatentHeat/(GasRUniv * z_Temp_tp(k)))) &
      &                                                           ) * Grav )
  end do

  ! 全光学的厚さの計算
  ! calculate optical depth
  !
  do k = ktp+1, kmax+1
    z_OptDep_tp(k) = z_OptDep_tp(k) + ( z_OptDep_tp(ktp) - preOptDep_tp )
  end do

  ! 放射平衡計算
  ! calculate radiative equilibrium in stratosphere
  !
  count = 0
  do kkk = 1, 10000
    call calc_RaiseTemp(                                     &
      & GasRUniv, StB, LatentHeat, Grav, MolWtDry, MolWtWet, &
      & AbsCoefDryCom, AbsCoefWetCom, p_v_0, CpWet, CpDry  , &
      & eps, DelTime, kmax+1, ktp, knum                    , &
      & z_Press_tp, z_Temp_tp, z_OptDep_tp                 , &
      & z_RadLUwFlux_tp, z_RadLDwFlux_tp, z_FluxConv_tp      )

    do k = 1, ktp-1 ! in kkk
      preSumFluxConv = preSumFluxConv + z_FluxConv_tp(k)
    end do ! k in kkk

    do kk = 1, 100000 ! in kkk
      call calc_RaiseTemp(                                     &
        & GasRUniv, StB, LatentHeat, Grav, MolWtDry, MolWtWet, &
        & AbsCoefDryCom, AbsCoefWetCom, p_v_0, CpWet, CpDry  , &
        & eps, DelTime, kmax+1, ktp, knum                    , &
        & z_Press_tp, z_Temp_tp, z_OptDep_tp                 , &
        & z_RadLUwFlux_tp, z_RadLDwFlux_tp, z_FluxConv_tp      )
      count = count + 1

      SumFluxConv = 0.0_8
      do k = 1, ktp-1 ! in kk, kkk
        SumFluxConv = SumFluxConv + z_FluxConv_tp(k)
      end do ! k in kk, kkk

      if ( OutText ) then
        write(fo, *) "z_Press_tp z_Temp_tp z_OptDep_tp z_RadLUwFlux_tp z_RadLDwFlux_tp z_FluxConv_tp"
        do k=1, kmax+1 ! in kk, kkk
          write(fo, *) z_Press_tp(k), z_Temp_tp(k), z_OptDep_tp(k), &
            & z_RadLUwFlux_tp(k), z_RadLDwFlux_tp(k), z_FluxConv_tp(k)
        end do ! k in kk, kkk
      end if

      ! netCDF書き出し(放射平衡計算途中結果)
      ! data output in netCDF(running history of radiative eq. calc. )
      !
      if( OutNetCDF_run ) then
        call HistoryPut('press'   ,z_Press_tp     )
        call HistoryPut('temp'    ,z_Temp_tp      )
        call HistoryPut('optdep'  ,z_OptDep_tp    )
        call HistoryPut('fluxup'  ,z_RadLUwFlux_tp)
        call HistoryPut('fluxdown',z_RadLDwFlux_tp)
        call HistoryPut('fluxconv',z_FluxConv_tp  )
      end if

      DelSumFluxConv = (SumFluxConv - preSumFluxConv) / SumFluxConv
      if(abs( DelSumFluxConv ) < eps_FluxConv ) then
        if ( OutText ) then
          write(fo, *) count, "kk=", kk, "DelSumFluxConv=", DelSumFluxConv, &
            & "SumFluxConv=", SumFluxConv                                 , &
            & "differencital FluxConvergence is enough small"
        end if
        exit
      else
        if ( OutText ) then
          write(fo, *) count, "kk=", kk, "DelSumFluxConv=", DelSumFluxConv, &
            & "SumFluxConv=", SumFluxConv                                 , &
            & "one more loop!"
        end if
      end if

         preSumFluxConv = SumFluxConv

    end do ! kk in kkk

    ! 放射平衡計算の収束判定
    ! convergence test of radiative equilibrium
    !
    if( abs( SumFluxConv ) < eps_FluxConv ) then
      write(fo, *), "End of Calculation, SumFluxConv=", SumFluxConv
      exit
    end if

    ! 圏界面移動
    ! move tropopause
    !
    call calc_Tropopause(                                    &
      & GasRUniv, StB, LatentHeat, Grav, MolWtDry, MolWtWet, &
      & AbsCoefDryCom, AbsCoefWetCom, p_v_0, eps, kmax+1   , &
      & z_Press_tp, z_Temp_tp, z_OptDep_tp                 , &
      & z_RadLUwFlux_tp, z_RadLDwFlux_tp, z_FluxConv_tp    , &
      & val_tropo                                            )

    z_Press_tp      (ktp  ) = val_tropo (1)
    z_Temp_tp       (ktp  ) = val_tropo (2)
    z_OptDep_tp     (ktp  ) = val_tropo (3)
    z_RadLUwFlux_tp (ktp  ) = val_tropo (4)
    z_RadLDwFlux_tp (ktp  ) = val_tropo (5)
    z_FluxConv_tp   (ktp-1) = val_tropo (6)
    z_FluxConv_tp   (ktp  ) = val_tropo (7)

    do while(val_tropo(1) > z_Press_tp(ktp+1))
      ktp = ktp + 1
      z_Press_tp(ktp-1) = z_Press_tp(ktp)
      z_Temp_tp(ktp-1) = z_Temp_tp(ktp)
      z_OptDep_tp(ktp-1) = z_OptDep_tp(ktp-1)
      z_RadLUwFlux_tp(ktp-1) = z_RadLUwFlux_tp(ktp)
      z_RadLDwFlux_tp(ktp-1) = z_RadLDwFlux_tp(ktp)
      z_FluxConv_tp(ktp-2) = z_FluxConv_tp(ktp-1)
   
      z_Press_tp(ktp) = val_tropo(1)
      z_Temp_tp(ktp) = val_tropo(2)
      z_OptDep_tp(ktp) = val_tropo(3)
      z_RadLUwFlux_tp(ktp) = val_tropo(4)
      z_RadLDwFlux_tp(ktp) = val_tropo(5)
      z_FluxConv_tp(ktp-1) = val_tropo(6)
      z_FluxConv_tp(ktp) = val_tropo(7)  
    end do

    do while( val_tropo(1) < z_Press_tp(ktp-1) )
      ktp = ktp -1
      z_Press_tp      (ktp+1) = z_Press_tp      (ktp  )
      z_Temp_tp       (ktp+1) = z_Temp_tp       (ktp  )
      z_OptDep_tp     (ktp+1) = z_OptDep_tp     (ktp-1)
      z_RadLUwFlux_tp (ktp+1) = z_RadLUwFlux_tp (ktp  )
      z_RadLDwFlux_tp (ktp+1) = z_RadLDwFlux_tp (ktp  )
      z_FluxConv_tp   (ktp+1) = z_FluxConv_tp   (ktp  )

      z_Press_tp      (ktp  ) = val_tropo (1)
      z_Temp_tp       (ktp  ) = val_tropo (2)
      z_OptDep_tp     (ktp  ) = val_tropo (3)
      z_RadLUwFlux_tp (ktp  ) = val_tropo (4)
      z_RadLDwFlux_tp (ktp  ) = val_tropo (5)
      z_FluxConv_tp   (ktp-1) = val_tropo (6)
      z_FluxConv_tp   (ktp  ) = val_tropo (7)
    end do
  end do ! kkk

  ! netCDFファイルの終了(放射平衡計算途中結果)
  ! file close of output data in netCDF(running history of radiative eq. calc. )
  !
  if( OutNetCDF_run ) then
    call HistoryClose
  end if

  ! netCDF書き出し(最終計算結果)
  ! data output in netCDF(final result of radiative eq. calc. )
  !
  if( OutNetCDF_final ) then
    call HistoryCreate(                              &
      & file='RCEMnetCDF_final.nc'                 , &
      & title='RCEMnetCDF_final'                   , &
      & source='Nakajima et al, 1992'              , &
      & institution='GFD_Dennou Club davis project', &
      & dims=(/'p'/), dimsizes=(/kmax+1/)          , &
      & longnames=(/'p-coordinate'/)               , &
      & units=(/'p'/)                                )
    
    call HistoryAddVariable(               &
      & varname='temp', dims=(/'p'/)     , &
      & longname='temperature'           , &
      & units='K', xtype='double'          )
    call HistoryAddVariable(               &
      & varname='optdep', dims=(/'p'/)   , &
      & longname='optical depth'         , &
      & units='1', xtype='double'          )
    call HistoryAddVariable(               &
      & varname='fluxup', dims=(/'p'/)   , &
      & longname='longwave upward flux'  , &
      & units='W m-2', xtype='double'      )
    call HistoryAddVariable(               &
      & varname='fluxdown', dims=(/'p'/) , &
      & longname='longwave downward flux', &
      & units='W m-2', xtype='double'      )
    call HistoryAddVariable(               &
      & varname='fluxconv', dims=(/'p'/) , &
      & longname='flux convergence'      , &
      & units='W m-2', xtype='double'      )

    call HistoryPut('p'       ,z_Press_tp     )
    call HistoryPut('temp'    ,z_Temp_tp      )
    call HistoryPut('optdep'  ,z_OptDep_tp    )
    call HistoryPut('fluxup'  ,z_RadLUwFlux_tp)
    call HistoryPut('fluxdown',z_RadLDwFlux_tp)
    call HistoryPut('fluxconv',z_FluxConv_tp  )

    call HistoryClose
  end if
  
  ! textデータファイルの終了
  ! file close of output data in text format
  !
  if( OutText ) then
    close(fo)
  end if

end program NakajimaTest
