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  :: c = 3.00d8 !光速 !
  real(8), parameter  :: rho_s = 1.3d0 !大気地表面密度
  real(8), parameter  :: rho_abs = 1.0d-2 !水蒸気地表面密度
  real(8), parameter  :: dt = 2.88d4 !時間ステップの間隔
  integer, parameter :: sn = 100 !層の数
end module globals


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

  function rho(z) result(r)
    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

  write(*, *) '#高度, 光学的厚さ, sigma * T ** 4, 温度'
  write(*, *) 0.0, z_t_s(0.0d0)* 2.0d0 / 3.0d0, 0.5d0 * f_s *(z_t_s(dble(0)) + 2.0d0), &
              &sqrt(sqrt(0.5d0 * f_s *(z_t_s(dble(0)) + 2.0d0) / sigma))
  do z = 100, 50000, 100
  write(*, *) z, z_t_s(dble(z)) * 2.0d0 / 3.0d0, 0.5d0 * f_s *(z_t_s(dble(z)) + 1.0d0), &
              &sqrt(sqrt(0.5d0 * f_s *(z_t_s(dble(z)) + 1.0d0) / sigma))
  enddo
end program main