
module globals
  implicit none
  real(8), parameter :: height = 5.0d0 * 10.0 ** 4.0              ! 大気の厚さ
  real(8), parameter :: scaleh = 8.0d0 * 10.0d0 ** 3.0            ! スケールハイト
  real(8), parameter :: rho0 = 1.23d0                              ! 地表面における大気の密度
  real(8), parameter :: arho0 = 1.0d0 * 10.0d0 ** (-2.0)          ! 地表面における吸収物質の密度
  real(8), parameter :: sigma = 5.67d0 * 10.0d0 ** (-8.0)         ! ステファンボルツマン定数
!  real(8), parameter :: k = 1.25d0 * 10.0d0 ** (-3.0)             ! 吸収係数1
  real(8), parameter :: k = 1.15 * 1.25d0 * 10.0d0 ** (-2.0)             ! 吸収係数2
!  real(8), parameter :: k = 3.75d0 * 10.0d0 ** (-2.0)             ! 吸収係数3
  real(8), parameter :: cp = 1.0d0 * 10.0d0 ** 3.0                ! 定圧比熱
  real(8), parameter :: gamma = 6.5d0 * 10.0d0 ** (-3.0)          ! 対流圏の気温減率
  real(8), parameter :: g = 9.8                                   ! 重力加速度
  real(8), parameter :: dt = 2.88d0 * 10.0d0 ** 4.0               ! 時間ステップ 8 時間
  real(8), parameter :: er0 = 1.0d-6                              ! 誤差
  integer, parameter :: n = 500                                   ! 分割した大気層の層数
  real(8), save :: dz                                             ! 大気層一枚当たりの厚さ
  real(8), save :: rho(0:n)                                       ! 大気の密度
  real(8), save :: arho(0:n)                                      ! 吸収物質の密度
  real(8), save :: tau(0:n)                                       ! 光学的厚さ
  real(8), save :: dtau(1:n)                                      ! 大気層一枚当たりの光学的厚さ
  real(8), save :: fu(0:n)                                        ! 上向き放射フラックス
  real(8), save :: fd(0:n)                                        ! 下向き放射フラックス
  real(8), save :: fn(0:n)                                        ! 正味の放射フラックス
  real(8) z(1:n)                                                  ! 高度
  real(8), save :: dp(1:n)                                        ! 気圧差
  real(8), save :: ts                                             ! 地表面温度
  real(8), save :: tsbk                                           ! 放射後の地表面温度のバックアップ
  real(8), save :: p(0:n)
  real(8), save :: pc(1:n)
end module globals

module value 
  use globals
  implicit none
contains
  subroutine valtau
    integer i
  
    dz = height / n

! 各層における大気と吸収物質の密度を計算
    
    do i = n, 0, -1
      arho(i) = arho0 * exp(-(n - i) * dz / scaleh )
      rho(i) = rho0 * exp(-(n - i) * dz / scaleh )
    enddo

    p(0) = 0.0d0
    do i = 1, n
       dp(i) = ((rho(i) + rho(i-1)) / 2.0) * g * dz
       p(i) = (p(i-1) + dp(i)) 
    enddo

    p(:) = p(:) * 0.01

    do i = 1, n
       pc(i) = 0.5 * (p(i-1) + p(i))
    enddo

! 一層当たりの光学的厚さと各層における光学的厚さを計算
    
    tau(0) = 0.0d0
    do i = 1, n
      dtau(i) = k * ((arho(i) + arho(i - 1))/2) * dz
!      dtau(i) = k * dp(i) / g
      tau(i) = tau(i-1) + dtau(i)
    enddo
    
  end subroutine valtau
end module value

module flux
  use globals
  use value
  implicit none
contains
  subroutine netflux(t)
    integer i
    real(8) fd0(0:n)                                  ! ある値
    real(8) fu2(0:n)                                  ! 上向き放射フラックスの第 2 項
    real(8) t(0:n)                                    ! 各層における温度
    
    call valtau

! 下向き放射フラックスを計算
    
    fd0(0) = 0.0d0                                    ! 大気上端における下向き放射フラックスはない
    do i = 1, n
      fd0(i) = fd0(i-1) + sigma * t(i) ** 4.0 * exp(1.5 * tau(i)) * dtau(i)
    enddo
    
    do i = 0, n
      fd(i) = 1.5 * fd0(i) * exp(-1.5 * tau(i)) 
    enddo

! 放射による地表面温度を計算
    
    ts = sqrt(sqrt((fd(n) + sigma * t(0) ** 4.0) / sigma ))
    tsbk = ts
    
