!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
!
! Note that Japanese and English are described in parallel.
!
!
!== 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
!
! Note that Japanese and English are described in parallel.
!
! 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