!= Module ECCM ! ! Authors:: SUGIYAMA Koichiro, ODAKA Masatsugu ! Version:: $Id: eccm.f90,v 1.28 2010-08-11 07:34:19 sugiyama Exp $ ! Tag Name:: $Name: arare4-20100910 $ ! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved. ! License:: See COPYRIGHT[link:../../COPYRIGHT] ! !== Overview ! !断熱的に上昇する気塊の温度減率を計算し, 静水圧平衡から圧力を決める ! !== Error Handling ! !== Known Bugs ! !== Note ! ! * 比熱は定数, 平均分子量は配列 ! * オイラースキームだと精度が足りないので, ルンゲクッタスキームを利用. ! !== Future Plans ! module ECCM !モジュール読み込み use gtool_history !外部にあるnetCDFからデータを読み込むための !モジュール(for ECC_Ima) use dc_message, only: MessageNotify use gridset, only: DimZMin, &! 配列の Z 方向の下限 & DimZMax, &! 配列の Z 方向の上限 & RegZMin, &! & SpcNum, &! & s_Z, &! & DelZ ! use basicset,only: MolWtDry, &! & MolWtWet, &! & CpDryMol, &! & SpcWetID, &! & TempSfc, &! & PressSfc, &! & Grav, &! & GasRDry, &!追加(for ECC_Ima) & CpDry !追加(for ECC_Ima) use chemcalc, only: SvapPress, &! & LatentHeatPerMol, &! & ReactHeatNH4SHPerMol use moistset, only: CondNum, &!凝結過程の数 & IdxCG, &!凝結過程(蒸気)の配列添え字 & IdxCC, &!凝結過程(雲)の配列添え字 & GasNum, &!気体の数 & IdxNH3, &!NH3(蒸気)の配列添え字 & IdxH2S !H2S(蒸気)の配列添え字 use ChemData, only: GasRUniv use MoistFunc,only: DelMolFrNH4SH use StoreStab,only: StoreStabTemp, StoreStabMolWt !暗黙の型宣言禁止 implicit none !属性の指定 private !関数の公開 public ECCM_MolFr public ECCM_Stab public ECCM_Dry public ECCM_Wet public ECCM_Wet2 public HUM !ECCM_Imaで用いる相対湿度の基本場を作成 public ECCM_Ima !Yamasaki(1983)の温度と相対湿度で基本場を作成 public HUM_Takemi2007 !Takemi(2007)の相対湿度の基本場(一部)を作成 public HUM_Takemi2007_TDRY !Takemi(2007)の相対湿度の基本場(一部)を作成 public HUM_Takemi2007_MDRY !Takemi(2007)の相対湿度の基本場(一部)を作成 public ECCM_Takemi2007 !Takemi(2007)の温度と湿度で基本場を作成 public HUM_Takemi2007ALL !Takemi(2007)の相対湿度の基本場を作成(for ECCM_MolFr) public Height !Takemi(2007)基本場に関するチェック用(不要) contains !!!------------------------------------------------------------------------------!!! ! subroutine ECCM_Dry( a_MolFrIni, Humidity, z_Temp, z_Press, z_MolWtMean, za_MolFr ) 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 ECCM_Wet( a_MolFrIni, Humidity, z_Temp, z_Press, z_MolWtMean, za_MolFr ) 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 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 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) :: z_humid(DimZMin:DimZMax) !線形補間された相対湿度 real(8),intent(out) :: kansoku_altitude(1:88) !観測データに対応する高度 real(8),intent(out) :: kansoku_HumZ(1:88) !相対湿度の観測値(1次元) integer :: k,m z_humid = 0.0d0 kansoku_altitude = 0.0d0 kansoku_HumZ = 0.0d0 ! 観測データに対応する高度(250m間隔) ! "88"はデータの数 ! 線形補間の際に使用する do m = 1,88 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 ECCM_Ima( a_MolFrIni, Humidity, z_Temp, z_Press, z_MolWtMean, za_MolFr ) ! !== 概要 ! * Yamasaki, 1983 の温度と相対湿度の観測値で基本場を作成する ! * 温度の基本場 ! * 観測データをNetCDFファイル化したものから値を読み込む ! * 読み込んだ値を線形補間して温度の基本場を作成する ! * 湿度の基本場 ! * 相対湿度は subroutine HUM で作成済み ! * 相対湿度をモル比に変換して湿度の基本場を作成する ! ! * 相対湿度だけ先に計算した目的は, "ECCM_Ima"の引数を ! 他の"ECCM_*"と同じにしたまま, basicenv.f90 で write し ! 値を確認するため ! * 気圧の基本場 ! * 湿潤成分と乾燥成分の分子量差を考慮した静水圧平衡から計算する ! implicit none real(8), intent(in):: a_MolFrIni(1:SpcNum) !下部境界でのモル比 real(8), intent(in) :: Humidity !相対湿度 ( Humidity <= 1.0 ) real(8), intent(out):: z_MolWtMean(DimZMin:DimZMax) 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) :: MolWtMeanDry !乾燥気塊の分子量の作業変数 real(8) :: MolWtMeanWet !湿潤気塊の分子量の作業変数 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 call HUM(z_HumZ, Ob_altitude, Ob_HumZ) ! 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) ! 地表面での乾燥大気の平均分子量を計算する MolWtMeanDry = MolWtDry * 1.0d0 - za_MolFr(RegZMin+1,1) ! 地表面での湿潤大気の平均分子量を計算する MolWtMeanWet = MolWtWet(1) * za_MolFr(RegZMin+1,1) z_MolWtMean(RegZMin+1) = MolWtMeanDry + MolWtMeanWet ! テスト計算で 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 ! 大気最下層の値から, 気圧と湿度の基本場を求める ! 気圧は湿潤成分と乾燥成分の分子量差を考慮した静水圧平衡から計算する ! 湿度はモル比で算出する ! 湿度の基本場の計算 do k = RegZMin+1,DimZMax-1 za_MolFr(k+1,1) = SvapPress( 6,z_Temp(k)) * (z_HumZ(k)/100)/z_Press(k) MolWtMeanDry = MolWtDry * (1.0d0 - za_MolFr(k+1,1)) MolWtMeanWet = MolWtWet(1) * za_MolFr(k+1,1) z_MolWtMean(k+1) = MolWtMeanDry + MolWtMeanWet ! 気圧の基本場の計算 z_Press(k+1) = & & z_Press(k)-(Grav*z_Press(k)*DelZ) & & /((GasRUniv/ & & (MolWtMeanDry + (MolWtMeanWet - MolWtMeanDry) & & *SvapPress(6,z_Temp(k)) & & *(z_HumZ(k)/100)/z_Press(k))) & & *z_Temp(k)) end do end subroutine ECCM_Ima !!!------------------------------------------------------------------------------!!! subroutine Height(B,U) ! !== 概要 ! * Takemi(2007)の基本場のチェックに用いる ! * Takemi(2007)では湿度の基本場が高度によって従う式が異なるため ! 式が変わる高度が正しく出力されているかを確認するために以下の ! サブルーチンを作成した ! * このサブルーチンはなくても計算に支障は無い ! * 削除する場合は basicenv.f90 で引用している部分も消す必要がある ! implicit none integer,intent(out) :: B integer,intent(out) :: U integer :: k ! 2.5km と 5.0km が式の変わり目 do k = RegZMin+1, DimZMax if (s_Z(k).ge.2500.and.s_Z(k).le.(2500 + DelZ)) then B = k end if end do do k = RegZMin+1, DimZMax if (s_Z(k).ge.5000.and.s_Z(k).le.(5000 + DelZ)) then U = k end if end do end subroutine Height !!!------------------------------------------------------------------------------!!! subroutine HUM_Takemi2007(z_RH,Ztr,TR,LL) ! !== 概要 ! * deepconv が地球を仮定した場合に正しく動いているか確認するために, ! Takemi(2007)を参考に, 大気が不安定・安定な場合の積雲対流の振る舞いを ! 見るための湿度の基本場を作成する ! * Takemi(2007) の熱帯場と中緯度場の相対湿度の"BASE CASE"の ! 高度 1.5km 以上を計算している(Takemi(2007)の図2を参照のこと) ! 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.5 km以下では任意の一定の混合比(q_v0)を与えるため, ! 以下の式は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 HUM_Takemi2007_TDRY(z_RH,Ztr,TR,LL) ! !== 概要 ! * deepconv が地球を仮定した場合に正しく動いているか確認するために, ! Takemi(2007)を参考に, 大気が不安定・安定な場合の積雲対流の振る舞いを ! 見るための湿度の基本場を作成する ! * Takemi(2007) の熱帯場の相対湿度の"DRY CASE"の ! 高度 1.5km 以上を計算している(Takemi(2007)の図2を参照のこと) ! 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.5 km以下では任意の一定の混合比(q_v0)を与えるため, ! 以下の式は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 !!!------------------------------------------------------------------------------!!! subroutine HUM_Takemi2007_MDRY(z_RH,Ztr,TR,LL) ! !== 概要 ! * deepconv が地球を仮定した場合に正しく動いているか確認するために, ! Takemi(2007)を参考に, 大気が不安定・安定な場合の積雲対流の振る舞いを ! 見るための湿度の基本場を作成する ! * Takemi(2007) の中緯度場の相対湿度の"DRY CASE"の ! 高度 1.5km 以上を計算している(Takemi(2007)の図2を参照のこと) ! 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.5 km以下では任意の一定の混合比(q_v0)を与えるため, ! 以下の式は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 ECCM_Takemi2007( a_MolFrIni, Humidity, z_Temp, z_Press, z_MolWtMean, & &za_MolFr, z_LLhumidity ) ! !== 概要 ! * deepconv が地球を仮定した場合に正しく動いているか確認するために, ! Takemi(2007)を参考に, 大気が不安定・安定な場合の積雲対流の振る舞いを ! 見るための湿度の基本場を作成する ! * 基本場の温度の式が温位で与えられているため, 温度に変換する必要がある implicit none real(8), intent(in):: a_MolFrIni(1:SpcNum) !下部境界でのモル比 real(8), intent(in) :: Humidity !相対湿度 ( Humidity <= 1.0 ) real(8), intent(out):: z_MolWtMean(DimZMin:DimZMax) 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 までの相対湿度 real(8) :: MolWtMeanDry !乾燥気塊の分子量の作業変数 real(8) :: MolWtMeanWet !湿潤気塊の分子量の作業変数 ! ! 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) ! 温度・湿度・気圧を求める 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 ! 大気最下層の気温・気圧・水蒸気混合比・分子量の計算 ! 地表面のモル比 ! 熱帯場は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 ! 地表面の大気の平均分子量 MolWtMeanDry = MolWtDry * (1 - a_MolFrSfc) MolWtMeanWet = MolWtWet(1) * a_MolFrSfc ! 地表面の相対湿度 ! モル比から逆算 HumSfc = a_MolFrSfc * PressSfc / SvapPress(6,TempSfc) ! 大気最下層の気圧 z_Press(RegZMin+1) = & & PressSfc - (Grav * PressSfc * DelZ) & & /((GasRUniv/ & & (MolWtMeanDry + (MolWtMeanWet - MolWtMeanDry) & & *SvapPress(6,TempSfc) & & *HumSfc / PressSfc)) & & *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 ! 大気最下層の大気の平均分子量 MolWtMeanDry = MolWtDry * (1 - za_MolFr(RegZMin+1,1)) MolWtMeanWet = MolWtWet(1) * za_MolFr(RegZMin+1,1) ! 大気最下層の相対湿度 ! モル比から逆算 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 ! 大気最下層から高度 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/ & & (MolWtMeanDry + (MolWtMeanWet - MolWtMeanDry) & & *SvapPress(6,z_Temp(k)) & & *z_LLhumidity(k)/z_Press(k))) & & *z_Temp(k)) ! 大気の平均分子量の計算 MolWtMeanDry = MolWtDry * (1 - za_MolFr(k+1,1)) MolWtMeanWet = MolWtWet(1) * za_MolFr(k+1,1) z_MolWtMean = MolWtMeanDry + MolWtMeanWet ! 温位の基本場から温度の基本場を逆算 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%を超える高度がある ! * Takemi では暗黙の内に 100% を超えたら前のステップの値にしているように思える ! * 相対湿度をプロットした図をTakemi(2007)の図2と比較するとよくわかる ! * MQ16 は 湿度 99.* % という値が出て, ほぼ 100% になってしまうので ! MQ18 とMQ14 の 100% を超える直前の値の間を取った値を入れておく ! * 苦肉の折衷策 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 ! 高度 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/ & & (MolWtMeanDry + (MolWtMeanWet - MolWtMeanDry) & & *SvapPress(6,z_Temp(k-1)) & & *z_TakemiRH(k-1)/z_Press(k-1))) & & *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) MolWtMeanDry = MolWtDry * (1.0d0 - za_MolFr(k,1)) MolWtMeanWet = MolWtWet(1) * za_MolFr(k,1) ! 対流圏界面より上 (k > T) の場 else if (k.gt.T) then ! 気圧の基本場 z_Press(k) = & & z_Press(k-1)-(Grav*z_Press(k-1)*DelZ) & & /((GasRUniv/ & & (MolWtMeanDry + (MolWtMeanWet - MolWtMeanDry) & & *SvapPress(6,z_Temp(k-1)) & & *z_TakemiRH(k-1)/z_Press(k-1))) & & *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) MolWtMeanDry = MolWtDry * (1.0d0 - za_MolFr(k,1)) MolWtMeanWet = MolWtWet(1) * za_MolFr(k,1) end if end do end subroutine ECCM_Takemi2007 !!!------------------------------------------------------------------------------!!! subroutine HUM_Takemi2007ALL(z_RHall) ! !== 概要 ! Takemi, 2007 の相対湿度の式は高度 1.5 km 以降の領域に適用するものであり, ! 地表から 1.5 km までの相対湿度は混合比の任意の一定値が与えられ, ! 混合比(モル比に変換)と温度に対する飽和蒸気圧・気圧から求めなくてはならない ! 前者は HUM_Takemi2007 で, 後者は ECCM_Takemi2007 の LLhumidity で求めている. ! そのため, 全高度での相対湿度 z_RHall を求めるために, このサブルーチンで ! 両者を結合している ! 値の確認が主な目的なため, 基本的に必要ないが, 水蒸気の初期擾乱を与えたい場合は ! ECCM_MolFr で相対湿度が必要なため, 一応求めている implicit none real(8),intent(out) :: z_RHall(DimZMin:DimZMax) ! real(8),intent(in) :: a_MolFrIn(1:SpcNum) ! real(8),intent(in) :: Humidity 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 call ECCM_Takemi2007(a_MolFrIn, Humidit, z_Tem, z_Pre, z_MolMea, za_MolF, z_humidity1) 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 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 real(8) :: Ob_altitude(1:88) !観測データの間隔 real(8) :: Ob_HumZ(1:88) !相対湿度の観測値(BasicZ の補間用) real(8) :: z_HumZ(DimZMin:DimZMax) integer :: k, s z_HumZ = 0.0d0 call HUM(z_HumZ, Ob_altitude, Ob_HumZ) ! call HUM_Takemi2007ALL(z_HumZ) !----------------------------------------------------------- ! 配列の初期化 !----------------------------------------------------------- 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)) = & za_MolFr(k,1) = & & min( & ! & za_MolFr(k-1,IdxCG(s)), & & za_MolFr(k-1,1), & & SvapPress(6, z_Temp(k) ) & ! & SvapPress( SpcWetID(IdxCC(s)), z_Temp(k) ) & ! & * Humidity / z_Press(k) & & * (z_HumZ(k)/100) / 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 ECCM_DTempDZ( Temp, Press, MolFr, DTempDZ ) !暗黙の型宣言禁止 implicit none !変数定義 real(8), intent(in) :: Temp real(8), intent(in) :: Press ! real(8), intent(inout) :: MolFr(0:SpcNum) !モル分率 real(8), intent(inout) :: MolFr(1:SpcNum) !モル分率 real(8), intent(out):: DTempDZ real(8) :: MolFrOrig(1:SpcNum) real(8) :: ReactHeat real(8) :: Heat(SpcNum) real(8) :: DelMolFr real(8) :: SatPress real(8) :: VapPress real(8) :: Humidity real(8) :: A, B integer :: s !初期化 DTempDZ = 0.0d0 ReactHeat = 0.0d0 Heat = 0.0d0 DelMolFr = 0.0d0 SatPress = 0.0d0 VapPress = 0.0d0 MolFrOrig = MolFr !------------------------------------------------------------ !NH4SH 以外の化学種の平衡条件 !------------------------------------------------------------ do s = 1, CondNum !飽和蒸気圧 SatPress = SvapPress( SpcWetID(IdxCC(s)), Temp ) !潜熱. Heat(IdxCG(s)) = LatentHeatPerMol( SpcWetID(IdxCC(s)), Temp ) !元々のモル分率を用いて現在の蒸気圧を計算 VapPress = MolFr(IdxCG(s)) * Press !飽和蒸気圧から凝結の有無を決める if ( VapPress < SatPress ) then !凝結していないので潜熱なし. Heat(IdxCG(s)) = 0.0d0 else !飽和蒸気圧と圧力から現在のモル比を計算 MolFr(IdxCG(s)) = max(SatPress / Press, 1.0d-16) end if end do !------------------------------------------------------------ !NH4SH の平衡条件 !------------------------------------------------------------ if ( IdxNH3 /= 0 ) then Humidity = 1.0d0 DelMolFr = & & max ( & & DelMolFrNH4SH( & & Temp, Press, MolFr(IdxNH3), MolFr(IdxH2S),& & Humidity & & ), & & 0.0d0 & & ) MolFr(IdxNH3) = MolFr(IdxNH3) - DelMolFr MolFr(IdxH2S) = MolFr(IdxH2S) - DelMolFr ReactHeat = ReactHeatNH4SHPerMol * DelMolFr end if !------------------------------------------------------------ !温度勾配を計算 !------------------------------------------------------------ !係数. 温度 Temp(i) で評価 A = dot_product( Heat(1:SpcNum), MolFrOrig(1:SpcNum)) & ! A = dot_product( Heat(1:SpcNum), MolFr(1:SpcNum)) & & / ( GasRUniv * Temp ) B = dot_product(( Heat(1:SpcNum) ** 2.0d0), MolFrOrig(1:SpcNum)) & ! B = dot_product(( Heat(1:SpcNum) ** 2.0d0), MolFr(1:SpcNum)) & & / ( CpDryMol * GasRUniv * ( Temp ** 2.0d0 ) ) !断熱温度減率 DTempDZ = - Grav * MolWtDry * (1.0d0 + A) / (CpDryMol * (1.0d0 + B)) & & + ReactHeat / (CpDryMol * DelZ) end subroutine ECCM_DTempDZ subroutine ECCM_Stab( xz_PotTemp, xz_Exner, xza_MixRt ) ! & xz_Stab, xz_StabTemp, xz_StabMolWt ) use gridset, only: DimXMin, &! 配列の X 方向の下限 & DimXMax, &! 配列の X 方向の上限 & DimZMin, &! 配列の Z 方向の下限 & DimZMax, &! 配列の Z 方向の上限 & 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, sg 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 * xz_EffMolWtBasicZ / CpDry ) xz_StabMolWt = & & - Grav * xz_avr_xr( xr_dz_xz( xz_MolWtWet ) ) & & / ( MolWtDry * xz_EffMolWtBasicZ ) xz_Stab = xz_StabTemp + xz_StabMolWt ! xz_Stab = & ! & Grav / xz_TempAll & ! & * ( xz_avr_xr( xr_dz_xz( xz_TempAll ) ) & ! & + Grav * xz_EffMolWtBasicZ / CpDry ) & ! & - Grav * xz_avr_xr( xr_dz_xz( xz_MolWtWet ) ) & ! & / ( MolWtDry * xz_EffMolWtBasicZ ) 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 end module ECCM