! 上向き放射フラックスを計算
    
    fu2(n) = 0.0d0
    do i = n-1 , 0 , -1
       fu2(i) = fu2(i+1) + sigma * t(i+1) ** 4.0 * exp(-1.5 * tau(i)) * dtau(i+1)
    enddo

    do i = 0, n
       fu(i) = sigma * ts ** 4.0 * exp(-1.5 * (tau(n) - tau(i))) + 1.5 * exp(1.5 * tau(i)) * fu2(i)
    enddo
    
! 正味の放射フラックスを計算
    
    do i = 0, n
       fn(i) = fu(i) - fd(i)
    enddo

  end subroutine netflux

end module flux

program main
  use globals
  use value
  use flux
  implicit none

  real(8) t(0:n), tb(0:n)                                 ! 現時間とひとつ前の時間の, 大気上端の輝度温度と各層における温度(放射のみ)
  real(8) tc0(1:n)                                        ! 地表面温度が対流により変化したことによる影響を考慮した各層における温度
  real(8) tc(1:n), tac(1:n)                               ! n 回対流したことによる各層の温度, n+1 回対流したことによる各層の温度
  real(8) tday(1:n)                                       ! 1 day 前の温度
  real(8) ti, time, time1, tsc                            ! 初期温度, 経過時間(s), 経過時間(day)対流後の地表面温度
  real(8) const(1:n+1)                                    ! 係数 1
  real(8) const1(1:n)                                     ! 係数 2
  real(8) dfn(1:n)                                        ! 正味の放射フラックスの変化
 
  integer i, j, m
  integer day

  character(50) :: timename                               ! ファイル名

  ! 初期値を設定
  
!  ti = 1.7 * 10.0 ** 2.0                                  ! 170 K の等温大気
  ti = 3.0 * 10.0 ** 2.0                                  ! 300 K の等温大気  
  
  time = 0.0d0
  do j = 1, n
    tb(j) = ti
  enddo
  tb(0) = 2.55 * 10.0 ** 2.0                              ! 太陽放射の輝度温度
  t(0) = tb(0)

! 各層の高度を計算
  
  call valtau
  do i = n, 1, -1
    z(i) = ( n - i ) * dz + 0.5 * dz
 enddo

 open(11, file='file000000')
 do i = n, 1, -1
    write(11, *) tb(i), z(i)
 enddo
 close(11)

! open(12, file='ts.d')
 open(12, file='tsc.d')
 
1 continue

!  write(*,*) ts

  fn(:) = 0.0
  dfn(:) = 0.0

  call netflux(tb(0:n))
    
! 正味の放射フラックスの差を計算
  
  do i = 1, n
    dfn(i) = fn(i-1) - fn(i)
  enddo

! 各層における係数を計算
 
  do i = 1, n
    const(i) = dt / ( cp * dp(i) / g ) 
  enddo

! 次の時間ステップの各層における温度を計算(放射のみ)
    
  do i = 1, n
     t(i) = tb(i) - const(i) * dfn(i)
  enddo

  if(.true.) then

! 次の時間ステップの各層における温度を計算(放射 + 対流)     
     
3 continue

! 係数を計算
     
     do m = n, 1, -1

        const1(m) = 4.0 * sigma * g * dt / (cp * dp(m))

     enddo

! 地表面温度と接地層の不安定を解消する     

     if ( ts - t(n) .gt. gamma * 0.5 * dz) then

        tsc = ( t(n) + 0.5 * gamma * dz + const1(n) * exp(-1.5 * (tau(n) - tau(n-1))) * ts * tsbk ** 3.0 )&
             &/( 1.0 + const1(n) * exp(-1.5 * (tau(n) - tau(n-1))) * tsbk ** 3.0 ) 

        tc(n) = tsc - gamma * 0.5 * dz
        
!        write(*,*) tc(n) * cp * dp(n) / g , &
!                  & cp * t(n) * dp(n) / g - 4.0 * ts ** 3.0 * sigma * (tsc - ts) * exp(-(tau(n) - tau(n-1))) * dt

! 地表面温度が変化したことによる上向きフラックスの変化を補正する        
        
        do m = n-1, 1, -1
     
           tc0(m) = t(m)&
                &+ const1(m) * (tsc - ts) * tsbk ** 3.0 * ( exp(- 1.5 * (tau(n) - tau(m))) - exp(- 1.5 * (tau(n) - tau(m-1))) )

