!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