module globals !定数を設定
  real(8), parameter :: k = 3.64d-2 !吸収係数
  real(8), parameter  :: g = 9.8d0 !重力加速度
  real(8), parameter  :: sheight = 8.0d3 !大気のスケールハイト
  real(8), parameter  :: h = 6.63d-34 !プランク定数
  real(8), parameter  :: kappa_b = 1.38d-23 !ボルツマン定数
  real(8), parameter  :: sigma = 5.67d-8 !ステファン・ボルツマン定数
  real(8), parameter  :: c_p = 1.0d3 !大気の定圧比熱
  real(8), parameter  :: rho_s = 1.3d0 !大気地表面密度
  real(8), parameter  :: rho_abs = 1.0d-2 !水蒸気地表面密度
  real(8), parameter  :: dt = 2.88d4 !時間ステップの間隔
  integer, parameter :: sn = 500 !層の数
end module globals


module subprogs !モジュール副プログラムの開始
  use globals !定数を読み込む
  implicit none
contains

  function rho(z) result(r) !高度 z を入れると, その高度の密度 rho(z) を返す
    real(8), intent(in) :: z
    real(8) r
    r = rho_s * exp(-z / sheight)
  end function rho
  
  function z_t_s(z) result(tau_star) !zを入れるとスター付きの光学的厚さを返す(分割数10で数値積分)
    real(8) tau, dzprime, f, tau_star
    integer i
    real(8), intent(in) :: z
      tau = 0.0d0
      dzprime = ( 5.0d4 - z ) * 1.0d-1
      do i = 0, 10
        f = (z + dble(i) * dzprime) / sheight
        if(i == 0 .or. i == 10) then
          tau = tau + 0.5d0 * k * rho_abs * exp(-f) * dzprime
        else
          tau = tau + k * rho_abs * exp(-f) * dzprime
        endif
      enddo
      tau_star = 1.5d0 * tau
  end function z_t_s
  
end module subprogs

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!放射平衡の数値解を出力するためのプログラム
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
program main !主プログラム
use subprogs
  implicit none
  integer z , x, eighthours
  real(8) p0, height(sn), heightb(sn + 1), T(sn), T_a(sn), f_up(sn + 1), f_down(sn + 1), f(sn + 1),&
        & tau(sn + 1), dtau(sn), dtdt(sn), T_s, f_s
  p0 = rho_s * g *sheight !地表面気圧
  f_s = 0.7d0 * 1.37d3 * 0.25d0 !(1 - アルベド) * 太陽定数 / 4
  T(:) = 3.0d2 !初期温度を代入
  T_s = 3.0d2 !初期地表面温度を代入
  
  do x = 1, sn
    heightb(x) = -sheight * log(1.0d0 - ((dble(x) - 1.0d0)/(dble(sn))))
  enddo
  heightb(sn + 1) = 5.0d4
  
  
  do x = 1, sn + 1
    tau(x) = z_t_s(heightb(x))
  enddo
  do x = 1, sn
    height(x) =  0.5d0 * (heightb(x) + heightb(x + 1))
  enddo
  do x = 1, sn
    dtau(x) = tau(x) - tau(x + 1)
  enddo

    do eighthours = 1, 10000
      f_up(1) = sigma * (T_s ** 4)
    do x = 2, sn + 1 !上向きフラックスを計算
      f_up(x) = f_up(x - 1) * exp(-dtau(x - 1)) + sigma * (T(x - 1) ** 4) * exp(-dtau(x - 1)) * dtau(x - 1)
    enddo
      f_down(sn + 1) = 0.0d0
    do x = sn + 1, 2, -1!下向きフラックスを計算
      f_down(x - 1) = f_down(x) * exp(-dtau(x - 1)) + sigma * (T(x - 1) ** 4) * exp(-dtau(x - 1)) * dtau(x - 1) 
    enddo
    do x = 1, sn + 1!正味フラックスを計算
      f(x) = f_up(x) - f_down(x)
    enddo
    do x = 1, sn
      dtdt(x) = (f(x + 1) - f(x)) / ( -(heightb(x + 1) - heightb(x)) * c_p * rho(height(x)))
    enddo
    do x = 1, sn
      T_a(x) = T(x) + dtdt(x) * dt
    enddo
    T_s = sqrt(sqrt( (f_s + f_down(1)) / sigma)) ! 地表面温度
    do x = 1, sn
      T(x) = T_a(x)
    enddo
    
    if(abs(f(1) - f(sn + 1)) < 1.0d-2) exit !平衡条件をみたせば, ループから抜ける
  enddo

  write(*, *)'# 高度     , 光学的厚さ    , sigma * T(x)^4  , 温度 , sigma * T(x)^4の解析解, 温度の解析解'
  write(*, *) 0.0, z_t_s(0.0d0) * 2.0d0 / 3.0d0, sigma * (T_s ** 4), T_s, &
              & 0.5d0 * f_s *(z_t_s(0.0d0) + 2), sqrt(sqrt(0.5d0 * f_s *(z_t_s(0.0d0) + 2) / sigma))
  do x = 1, sn
    write(*, *) height(x), z_t_s(height(x)) * 2.0d0 / 3.0d0, sigma * (T(x) **4 ), T(x), &
              & 0.5d0 * f_s *(z_t_s(height(x)) + 1), sqrt(sqrt(0.5d0 * f_s *(z_t_s(height(x)) + 1) / sigma))
  enddo
  write(*, *) '#地表面温度', T_s
  !write(*, *) '#', f_up(sn + 1), f_s, eighthours, heightb(sn), height(sn)
end program main