!           write(*,*) tc0(m) * cp * dp(m) / g , &
!                & cp * t(m) * dp(m) / g - 4.0 * ts ** 3.0 * sigma * (tsc - ts) * exp(-(tau(n) - tau(m-1))) * dt &
!                &+ 4.0 * ts ** 3.0 * sigma * (tsc - ts) * exp(-(tau(n) - tau(m))) * dt 
          
        enddo

     else

        tsc = ts
        tc(n) = t(n)
        do m = n-1, 1, -1
           tc0(m) = t(m)
        enddo
     
     endif

     do m = n-1, n-2, -1

     tc(m) = ((tc(m+1) - gamma * dz) * dp(m+1) + tc0(m) * dp(m))/ (dp(m+1) + dp(m))
                
     tac(m+1) = tc(m) + gamma * dz

  enddo


! 各層間の不安定を調べる
     
     do m = n-1, 1, -1

        if ( tc(m+1) - tc0(m) .gt. gamma * dz ) then

! 2 層を混ぜる              

           tc(m) = ((tc(m+1) - gamma * dz) * dp(m+1) + tc0(m) * dp(m))/ (dp(m+1) + dp(m))
     
           tac(m+1) = tc(m) + gamma * dz

!            write(*,*) m, cp * tac(m+1) * dp(m+1) / g + cp * tc(m) * dp(m) /g,&
!                 & cp * tc(m+1) * dp(m+1) / g + cp * tc0(m) * dp(m) / g
        
             
        else

           tc(m) = tc0(m)
           tac(m+1) = tc(m+1)
           
        endif
     enddo
                                                                                                                      
     tac(1) = tc(1)
        
! 再度不安定を調べる
        
     do m = n, 2, -1
   
        if (tac(m) - tac(m-1) .gt. gamma * dz + er0 ) then
     
           t(1:n) = tac(1:n)
           ts = tsc 
!           write(*,*) m
! 不安定がなくなるまでループさせる
              
           goto 3    

        endif
     enddo

endif


! 経過時間を計算

  time = time + dt
  time1 = time / (8.64d0 * 10.0 ** 4.0)

! 平衡条件を満たしているかの判断

  if( time1 .eq. real(time1)) then
!     write(12, *) time1, ts
     write(12, *) time1, tsc
     
     do m = n, 1, -1
        if (abs(tac(m) - tday(m)) .gt. 1.0d-4) then
!        if (abs(t(m) - tday(m)) .gt. 1.0d-4) then
!           tday(1:n) = t(1:n)
           tday(1:n) = tac(1:n)           
           goto 4
        endif
     enddo

!     do m = n, 1, -1
!        if(abs(tac(m) - tb(m)) .gt. er0) then
!     if(abs(t(m) - tb(m)) .gt. er0) then   
!           goto 4
!        endif
!     enddo
     write(*,*) sigma * t(0) ** 4.0&
          &- (fu(0) + 4.0 * sigma * tsbk ** 3.0 * (tsc-tsbk) * exp(-1.5 * (tau(n)-tau(0))))
    
  write(*,*) 'seconds = ', time
  write(*,*) 'days = ', time1
!  write(*,*) 'ts = ', ts
  write(*,*) 'tsc = ', tsc
  open(10, file = 'conv.d')
! open(10, file = 'radiation.d')
  write(10,*) tsc, 0.0
!  write(10,*) tsc, p(n)
!  write(10,*) ts, 0.0
  do i = n, 1, -1
!     write(10,*) tac(i), pc(i)
     write(10,*) tac(i), z(i)
!    write(10,*) t(i), z(i)
  enddo
  close(10)
  close(12)

  goto 2

endif
  
4 continue  
  
! ある時間における各層の温度を出力

  if ( mod(time1,50.0d0) .eq. 0) then
     write(*,*) time1
     day = real(time1)
    write(timename, '(a,i6.6)') 'file',day
    open(10, file=timename)
    write(10, *) tsc, 0.0
!    write(10, *) ts, 0.0
    do i = n, 1, -1
       write(10, *) tac(i), z(i)
!       write(10, *) t(i), z(i)
    enddo
    close(10)
  endif
  
! ループの上限

  if(time1 .gt. 3000) then

     write(*,*) sigma * t(0) ** 4.0&
          &- (fu(0) + 4.0 * sigma * tsbk ** 3.0 * (tsc-tsbk) * exp(-1.5 * (tau(n)-tau(0))))

     stop 'over time'

  endif
  
! 値の入れ替え
 
  do i = 1, n

     tb(i) = tac(i)
!     tb(i) = t(i)
     
  enddo

  write(*,*) time1
  
  goto 1
  
2 continue
  
end program main
