| Class | ECCM |
| In: |
moist/eccm.f90
|
| Subroutine : | |||
| a_MolFrIni(1:SpcNum) : | real(8), intent(in)
| ||
| Humidity : | real(8), intent(in)
| ||
| z_Temp(DimZMin:DimZMax) : | real(8), intent(out)
| ||
| z_Press(DimZMin:DimZMax) : | real(8), intent(out)
| ||
| za_MolFr(DimZMin:DimZMax, 1:SpcNum) : | real(8), intent(out)
|
* 乾燥断熱温度減率に沿った温度、圧力を求め、それらに対して指定された相対湿度となるように凝結成分のモル比を決める * 比熱は乾燥気塊のもので代表させる * 流体の方程式において比熱は乾燥成分のもので代表させているため * 大気の平均分子量には湿潤成分の分子量を効果 * 流体の方程式において, 湿潤成分の分子量は考慮しているため
subroutine ECCM_Dry( a_MolFrIni, Humidity, z_Temp, z_Press, za_MolFr )
!
!== 概要
! * 乾燥断熱温度減率に沿った温度、圧力を求め、それらに対して指定された相対湿度となるように凝結成分のモル比を決める
! * 比熱は乾燥気塊のもので代表させる
! * 流体の方程式において比熱は乾燥成分のもので代表させているため
! * 大気の平均分子量には湿潤成分の分子量を効果
! * 流体の方程式において, 湿潤成分の分子量は考慮しているため
!
!暗黙の型宣言禁止
implicit none
real(8), intent(in) :: a_MolFrIni(1:SpcNum) !下部境界でのモル比
real(8), intent(in) :: Humidity !相対湿度 ( Humidity <= 1.0 )
real(8), intent(out):: z_Temp(DimZMin:DimZMax) !温度
real(8), intent(out):: z_Press(DimZMin:DimZMax)!圧力
! real(8), intent(out):: z_MolWtMean(DimZMin:DimZMax)
real(8) :: z_MolWtMean(DimZMin:DimZMax)
!平均分子量
real(8), intent(out):: za_MolFr(DimZMin:DimZMax, 1:SpcNum)
!モル分率
real(8) :: MolWtMeanDry !乾燥気塊の分子量の作業変数
real(8) :: MolWtMeanWet !湿潤気塊の分子量の作業変数
real(8) :: SatPress !飽和蒸気圧
real(8) :: VapPress !蒸気圧
real(8) :: DelMolFr
integer :: k, s
!-------------------------------------------------------------
! 配列の初期化
!-------------------------------------------------------------
!地表面の分子量を決める
za_MolFr(RegZMin+1, 1:SpcNum) = a_MolFrIni(1:SpcNum)
!地表面での平均分子量を決める
MolWtMeanDry = MolWtDry * (1.0d0 - sum(a_MolFrIni))
MolWtMeanWet = dot_product(MolWtWet, a_MolFrIni)
z_MolWtMean(RegZMin+1) = MolWtMeanDry + MolWtMeanWet
!地表面での温度(RegZMin+1 は, 高度 DelZ / 2 に相当)
z_Temp = 1.0d-60
z_Temp(RegZMin+1) = TempSfc - Grav * MolWtDry / CpDryMol * ( DelZ * 5.0d-1 )
!地表面での圧力(RegZMin+1 は, 高度 DelZ / 2 に相当)
z_Press = 1.0d-60
z_Press(RegZMin+1) = PressSfc *((TempSfc / z_Temp(RegZMin+1)) ** (- CpDryMol / GasRUniv))
!-----------------------------------------------------------
! (1) 乾燥断熱線に沿った温度を決める
! (2) 静水圧平衡から圧力を求める
! (3) (1),(2) の温度圧力に対して, とある相対湿度となるモル比を決める
!-----------------------------------------------------------
DtDz: do k = RegZMin+1, DimZMax-1
!(1)乾燥断熱線に沿って k+1 での温度を計算
z_Temp(k+1) = z_Temp(k) - Grav * MolWtDry / CpDryMol * DelZ
!念為
if (z_Temp(k+1) <= 0.0d0 ) z_Temp(k+1) = z_Temp(k)
!(2)圧力を静水圧平衡から計算
z_Press(k+1) = z_Press(k) * ((z_Temp(k) / z_Temp(k+1)) ** (- CpDryMol / GasRUniv))
!(3)モル比の計算
! まずはモル比は変化しないものとしてモル比を与える
! 飽和蒸気圧と平衡定数との平衡条件の前に適用しておく
za_MolFr(k+1,:) = za_MolFr(k,:)
do s = 1, CondNum
!飽和蒸気圧
SatPress = SvapPress( SpcWetID(IdxCC(s)), z_Temp(k+1) )
!元々のモル分率を用いて現在の蒸気圧を計算
VapPress = za_MolFr(k,IdxCG(s)) * z_Press(k+1)
!飽和蒸気圧と圧力から現在のモル比を計算
if ( VapPress > SatPress ) then
za_MolFr(k+1,IdxCG(s)) = max(SatPress * Humidity / z_Press(k+1), 1.0d-16)
end if
! write(*,*) k, s, z_Temp(k), za_MolFr(k,IdxCG(s)), VapPress, SatPress
end do
!NH4SH の平衡条件
if ( IdxNH3 /= 0 ) then
DelMolFr = max ( DelMolFrNH4SH( z_Temp(k+1), z_Press(k+1), za_MolFr(k+1,IdxNH3), za_MolFr(k+1,IdxH2S), Humidity ), 0.0d0 )
za_MolFr(k+1,IdxNH3) = za_MolFr(k+1,IdxNH3) - DelMolFr
za_MolFr(k+1,IdxH2S) = za_MolFr(k+1,IdxH2S) - DelMolFr
end if
!------------------------------------------------------------
!温度勾配を計算
!------------------------------------------------------------
!地表面での平均分子量を決める
MolWtMeanDry = MolWtDry * (1.0d0 - sum(za_MolFr(k+1, 1:SpcNum)))
MolWtMeanWet = dot_product(MolWtWet(1:SpcNum), za_MolFr(k+1, 1:SpcNum))
z_MolWtMean(k+1) = MolWtMeanDry + MolWtMeanWet
end do DtDz
end subroutine ECCM_Dry
| Subroutine : | |
| a_MolFrIni(1:SpcNum) : | real(8), intent(in) |
| Humidity : | real(8), intent(in) |
| z_Temp(DimZMin:DimZMax) : | real(8), intent(in) |
| z_Press(DimZMin:DimZMax) : | real(8), intent(in) |
| za_MolFr(DimZMin:DimZMax, 1:SpcNum) : | real(8), intent(out) |
与えられた温度に対し, 気塊が断熱的に上昇した時に実現される モル比のプロファイルを求める
subroutine ECCM_MolFr( a_MolFrIni, Humidity, z_Temp, z_Press, za_MolFr )
!
! 与えられた温度に対し, 気塊が断熱的に上昇した時に実現される
! モル比のプロファイルを求める
!
!暗黙の型宣言禁止
implicit none
real(8), intent(in) :: a_MolFrIni(1:SpcNum)
real(8), intent(in) :: Humidity
real(8), intent(in) :: z_Temp(DimZMin:DimZMax)
real(8), intent(in) :: z_Press(DimZMin:DimZMax)
real(8), intent(out):: za_MolFr(DimZMin:DimZMax, 1:SpcNum)
real(8) :: DelMolFr
integer :: k, s
!-----------------------------------------------------------
! 配列の初期化
!-----------------------------------------------------------
do s = 1, SpcNum
za_MolFr(:,s) = a_MolFrIni(s)
end do
!-----------------------------------------------------------
! 断熱減率 dT/dz の計算.
!-----------------------------------------------------------
do k = RegZMin, DimZMax
za_MolFr(k,:) = za_MolFr(k-1,:)
!------------------------------------------------------------
!NH4SH 以外の化学種の平衡条件
!------------------------------------------------------------
do s = 1, CondNum
!モル比を求める
!モル比は前のステップでのモル比を超えることはない
za_MolFr(k,IdxCG(s)) = min( za_MolFr(k-1,IdxCG(s)), SvapPress( SpcWetID(IdxCC(s)), z_Temp(k) ) * Humidity / z_Press(k) )
end do
!------------------------------------------------------------
!NH4SH の平衡条件
!------------------------------------------------------------
if ( IdxNH3 /= 0 ) then
!モル比の変化.
!とりあえず NH4SH に対する飽和比は 1.0 とする(手抜き...).
DelMolFr = max ( DelMolFrNH4SH( z_Temp(k), z_Press(k), za_MolFr(k,IdxNH3), za_MolFr(k,IdxH2S), Humidity ), 0.0d0 )
za_MolFr(k,IdxNH3) = za_MolFr(k,IdxNH3) - DelMolFr
za_MolFr(k,IdxH2S) = za_MolFr(k,IdxH2S) - DelMolFr
end if
end do
end subroutine ECCM_MolFr
| Subroutine : | |||
| z_Temp(DimZMin:DimZMax) : | real(8), intent(out)
| ||
| z_Press(DimZMin:DimZMax) : | real(8), intent(out)
| ||
| za_MolFr(DimZMin:DimZMax, 1:SpcNum) : | real(8), intent(out)
|
* Yamasaki, 1983 の温度と相対湿度の観測値で基本場を作成する
* 温度の基本場
* 観測データをNetCDFファイル化したものから値を読み込む
* 読み込んだ値を線形補間して温度の基本場を作成する
* 湿度の基本場
* 相対湿度は subroutine HUM で作成済み
* ここでは相対湿度をモル比に変換して湿度の基本場を作成する
* 気圧の基本場
* 湿潤成分と乾燥成分の分子量差を考慮した静水圧平衡から計算する
subroutine ECCM_N1994( z_Temp, z_Press, za_MolFr )
!
!== 概要
! * Yamasaki, 1983 の温度と相対湿度の観測値で基本場を作成する
! * 温度の基本場
! * 観測データをNetCDFファイル化したものから値を読み込む
! * 読み込んだ値を線形補間して温度の基本場を作成する
! * 湿度の基本場
! * 相対湿度は subroutine HUM で作成済み
! * ここでは相対湿度をモル比に変換して湿度の基本場を作成する
! * 気圧の基本場
! * 湿潤成分と乾燥成分の分子量差を考慮した静水圧平衡から計算する
!
implicit none
real(8), intent(out):: z_Press(DimZMin:DimZMax)!圧力
real(8), intent(out):: z_Temp(DimZMin:DimZMax) !温度
real(8), intent(out):: za_MolFr(DimZMin:DimZMax, 1:SpcNum) !モル比
real(8) :: Ob_altitude(1:89) !観測データに対応する高度
real(8) :: Ob_HumZ(1:88) !相対湿度の観測値(1次元)
real(8) :: z_HumZ(DimZMin:DimZMax) !相対湿度の基本場の値(%)
real(8) :: Observed_Temp(1:100,1:89)!NetCDF から読み込んだ気温の観測値(2次元)
real(8) :: Ob_TempZ(1:89) !NetCDF から読み込んだ気温の観測値(1次元)
real(8) :: z_PressSaturated(DimZMin:DimZMax)!Tetens(1930) の飽和蒸気圧
integer :: k, s, i, m, b
z_HumZ = 0.0d0
! subroutine HUMから相対湿度の基本場を取得
call HUM(z_HumZ, Ob_altitude, Ob_HumZ)
! 温度データをNetCDFから取得し線形補間する
! NetCDF ファイルからの気温の観測値の読み込み
call HistoryGet('Ob_Temp.nc', 'Ob_Temp', array=Observed_Temp)
! 温度の観測値のファイル"Ob_Temp.nc"は2次元データ(x,z)なので
! 1次元データ(z)にする
do m=1,89
Ob_TempZ(m) = Observed_Temp(1,m)
end do
! Ob_TempZ を線形補間して温度の基本場 z_Temp を作成する
! kansoku_altitude(m) < s_Z(k) < kansoku_altitude(m+1),
! kansoku_altitude = s_Z(k),
! kansoku_altitude < s_Z(k)
! の3つで場合分け
! 初期化
z_Temp = 1.0d-60
do k = RegZMin+1,DimZMax
do m = 1,88
! s_Z が Ob_altitude のある区間に挟まれるとき,
! その区間の Obaltitude(m), Obaltitude(m+1) を結ぶ
! 直線で Ob_TempZ を線形補間する
if (Ob_altitude(m).ne.s_Z(k).and.s_Z(k).gt.Ob_altitude(m) .and.s_Z(k).lt.Ob_altitude(m+1)) then
z_Temp(k) = Ob_TempZ(m)+((Ob_TempZ(m+1)-Ob_TempZ(m)) /(Ob_altitude(m+1)-Ob_altitude(m)))*(s_Z(k)-Ob_altitude(m))
! s_Z(k) = Obaltitude(m) の時は, Ob_TempZ(m) = z_Temp(k) である
else if (Ob_altitude(m).eq.s_Z(k)) then
z_Temp(k) = Ob_TempZ(m)
! s_Z(k) > Ob_altitude では観測データが無いので等温大気にする
else if (Ob_altitude(m).lt.s_Z(k)) then
z_Temp(k) = z_Temp(k-1)
end if
end do
end do
! 気圧とモル比を計算する
! 初期化
z_Press = 1.0d-60
! 地表面気圧と温度から大気最下層の気圧を求める
! 静水圧の式 dP/dz = - \rho * g を使用
! 簡単のため, \rho = 1 kg/m^3 の乾燥大気を仮定している
z_Press(RegZMin+1) = PressSfc - (Grav * DelZ)
! 大気最下層のモル比を計算する
za_MolFr(RegZMin+1,1) = SvapPress( 6,z_Temp(RegZMin+1)) * (z_HumZ(RegZMin+1)/100)/z_Press(RegZMin+1)
! Tetens(1930) の飽和蒸気圧も使ってみたいので以下で計算
do k=RegZMin+1,DimZMax
z_PressSaturated(k) = 6.11*(10**(7.5*(z_Temp(k)-273)/(z_Temp(k)-35.7)))*100
end do
! 大気最下層の値を元に, 気圧と湿度の基本場を求める
! 気圧は湿潤成分と乾燥成分の分子量差を考慮した静水圧平衡から計算する
! 相対湿度はモル比に変換する
! 以下のif文のDryHeightは基本場に乾燥した層を入れる場合に必要な値
! NAMELISTで鉛直領域以上の値を指定しておけば乾燥した層は入らない
! 湿度の基本場の計算
do k = RegZMin+1,DimZMax-1
if (s_Z(k+1).lt.DryHeight) then
za_MolFr(k+1,1) = SvapPress( 6,z_Temp(k)) * (z_HumZ(k)/100)/z_Press(k)
! 気圧の基本場の計算
z_Press(k+1) = z_Press(k)-(Grav*z_Press(k)*DelZ) /((GasRUniv/ (MolWtDry + za_MolFr(k,1) * (MolWtWet(1) - MolWtDry))) *z_Temp(k))
!DryHeightより上に領域があるときは上で求めたモル比をゼロにして
!乾燥大気における静水圧の式から気圧を求める
else if (s_Z(k+1).ge.DryHeight) then
za_MolFr(k+1,1) = 0.0d0
z_Press(k+1) = z_Press(k)-(Grav*z_Press(k)*DelZ) /((GasRUniv/MolWtDry)*z_Temp(k))
end if
end do
end subroutine ECCM_N1994
| Subroutine : | |||||
| xz_PotTemp(DimXMin:DimXMax,DimZMin:DimZMax) : | real(8), intent(in) | ||||
| xz_Exner(DimXMin:DimXMax, DimZMin:DimZMax) : | real(8), intent(in) | ||||
| xza_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) : | real(8), intent(in)
& xz_Stab, xz_StabTemp, xz_StabMolWt ) |
subroutine ECCM_Stab( xz_PotTemp, xz_Exner, xza_MixRt )
! & xz_Stab, xz_StabTemp, xz_StabMolWt )
use gridset, only: DimXMin, DimXMax, DimZMin, DimZMax, SpcNum !
use basicset,only: MolWtDry, MolWtWet, CpDry, Grav, xz_ExnerBasicZ, xz_PotTempBasicZ, xz_EffMolWtBasicZ, xza_MixRtBasicZ
use average, only: xz_avr_xr
use differentiate_center2, only: xr_dz_xz
implicit none
real(8), intent(in) :: xz_PotTemp(DimXMin:DimXMax,DimZMin:DimZMax)
real(8), intent(in) :: xz_Exner(DimXMin:DimXMax, DimZMin:DimZMax)
real(8), intent(in) :: xza_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
! real(8), intent(out) :: xz_Stab(DimXMin:DimXMax, DimZMin:DimZMax)
! real(8), intent(out) :: xz_StabTemp(DimXMin:DimXMax, DimZMin:DimZMax)
! real(8), intent(out) :: xz_StabMolWt(DimXMin:DimXMax, DimZMin:DimZMax)
real(8) :: xz_Stab(DimXMin:DimXMax, DimZMin:DimZMax)
real(8) :: xz_StabTemp(DimXMin:DimXMax, DimZMin:DimZMax)
real(8) :: xz_StabMolWt(DimXMin:DimXMax, DimZMin:DimZMax)
real(8) :: xza_MolFrAll(DimXMin:DimXMax,DimZMin:DimZMax,SpcNum)
real(8) :: xz_TempAll(DimXMin:DimXMax, DimZMin:DimZMax)
real(8) :: xz_MolWtWet(DimXMin:DimXMax, DimZMin:DimZMax)
integer :: i, k, s
xz_TempAll = (xz_PotTemp + xz_PotTempBasicZ) * (xz_Exner + xz_ExnerBasicZ)
do s = 1, SpcNum
xza_MolFrAll(:,:,s) = (xza_MixRt(:,:,s) + xza_MixRtBasicZ(:,:,s)) * MolWtDry / MolWtWet(s)
end do
do k = DimZMin, DimZMax
do i = DimXMin, DimXMax
xz_MolWtWet(i,k) = dot_product( MolWtWet(1:GasNum), xza_MolFrAll(i,k,1:GasNum) )
end do
end do
xz_StabTemp = Grav / xz_TempAll * ( xz_avr_xr( xr_dz_xz( xz_TempAll ) ) + Grav / CpDry )
xz_StabMolWt = - Grav * xz_avr_xr( xr_dz_xz( xz_MolWtWet ) ) / MolWtDry
xz_Stab = xz_StabTemp + xz_StabMolWt
call StoreStabTemp( xz_StabTemp )
call StoreStabMolWt( xz_StabMolWt )
where (xz_Stab < 1.0d-7)
xz_Stab = 1.0d-7
end where
end subroutine ECCM_Stab
| Subroutine : | |||
| z_Temp(DimZMin:DimZMax) : | real(8), intent(out)
| ||
| z_Press(DimZMin:DimZMax) : | real(8), intent(out)
| ||
| za_MolFr(DimZMin:DimZMax, 1:SpcNum) : | real(8), intent(out)
| ||
| z_LLhumidity(DimZMin:DimZMax) : | real(8), intent(out)
|
subroutine ECCM_Takemi2007( z_Temp, z_Press, za_MolFr, z_LLhumidity )
!
!== 概要
! * deepconv の地球用のテスト計算としてTakemi(2007)の再現計算を
! するための湿度の基本場を作成する
! * 基本場の温度の式が温位で与えられているため, 温度に変換する必要がある
!
implicit none
real(8), intent(out):: z_Press(DimZMin:DimZMax)!圧力
real(8), intent(out):: z_Temp(DimZMin:DimZMax) !温度
real(8), intent(out):: za_MolFr(DimZMin:DimZMax, 1:SpcNum) !モル比
real(8), intent(out) :: z_LLhumidity(DimZMin:DimZMax) ! 1.5 km までの相対湿度
!
! ECCM_Takemiでのみ必要な変数
!
real(8) :: z_TakemiRH(DimZMin:DimZMax) !高度 1.5km 以上の相対湿度
real(8) :: z_TakemiTheta(DimZMin:DimZMax) !基本場の温位
real(8) :: Theta_0 ! 地表面の温位
real(8) :: Theta_tr ! 対流圏界面の温位
real(8) :: TakemiZtr ! 対流圏界面の高度
real(8) :: HumSfc ! 地表面の相対湿度
real(8) :: a_MolFrSfc ! 地表面のモル比
!!--------------------------------------------------
integer :: k, s, i, m, b, T, L
z_TakemiRH = 0.0d0
z_TakemiTheta = 0.0d0
! 湿度場の選択
call HUM_Takemi2007(z_TakemiRH,TakemiZtr,T,L)
! call HUM_Takemi2007_TDRY(z_TakemiRH,TakemiZtr,T,L)
! call HUM_Takemi2007_MDRY(z_TakemiRH,TakemiZtr,T,L)
! 温度・湿度・気圧を求める
!地表面の温位は300K
Theta_0 = 300.0d0
! 中緯度場と熱帯場の選択
! 対流圏界面の温位で区別している
! Theta_tr = 343.0d0 ! 中緯度
Theta_tr = 358.0d0 ! 熱帯
! 大気最下層の温位
z_TakemiTheta(RegZMin+1) = Theta_0 + (Theta_tr - Theta_0)*(s_Z(RegZMin+1) / TakemiZtr)**1.25
z_Temp = 1.0d-60
z_Press = 1.0d-60
z_LLhumidity = 0.0d0
! 大気最下層の気温・気圧・水蒸気混合比・分子量の計算
! 地表面のモル比
! 10g/kg〜18g/kgの混合比で与えられているためモル比に直す
! なお熱帯場はQ18のみ使用
a_MolFrSfc = 0.0180d0 * MolWtDry / MolWtWet(1) ! Q18
! a_MolFrSfc = 0.0160d0 * MolWtDry / MolWtWet(1) ! Q16
! a_MolFrSfc = 0.0140d0 * MolWtDry / MolWtWet(1) ! Q14
! a_MolFrSfc = 0.0120d0 * MolWtDry / MolWtWet(1) ! Q12
! a_MolFrSfc = 0.0100d0 * MolWtDry / MolWtWet(1) ! Q10
! 地表面の相対湿度
! モル比から逆算. 気圧を求めるのに必要
HumSfc = a_MolFrSfc * PressSfc / SvapPress(6,TempSfc)
! 大気最下層の気圧
z_Press(RegZMin+1) = PressSfc - (Grav * PressSfc * DelZ) /((GasRUniv/ (MolWtDry + a_MolFrSfc * (MolWtWet(1) - MolWtDry))) *TempSfc)
! 大気最下層の温度. 温位を元に計算
! 地表面気圧は 1000hPa と仮定する
z_Temp(RegZMin+1) = z_TakemiTheta(RegZMin+1) * (z_Press(RegZMin+1) / 100000.0d0)**(GasRDry / CpDry)
! 大気最下層のモル比
! 熱帯場はQ18のみ使用
za_MolFr(RegZMin+1,1) = 0.0180d0 * MolWtDry / MolWtWet(1) ! Q18
! za_MolFr(RegZMin+1,1) = 0.0160d0 * MolWtDry / MolWtWet(1) ! Q16
! za_MolFr(RegZMin+1,1) = 0.0140d0 * MolWtDry / MolWtWet(1) ! Q14
! za_MolFr(RegZMin+1,1) = 0.0120d0 * MolWtDry / MolWtWet(1) ! Q12
! za_MolFr(RegZMin+1,1) = 0.0100d0 * MolWtDry / MolWtWet(1) ! Q10
! 大気最下層の相対湿度
! モル比から逆算. 気圧を求めるのに必要
z_LLhumidity(RegZMin+1) = za_MolFr(RegZMin+1,1) * z_Press(RegZMin+1) / SvapPress(6,z_Temp(RegZMin+1))
! 大気最下層の種々の物理量をもとに温度・湿度気圧の基本場を作成する
! 高度 0 - 1.5km と 1.5km - で湿度の与え方が異なり,
! さらに, 対流圏界面以上では温度の式が異なるので, 以下の 3パターンに分けて
! 基本場を求める
! 高度 0 - 1.5km, 1.5km - 対流圏界面(Ztr), Ztr以上
do k = RegZMin+1, DimZMax
!
! 基本場(I)
! 大気最下層から高度 1.5 km 以下(k が L-1 以下)の基本場
!
if (k.le.L-1) then
! モル比
! 高度 1.5km まではモル比が任意の値で固定されている
! ただし, 熱帯場は 18 g/kg のみ
za_MolFr(k+1,1) = 0.0180d0 * MolWtDry / MolWtWet(1)
! za_MolFr(k+1,1) = 0.0160d0 * MolWtDry / MolWtWet(1)
! za_MolFr(k+1,1) = 0.0140d0 * MolWtDry / MolWtWet(1)
! za_MolFr(k+1,1) = 0.0120d0 * MolWtDry / MolWtWet(1)
! za_MolFr(k+1,1) = 0.0100d0 * MolWtDry / MolWtWet(1)
! 気圧
z_Press(k+1) = z_Press(k)-(Grav*z_Press(k)*DelZ) /((GasRUniv/ (MolWtDry + za_MolFr(k,1) * (MolWtWet(1) - MolWtDry))) *z_Temp(k))
! 温位から温度を計算
z_TakemiTheta(k+1) = Theta_0 + (Theta_tr - Theta_0)*(s_Z(k+1) / TakemiZtr)**1.25
z_Temp(k+1) = z_TakemiTheta(k+1) * (z_Press(k+1) / 100000.0d0)**(GasRDry / CpDry)
! モル比から相対湿度を計算
z_LLhumidity(k+1) = za_MolFr(k+1,1) * z_Press(k+1) / SvapPress(6,z_Temp(k+1))
!
!== Takemi(2007)の基本場の使用上の注意
! * 低層で混合比を固定し続けると高度 1.5km までに相対湿度が100%を超える高度がある
! * 100% を超えてしまった高度の湿度は, その直下の高度の湿度にすると
! Takemi(2007)と同じ湿度の鉛直プロファイルが描ける
! * なぞ
! * MQ16 は 湿度 99.・・・% という値が出て, 100% を超えなくてもほぼ100%
! な値が出てしまう
! * Takemi(2007)ではこのような値があるようには見えない
! しょうがないのでMQ18 とMQ14 の 100% を超える直前の値の間を取った値を入れておく
! * 0.94511507230d0という怪しい数字が上記の値
if (z_LLhumidity(k+1).gt.1) then
z_LLhumidity(k+1) = z_LLhumidity(k)
! if (z_LLhumidity(k+1).gt.0.99) then !MQ16 用
! z_LLhumidity(k+1) = 0.94511507230d0 !MQ16 用
end if
!
! 基本場(II)
! 高度 1.5km より上で, 対流圏界面以下(L < k < T)の基本場
!
else if (k.le.T.and.k.gt.L) then
! k = L での相対湿度
! この if 文の初めの k は L+1 から始まり z_TakemiRH(L) が区間外なので値を与える
z_TakemiRH(L) = z_LLhumidity(L)
! 気圧の基本場
z_Press(k) = z_Press(k-1)-(Grav*z_Press(k-1)*DelZ) /((GasRUniv/ (MolWtDry + za_MolFr(k+1,1) * (MolWtWet(1) - MolWtDry))) *z_Temp(k-1))
! 温位から温度を計算
z_TakemiTheta(k) = Theta_0 + (Theta_tr - Theta_0)*(s_Z(k) / TakemiZtr)**1.25
z_Temp(k) = z_TakemiTheta(k) * (z_Press(k) / 100000.0d0)**(GasRDry / CpDry)
! モル比と大気の平均分子量
za_MolFr(k,1) = SvapPress( 6,z_Temp(k))* z_TakemiRH(k)/z_Press(k)
!
! 基本場(III)
! 対流圏界面より上 (k > T) の場
!
else if (k.gt.T) then
! 気圧の基本場
z_Press(k) = z_Press(k-1)-(Grav*z_Press(k-1)*DelZ) /((GasRUniv/ (MolWtDry + za_MolFr(k,1) * (MolWtWet(1) - MolWtDry))) *z_Temp(k-1))
! 温位から温度を計算
z_TakemiTheta(k) = Theta_tr * exp(Grav * (s_Z(k) - TakemiZtr) / (CpDry * z_Temp(T)))
z_Temp(k) = z_TakemiTheta(k) * (z_Press(k) / 100000.0d0)**(GasRDry / CpDry)
! モル比と大気の平均分子量
za_MolFr(k,1) = SvapPress( 6,z_Temp(k))* z_TakemiRH(k)/z_Press(k)
end if
end do
end subroutine ECCM_Takemi2007
| Subroutine : | |||
| a_MolFrIni(1:SpcNum) : | real(8), intent(in)
| ||
| Humidity : | real(8), intent(in)
| ||
| z_Temp(DimZMin:DimZMax) : | real(8), intent(out)
| ||
| z_Press(DimZMin:DimZMax) : | real(8), intent(out)
| ||
| za_MolFr(DimZMin:DimZMax, 1:SpcNum) : | real(8), intent(out)
|
* 湿潤断熱温度減率に沿った温度、圧力を求め、それらに対して指定された相対湿度となるように凝結成分のモル比を決める * 比熱は乾燥気塊のもので代表させる * 流体の方程式において比熱は乾燥成分のもので代表させているため * 大気の平均分子量には湿潤成分の分子量を効果 * 流体の方程式において, 湿潤成分の分子量は考慮しているため
subroutine ECCM_Wet( a_MolFrIni, Humidity, z_Temp, z_Press, za_MolFr )
!
!== 概要
! * 湿潤断熱温度減率に沿った温度、圧力を求め、それらに対して指定された相対湿度となるように凝結成分のモル比を決める
! * 比熱は乾燥気塊のもので代表させる
! * 流体の方程式において比熱は乾燥成分のもので代表させているため
! * 大気の平均分子量には湿潤成分の分子量を効果
! * 流体の方程式において, 湿潤成分の分子量は考慮しているため
!
!暗黙の型宣言禁止
implicit none
real(8), intent(in) :: a_MolFrIni(1:SpcNum) !下部境界でのモル比
real(8), intent(in) :: Humidity !相対湿度 ( Humidity <= 1.0 )
real(8), intent(out):: z_Temp(DimZMin:DimZMax) !温度
real(8), intent(out):: z_Press(DimZMin:DimZMax)!圧力
! real(8), intent(out):: z_MolWtMean(DimZMin:DimZMax)
real(8) :: z_MolWtMean(DimZMin:DimZMax)
!平均分子量
real(8), intent(out):: za_MolFr(DimZMin:DimZMax, 1:SpcNum)
!モル分率
real(8) :: MolWtMeanDry !乾燥気塊の分子量の作業変数
real(8) :: MolWtMeanWet !湿潤気塊の分子量の作業変数
real(8) :: SatPress !飽和蒸気圧
real(8) :: VapPress !蒸気圧
real(8) :: DelMolFr
integer :: k, s, j
real(8) :: TempA, PressA, a_MolFrA(SpcNum)
real(8) :: TempN, PressN, a_MolFrN(SpcNum)
real(8) :: DZ
real(8) :: A, B
real(8) :: Heat(SpcNum)
real(8) :: ReactHeat
logical :: MoistFlag(SpcNum), ReactFlag
!-------------------------------------------------------------
! 配列の初期化
!-------------------------------------------------------------
!地表面の分子量を決める
za_MolFr(RegZMin+1, 1:SpcNum) = a_MolFrIni(1:SpcNum)
!地表面での平均分子量を決める
MolWtMeanDry = MolWtDry * (1.0d0 - sum(a_MolFrIni))
MolWtMeanWet = dot_product(MolWtWet, a_MolFrIni)
z_MolWtMean(RegZMin+1) = MolWtMeanDry + MolWtMeanWet
!地表面での温度(RegZMin+1 は, 高度 DelZ / 2 に相当)
z_Temp = 1.0d-60
z_Temp(RegZMin+1) = TempSfc - Grav * MolWtDry / CpDryMol * ( DelZ * 5.0d-1 )
!地表面での圧力(RegZMin は, 高度 DelZ / 2 に相当)
z_Press = 1.0d-60
z_Press(RegZMin+1) = PressSfc *((TempSfc / z_Temp(RegZMin+1)) ** (- CpDryMol / GasRUniv))
!-----------------------------------------------------------
! 断熱減率 dT/dz の計算.
!-----------------------------------------------------------
DtDz: do k = RegZMin+1, DimZMax-1
TempN = z_Temp(k)
PressN = z_Press(k)
a_MolFrN = za_MolFr(k,:)
DZ = DelZ/100
do j = 1, 100
MoistFlag = .false.
Heat = 0.0d0
a_MolFrA = a_MolFrN !初期化
!乾燥断熱線に沿って k+1 での温度を計算
TempA = TempN - Grav * MolWtDry / CpDryMol * DZ
! write(*,*) "TempA", TempA
!念為
if (TempA <= 0.0d0 ) TempA = TempN
!圧力を静水圧平衡から計算
PressA = PressN * ((TempN / TempA) ** (- CpDryMol / GasRUniv))
! 凝結が生じるか判定
do s = 1, CondNum
!飽和蒸気圧
SatPress = SvapPress( SpcWetID(IdxCC(s)), TempA )
!元々のモル分率を用いて現在の蒸気圧を計算
VapPress = a_MolFrN(IdxCG(s)) * PressA
!飽和蒸気圧から凝結の有無を決める
if ( VapPress < SatPress ) then
!凝結していないので潜熱なし.
a_MolFrA(IdxCG(s)) = a_MolFrN(IdxCG(s))
!凝結しないので潜熱なし
Heat(IdxCG(s)) = 0.0d0
else
!潜熱. 現在の高度での値
Heat(IdxCG(s)) = LatentHeatPerMol( SpcWetID(IdxCC(s)), TempN )
!凝結のスイッチ
MoistFlag(s) = .true.
end if
end do
!NH4SH の平衡条件
if ( IdxNH3 /= 0 ) then
DelMolFr = max ( DelMolFrNH4SH( TempA, PressA, a_MolFrN(IdxNH3), a_MolFrN(IdxH2S), Humidity ), 0.0d0 )
! 反応熱
ReactHeat = ReactHeatNH4SHPerMol * DelMolFr
! 反応のスイッチ
if (DelMolFr > 0) ReactFlag = .true.
end if
! write(*,*) "MoistFlag", MoistFlag
! write(*,*) "ReactFlag", ReactFlag
if (count(MoistFlag) > 0 .or. ReactFlag) then
!係数. 高度 k での値で評価
A = dot_product( Heat(1:SpcNum), a_MolFrN(1:SpcNum)) / ( GasRUniv * TempN )
B = dot_product(( Heat(1:SpcNum) ** 2.0d0), a_MolFrN(1:SpcNum)) / ( CpDryMol * GasRUniv * ( TempN ** 2.0d0 ) )
!断熱温度減率
TempA = TempN - Grav * MolWtDry * (1.0d0 + A) / (CpDryMol * (1.0d0 + B)) * DZ + ReactHeat / (CpDryMol * DelZ)
! write(*,*) "Moist TempA", TempA
! write(*,*) "Moist Heat ", Heat
! write(*,*) "Moist MolFr ", a_MolFrN
! write(*,*) "Moist DTempDZ", - Grav * MolWtDry * (1.0d0 + A) / (CpDryMol * (1.0d0 + B))
!念為
if (TempA <= 0.0d0 ) TempA = TempN
!圧力を静水圧平衡から計算
PressA = PressN * ((TempN / TempA) ** (- CpDryMol / GasRUniv))
do s = 1, CondNum
if (MoistFlag(s)) then
!飽和蒸気圧
SatPress = SvapPress( SpcWetID(IdxCC(s)), TempA )
! 飽和モル比
a_MolFrA(IdxCG(s)) = max(SatPress / PressA, 1.0d-16)
else
! 前のモル比と同じ
a_MolFrA(IdxCG(s)) = a_MolFrN(IdxCG(s))
end if
end do
if ( ReactFlag ) then
DelMolFr = max ( DelMolFrNH4SH( TempA, PressA, a_MolFrA(IdxNH3), a_MolFrA(IdxH2S), Humidity ), 0.0d0 )
a_MolFrA(IdxNH3) = a_MolFrA(IdxNH3) - DelMolFr
a_MolFrA(IdxH2S) = a_MolFrA(IdxH2S) - DelMolFr
end if
end if
TempN = TempA
PressN = PressA
a_MolFrN = a_MolFrA
end do
z_Temp(k+1) = TempA
z_Press(k+1) = PressA
za_MolFr(k+1,:) = a_MolFrA(:)
!平均分子量を決める
MolWtMeanDry = MolWtDry * (1.0d0 - sum(za_MolFr(k+1, 1:SpcNum)))
MolWtMeanWet = dot_product(MolWtWet(1:SpcNum), za_MolFr(k+1, 1:SpcNum))
z_MolWtMean(k+1) = MolWtMeanDry + MolWtMeanWet
end do DtDz
end subroutine ECCM_Wet
| Subroutine : | |||
| Idx : | integer | ||
| Humidity : | real(8), intent(in)
| ||
| z_Temp(DimZMin:DimZMax) : | real(8), intent(inout)
| ||
| z_Press(DimZMin:DimZMax) : | real(8), intent(inout)
| ||
| za_MolFr(DimZMin:DimZMax, 1:SpcNum) : | real(8), intent(inout)
|
* 湿潤断熱温度減率に沿った温度、圧力を求め、それらに対して指定された相対湿度となるように凝結成分のモル比を決める * 比熱は乾燥気塊のもので代表させる * 流体の方程式において比熱は乾燥成分のもので代表させているため * 大気の平均分子量には湿潤成分の分子量を効果 * 流体の方程式において, 湿潤成分の分子量は考慮しているため
subroutine ECCM_Wet2( Idx, Humidity, z_Temp, z_Press, za_MolFr)
!
!== 概要
! * 湿潤断熱温度減率に沿った温度、圧力を求め、それらに対して指定された相対湿度となるように凝結成分のモル比を決める
! * 比熱は乾燥気塊のもので代表させる
! * 流体の方程式において比熱は乾燥成分のもので代表させているため
! * 大気の平均分子量には湿潤成分の分子量を効果
! * 流体の方程式において, 湿潤成分の分子量は考慮しているため
!
!暗黙の型宣言禁止
implicit none
integer :: Idx
real(8), intent(inout) :: z_Temp(DimZMin:DimZMax) !下部境界でのモル比
real(8), intent(inout) :: z_Press(DimZMin:DimZMax) !下部境界でのモル比
real(8), intent(inout) :: za_MolFr(DimZMin:DimZMax, 1:SpcNum) !下部境界でのモル比
real(8), intent(in) :: Humidity !相対湿度 ( Humidity <= 1.0 )
real(8) :: SatPress !飽和蒸気圧
real(8) :: VapPress !蒸気圧
real(8) :: DelMolFr
integer :: k, s, j
real(8) :: TempA, PressA, a_MolFrA(SpcNum)
real(8) :: TempN, PressN, a_MolFrN(SpcNum)
real(8) :: DZ
real(8) :: A, B
real(8) :: Heat(SpcNum)
real(8) :: ReactHeat
logical :: MoistFlag(SpcNum), ReactFlag
!-----------------------------------------------------------
! 断熱減率 dT/dz の計算.
!-----------------------------------------------------------
DtDz: do k = Idx, DimZMax-1
TempN = z_Temp(k)
PressN = z_Press(k)
a_MolFrN = za_MolFr(k,:)
DZ = DelZ/100
do j = 1, 100
MoistFlag = .false.
Heat = 0.0d0
a_MolFrA = a_MolFrN !初期化
!乾燥断熱線に沿って k+1 での温度を計算
TempA = TempN - Grav * MolWtDry / CpDryMol * DZ
! write(*,*) "TempA", TempA
!念為
if (TempA <= 0.0d0 ) TempA = TempN
!圧力を静水圧平衡から計算
PressA = PressN * ((TempN / TempA) ** (- CpDryMol / GasRUniv))
! 凝結が生じるか判定
do s = 1, CondNum
!飽和蒸気圧
SatPress = SvapPress( SpcWetID(IdxCC(s)), TempA )
!元々のモル分率を用いて現在の蒸気圧を計算
VapPress = a_MolFrN(IdxCG(s)) * PressA
!飽和蒸気圧から凝結の有無を決める
if ( VapPress < SatPress ) then
!凝結していないので潜熱なし.
a_MolFrA(IdxCG(s)) = a_MolFrN(IdxCG(s))
!凝結しないので潜熱なし
Heat(IdxCG(s)) = 0.0d0
else
!潜熱. 現在の高度での値
Heat(IdxCG(s)) = LatentHeatPerMol( SpcWetID(IdxCC(s)), TempN )
!凝結のスイッチ
MoistFlag(s) = .true.
end if
end do
!NH4SH の平衡条件
if ( IdxNH3 /= 0 ) then
DelMolFr = max ( DelMolFrNH4SH( TempA, PressA, a_MolFrN(IdxNH3), a_MolFrN(IdxH2S), Humidity ), 0.0d0 )
! 反応熱
ReactHeat = ReactHeatNH4SHPerMol * DelMolFr
! 反応のスイッチ
if (DelMolFr > 0) ReactFlag = .true.
end if
! write(*,*) "MoistFlag", MoistFlag
! write(*,*) "ReactFlag", ReactFlag
if (count(MoistFlag) > 0 .or. ReactFlag) then
!係数. 高度 k での値で評価
A = dot_product( Heat(1:SpcNum), a_MolFrN(1:SpcNum)) / ( GasRUniv * TempN )
B = dot_product(( Heat(1:SpcNum) ** 2.0d0), a_MolFrN(1:SpcNum)) / ( CpDryMol * GasRUniv * ( TempN ** 2.0d0 ) )
!断熱温度減率
TempA = TempN - Grav * MolWtDry * (1.0d0 + A) / (CpDryMol * (1.0d0 + B)) * DZ + ReactHeat / (CpDryMol * DZ)
! write(*,*) "Moist TempA", TempA
! write(*,*) "Moist Heat ", Heat
! write(*,*) "Moist MolFr ", a_MolFrN
! write(*,*) "Moist DTempDZ", - Grav * MolWtDry * (1.0d0 + A) / (CpDryMol * (1.0d0 + B))
!念為
if (TempA <= 0.0d0 ) TempA = TempN
!圧力を静水圧平衡から計算
PressA = PressN * ((TempN / TempA) ** (- CpDryMol / GasRUniv))
do s = 1, CondNum
if (MoistFlag(s)) then
!飽和蒸気圧
SatPress = SvapPress( SpcWetID(IdxCC(s)), TempA )
! 飽和モル比
a_MolFrA(IdxCG(s)) = max(SatPress / PressA, 1.0d-16)
else
! 前のモル比と同じ
a_MolFrA(IdxCG(s)) = a_MolFrN(IdxCG(s))
end if
end do
if ( ReactFlag ) then
DelMolFr = max ( DelMolFrNH4SH( TempA, PressA, a_MolFrA(IdxNH3), a_MolFrA(IdxH2S), Humidity ), 0.0d0 )
a_MolFrA(IdxNH3) = a_MolFrA(IdxNH3) - DelMolFr
a_MolFrA(IdxH2S) = a_MolFrA(IdxH2S) - DelMolFr
end if
end if
TempN = TempA
PressN = PressA
a_MolFrN = a_MolFrA
end do
z_Temp(k+1) = TempA
z_Press(k+1) = PressA
za_MolFr(k+1,:) = a_MolFrA(:)
end do DtDz
end subroutine ECCM_Wet2
| Subroutine : | |
| cfgfile : | character(*), intent(in) |
ECCM モジュールの初期化ルーチン
This procedure input/output NAMELIST#eccm .
subroutine ECCM_init(cfgfile)
!
! ECCM モジュールの初期化ルーチン
!
!暗黙の型宣言禁止
implicit none
!
!NAMELIST から放射強制の設定を取得
!
character(*), intent(in) :: cfgfile
! NAMELIST から情報を取得
NAMELIST /eccm/ DryHeight
open (10, FILE=cfgfile)
read(10, NML=eccm)
close(10)
!確認
call MessageNotify( "M", "ECCM_init", "DryHeight = %f", d=(/DryHeight/) )
end subroutine ECCM_init
| Subroutine : | |||
| z_humid(DimZMin:DimZMax) : | real(8),intent(out)
| ||
| kansoku_altitude(1:89) : | real(8),intent(out)
| ||
| kansoku_HumZ(1:88) : | real(8),intent(out)
|
* Yamasaki, 1983 の湿度の観測データから基本場を作成する * 相対湿度の観測データをNetCDFファイル化したものから値を読み込む * 読み込んだ値を線形補間して相対湿度の基本場を作成する
subroutine HUM(z_humid, kansoku_altitude, kansoku_HumZ)
!
!== 概要
!
! * Yamasaki, 1983 の湿度の観測データから基本場を作成する
! * 相対湿度の観測データをNetCDFファイル化したものから値を読み込む
! * 読み込んだ値を線形補間して相対湿度の基本場を作成する
!
implicit none
real(8) :: Observed_Humidity(1:100,1:88) !NetCDF から読み込んだ相対湿度の観測値
!水平一様な2次元データ
real(8),intent(out) :: kansoku_HumZ(1:88) !鉛直1次元にした相対湿度の観測値の代入先
real(8),intent(out) :: z_humid(DimZMin:DimZMax) !線形補間した相対湿度
real(8),intent(out) :: kansoku_altitude(1:89) !観測データに対応する高度
integer :: k,m
z_humid = 0.0d0
kansoku_altitude = 0.0d0
kansoku_HumZ = 0.0d0
! 観測データに対応する高度を与える
! これは線形補間の際に使用する
! なお, Yamsaki(1983)にある元データはデータ間隔がばらばらであったので
! 250m間隔に直したものを使用している
do m = 1,89
kansoku_altitude(m) = 25.0d1*(m-1)
end do
! NetCDF ファイルから相対湿度を読み込む
call HistoryGet('Ob_Humidity.nc', 'Ob_Humidity', array=Observed_Humidity)
! 相対湿度の観測値のファイル"Ob_Humidity.nc"は2次元データ(x,z)であるため
! 1次元データにする
do m=1,88
kansoku_HumZ(m) = Observed_Humidity(1,m)
end do
! kansoku_HumZ を線形補間してz_humidを作成する
! kansoku_altitude(m) < s_Z(k) < kansoku_altitude(m+1),
! kansoku_altitude = s_Z(k),
! kansoku_altitude < s_Z(k)
! の3つで場合分け
do k = RegZMin+1,DimZMax
do m = 1,87
if (kansoku_altitude(m).ne.s_Z(k).and.s_Z(k).gt.kansoku_altitude(m) .and.s_Z(k).lt.kansoku_altitude(m+1)) then
z_humid(k) = kansoku_HumZ(m)+((kansoku_HumZ(m+1) - kansoku_HumZ(m)) /(kansoku_altitude(m+1)-kansoku_altitude(m)))*(s_Z(k)-kansoku_altitude(m))
else if (kansoku_altitude(m).eq.s_Z(k)) then
z_humid(k) = kansoku_HumZ(m)
else if (kansoku_altitude(m).lt.s_Z(k)) then
z_humid(k) = z_humid(k-1)
end if
end do
end do
end subroutine HUM
| Subroutine : | |||
| z_RH(DimZMin:DimZMax) : | real(8),intent(out)
| ||
| Ztr : | real(8),intent(out)
| ||
| TR : | integer,intent(out)
| ||
| LL : | integer,intent(out)
|
subroutine HUM_Takemi2007(z_RH,Ztr,TR,LL)
!
!== 概要
! * deepconv の地球用のテスト計算としてTakemi(2007)の再現計算を
! するための湿度の基本場を作成する
! * Takemi(2007) の熱帯場と中緯度場の相対湿度の"BASE CASE"の
! 高度 1.5km 以上の湿度を計算している(詳細はTakemi(2007)の図2を参照のこと)
! * Takemi(2007)では高度1.5km以上はある関数で与えているが,
! 高度1.5km以下は一定の混合比で与えているため,
! 高度1.5km以下の湿度は温度や気圧が求まらないと得られないので後で求める
!
implicit none
real(8),intent(out) :: z_RH(DimZMin:DimZMax) !Relative Humidity(相対湿度)
real(8),intent(out) :: Ztr !対流圏界面の高度
integer,intent(out) :: TR !対流圏界面の高度の格子番号
integer,intent(out) :: LL !高度 1.5 km 直下の格子番号 (Low Level)
integer :: k
z_RH = 0.0d0
Ztr = 0.0d0
! Ztrに対流圏界面高度の値を与える
! Takemi(2007)では12kmだが, 必ずしも格子点がちょうど12kmには乗らないため,
! 12kmか,12kmの直上の位置を対流圏界面にする
do k = RegZMin+1,DimZMax
! s_Z が ちょうど 12 km なら 12 km が対流圏界面
if (s_Z(k).eq.12000) then
Ztr = s_Z(k)
! s_Z がちょうど12kmではなかった場合,12kmの直上のs_Zが対流圏界面Ztrになる
else if (s_Z(k).gt.12000) then
Ztr = minval( s_Z, 1, s_Z > 12000 )
end if
end do
! 対流圏界面高度の格子点を求める (Ztr = s_Z(TR))
do k = RegZMin+1,DimZMax
if (s_Z(k).ge.12000.and.s_Z(k).le.(12000 + DelZ)) then
TR = k
end if
end do
! 高度1.5kmもしくはその直下の格子番号を求める
do k = RegZMin+1,DimZMax
if (s_Z(k).le.1500.and.s_Z(k).ge.(1500 - DelZ)) then
LL = k
end if
end do
! Takemi(2007)の式を用いて相対湿度を求める
! とりあえず高度1.5km以下の湿度はゼロにしておく
do k = RegZMin+1, DimZMax
if (k.le.LL) then
z_RH(k) = 0.0d0
else if (k.le.TR.and.k.gt.LL) then
z_RH(k) = 1.0d0 - 0.750d0 * (s_Z(k) / Ztr)**1.25
else if (k.gt.TR) then
z_RH(k) = 0.25
end if
end do
end subroutine HUM_Takemi2007
| Subroutine : | |
| z_RHall(DimZMin:DimZMax) : | real(8),intent(out) |
Takemi(2007) での湿度場の与え方の都合上, 高度1.5km以下と高度1.5km以上で 相対湿度を別々に求めた 前者は HUM_Takemi2007 or HUM_Takemi2007_TDRY or HUM_Takemi2007_MDRY で求め 後者は ECCM_Takemi2007 で求めている. そのため, 全高度での相対湿度 z_RHall を確認のためにも求めるために, このサブルーチンで両者を結合している 値の確認が主な目的なため, 基本的に必要ないが, 仮にECCM_MolFrが必要なときに 相対湿度が必要なこともあり, 一応求めている
subroutine HUM_Takemi2007ALL(z_RHall)
!
!== 概要
! Takemi(2007) での湿度場の与え方の都合上, 高度1.5km以下と高度1.5km以上で
! 相対湿度を別々に求めた
! 前者は HUM_Takemi2007 or HUM_Takemi2007_TDRY or HUM_Takemi2007_MDRY で求め
! 後者は ECCM_Takemi2007 で求めている.
! そのため, 全高度での相対湿度 z_RHall を確認のためにも求めるために,
! このサブルーチンで両者を結合している
! 値の確認が主な目的なため, 基本的に必要ないが, 仮にECCM_MolFrが必要なときに
! 相対湿度が必要なこともあり, 一応求めている
implicit none
real(8),intent(out) :: z_RHall(DimZMin:DimZMax)
real(8) :: a_MolFrIn(1:SpcNum)
real(8) :: Humidit
real(8) :: z_Tem(DimZMin:DimZMax)
real(8) :: z_Pre(DimZMin:DimZMax)
real(8) :: z_MolMea(DimZMin:DimZMax)
real(8) :: za_MolF(DimZMin:DimZMax)
real(8) :: ZTR
real(8) :: z_humidity1(DimZMin:DimZMax) !高度0 - 1.5 km の相対湿度
real(8) :: z_humidity2(DimZMin:DimZMax) !高度1.5 km 以上の相対湿度
integer :: k, T, L
!高度1.5km以下の湿度を取得
call ECCM_Takemi2007(z_Tem, z_Pre, za_MolF, z_humidity1)
!高度1.5km以上の湿度を取得
!三つから選択
call HUM_Takemi2007(z_humidity2, ZTR, T, L)
! call HUM_Takemi2007_TDRY(z_humidity2, ZTR, T, L)
! call HUM_Takemi2007_MDRY(z_humidity2, ZTR, T, L)
do k = RegZMin+1, DimZMax
if (k.le.L) then
z_RHall(k) = z_humidity1(k)
else if (k.gt.L) then
z_RHall(k) = z_humidity2(k)
end if
end do
end subroutine HUM_Takemi2007ALL
| Subroutine : | |||
| z_RH(DimZMin:DimZMax) : | real(8),intent(out)
| ||
| Ztr : | real(8),intent(out)
| ||
| TR : | integer,intent(out)
| ||
| LL : | integer,intent(out)
|
subroutine HUM_Takemi2007_MDRY(z_RH,Ztr,TR,LL)
!
!== 概要
! * deepconv の地球用のテスト計算としてTakemi(2007)の再現計算を
! するための湿度の基本場を作成する
! * Takemi(2007) の熱帯場と中緯度場の相対湿度の"DRY CASE"の
! 高度 1.5km 以上の湿度を計算している(詳細はTakemi(2007)の図2を参照のこと)
! * Takemi(2007)では高度1.5km以上はある関数で与えているが,
! 高度1.5km以下は一定の混合比で与えているため,
! 高度1.5km以下の湿度は温度や気圧が求まらないと得られないので後で求める
! * 熱帯場とは乾燥のさせ方が違う
implicit none
real(8),intent(out) :: z_RH(DimZMin:DimZMax) !Relative Humidity(相対湿度)
real(8),intent(out) :: Ztr !対流圏界面の高度
integer,intent(out) :: TR !対流圏界面の高度の格子番号
integer,intent(out) :: LL !高度 1.5 km 直下の格子番号 (Low Level)
integer :: k
z_RH = 0.0d0
Ztr = 0.0d0
! Ztrに対流圏界面高度の値を与える
! Takemi(2007)では12kmだが, 必ずしも格子点がちょうど12kmには乗らないため,
! 12kmか,12kmの直上の位置を対流圏界面にする
do k = RegZMin+1,DimZMax
! s_Z が ちょうど 12 km なら 12 km が対流圏界面
if (s_Z(k).eq.12000) then
Ztr = s_Z(k)
! s_Z がちょうど12kmではなかった場合,12kmの直上のs_Zが対流圏海面Ztrになる
else if (s_Z(k).gt.12000) then
Ztr = minval( s_Z, 1, s_Z > 12000 )
end if
end do
! 対流圏界面高度の格子点を求める (Ztr = s_Z(TR))
do k = RegZMin+1,DimZMax
if (s_Z(k).ge.12000.and.s_Z(k).le.(12000 + DelZ)) then
TR = k
end if
end do
! 高度1.5kmもしくはその直下の格子番号を求める
do k = RegZMin+1,DimZMax
if (s_Z(k).le.1500.and.s_Z(k).ge.(1500 - DelZ)) then
LL = k
end if
end do
! Takemi(2007)の式を用いて相対湿度を求める
! とりあえず高度1.5kmまでの湿度はゼロにしておく
! 乾燥大気のエントレインメントを考慮した中緯度場の"DRY CASE"では
! 高度 2.5 km で相対湿度を 13 % もしくは 30 % 減らす DRY1,2 ケースを用意する.
do k = RegZMin+1, DimZMax
if (k.le.LL) then
z_RH(k) = 0.0d0
! 1.5 km から2.5 km まではそのまま変化
else if (s_Z(k).lt.2500.and.k.gt.LL) then
z_RH(k) = 1.0d0 - 0.750d0 * (s_Z(k) / Ztr)**1.25
! 2.5 km 以降は 13 % or 30 % 減らす
else if (s_Z(k).ge.2500) then
! z_RH(k) = 1.0d0 - 0.750d0 * (s_Z(k) / Ztr)**1.25 - 0.13d0 ! DRY1
z_RH(k) = 1.0d0 - 0.750d0 * (s_Z(k) / Ztr)**1.25 - 0.30d0 ! DRY2
end if
end do
! 上空で相対湿度が 25% 以下になる場合は 25% で固定する
do k = RegZMin+1, DimZMax
if (z_RH(k).le.(0.25).and.k.gt.LL) then
z_RH(k) = 0.250d0
end if
end do
end subroutine HUM_Takemi2007_MDRY
| Subroutine : | |||
| z_RH(DimZMin:DimZMax) : | real(8),intent(out)
| ||
| Ztr : | real(8),intent(out)
| ||
| TR : | integer,intent(out)
| ||
| LL : | integer,intent(out)
|
subroutine HUM_Takemi2007_TDRY(z_RH,Ztr,TR,LL)
!
!== 概要
! * deepconv の地球用のテスト計算としてTakemi(2007)の再現計算を
! するための湿度の基本場を作成する
! * Takemi(2007) の熱帯場と中緯度場の相対湿度の"DRY CASE"の
! 高度 1.5km 以上の湿度を計算している(詳細はTakemi(2007)の図2を参照のこと)
! * Takemi(2007)では高度1.5km以上はある関数で与えているが,
! 高度1.5km以下は一定の混合比で与えているため,
! 高度1.5km以下の湿度は温度や気圧が求まらないと得られないので後で求める
! * 中緯度場とは乾燥のさせ方が違う
!
implicit none
real(8),intent(out) :: z_RH(DimZMin:DimZMax) !Relative Humidity(相対湿度)
real(8),intent(out) :: Ztr !対流圏界面の高度
integer,intent(out) :: TR !対流圏界面の高度の格子番号
integer,intent(out) :: LL !高度 1.5 km 直下の格子番号 (Low Level)
integer :: k
z_RH = 0.0d0
Ztr = 0.0d0
! Ztrに対流圏界面高度の値を与える
! Takemi(2007)では12kmだが, 必ずしも格子点がちょうど12kmには乗らないため,
! 12kmか,12kmの直上の位置を対流圏界面にする
do k = RegZMin+1,DimZMax
! s_Z が ちょうど 12 km なら 12 km が対流圏界面
if (s_Z(k).eq.12000) then
Ztr = s_Z(k)
! s_Z がちょうど12kmではなかった場合,12kmの直上のs_Zが対流圏界面Ztrになる
else if (s_Z(k).gt.12000) then
Ztr = minval( s_Z, 1, s_Z > 12000 )
end if
end do
! 対流圏界面高度の格子点を求める (Ztr = s_Z(TR))
do k = RegZMin+1,DimZMax
if (s_Z(k).ge.12000.and.s_Z(k).le.(12000 + DelZ)) then
TR = k
end if
end do
! 高度1.5kmもしくはその直下の格子番号を求める
do k = RegZMin+1,DimZMax
if (s_Z(k).le.1500.and.s_Z(k).ge.(1500 - DelZ)) then
LL = k
end if
end do
! Takemi(2007)の式を用いて相対湿度を求める
! とりあえず高度1.5km以下の湿度はゼロにしておく
! さらに, 乾燥大気のエントレインメントを考慮した熱帯の"DRY CASE"では
! 高度 2.5 km, 5.0 km, 7.5 km のどこかで相対湿度を 20 % 減らすDRY1,2,3ケースを用意する
do k = RegZMin+1, DimZMax
if (k.le.LL) then
z_RH(k) = 0.0d0
! 1.5 km から (2.5 km, 5.0 km, 7.5 km) まではそのまま変化
! else if (s_Z(k).lt.2500.and.k.gt.LL) then ! DRY1
! else if (s_Z(k).lt.5000.and.k.gt.LL) then ! DRY2
else if (s_Z(k).lt.7500.and.k.gt.LL) then ! DRY3
z_RH(k) = 1.0d0 - 0.750d0 * (s_Z(k) / Ztr)**1.25
! (2.5 km, 5.0 km, 7.5 km)で20 % 湿度を減少させる
! else if (s_Z(k).ge.2500) then ! DRY1
! else if (s_Z(k).ge.5000) then ! DRY2
else if (s_Z(k).ge.7500) then ! DRY3
z_RH(k) = 1.0d0 - 0.750d0 * (s_Z(k) / Ztr)**1.25 - 0.20d0
end if
end do
! 上空で相対湿度が 25% 以下になる場合は 25% で固定する
do k = RegZMin+1, DimZMax
if (z_RH(k).le.(0.25).and.k.gt.LL) then
z_RH(k) = 0.250d0
end if
end do
end subroutine HUM_Takemi2007_TDRY