!= 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 !case 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, &!追加(ima) & CpDry !追加(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 ECCM_Ima public HUM public HUM_Takemi2007 public HUM_Takemi2007_TDRY public HUM_Takemi2007_MDRY public ECCM_Takemi2007 public HUM_Takemi2007ALL public Height 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 から相対湿度を読み込んで代入する用 real(8),intent(out) :: z_humid(DimZMin:DimZMax) real(8),intent(out) :: kansoku_altitude(1:88) !観測データの間隔 real(8),intent(out) :: kansoku_HumZ(1:88) !相対湿度の観測値(BasicZ の補間用) integer :: k,m z_humid = 0.0d0 kansoku_altitude = 0.0d0 kansoku_HumZ = 0.0d0 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 は二次元データ (x,z) なので --!! !!-- 一次元にする --!! do m=1,88 kansoku_HumZ(m) = Observed_Humidity(1,m) end do !!-- Ob_HumZ から線形補完して HumZ に値を与える --!! ! z_HumZ(RegZMin) = 83.0d0 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 の基本場を再現するプログラム. ! !! 相対湿度は HUM から引いてきてモル比に変換する. ! !! 気温はこのプログラム内で, 外部の netCDF ファイルから持ってくる. ! !! このようにしたのは相対湿度が正しく入っているかを 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) :: SatPress !飽和蒸気圧 ! real(8) :: VapPress !蒸気圧 ! real(8) :: DelMolFr ! real(8) :: a_MolFr(SpcNum) !モル比の作業配列 real(8) :: Ob_altitude(1:88) !観測データの間隔 real(8) :: Ob_TempZ(1:89) !気温の観測値(BasicZ の補間用) real(8) :: Ob_HumZ(1:88) !相対湿度の観測値(BasicZ の補間用) real(8) :: z_PressSaturated(DimZMin:DimZMax) real(8) :: z_HumZ(DimZMin:DimZMax) real(8) :: Observed_Temp(1:100,1:89) !NetCDF から気温を読み込んで代入する用 real(8) :: Observed_Humidity(1:100,1:88)!NetCDF から相対湿度を読み込んで代入する用 integer :: k, s, i, m, b z_HumZ = 0.0d0 call HUM(z_HumZ, Ob_altitude, Ob_HumZ) ! do m = 1,88 ! Ob_altitude(m) = 25.0d1*(m-1) ! 観測データ内での高度 ! end do !!-- NetCDF ファイルからの気温と相対湿度の読み込み --!! call HistoryGet('Ob_Temp.nc', 'Ob_Temp', array=Observed_Temp) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!---- 気温を求める ----!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!-- 別プログラムで作成した Ob_Temp.nc は二次元データ (x,z) なので --! !!-- 一次元 (z) にする --! do m=1,89 Ob_TempZ(m) = Observed_Temp(1,m) end do !!-- Ob_TempZ から線形補完して Temp に値を与える --!! z_Temp = 1.0d-60 ! 初期化 do k = RegZMin+1,DimZMax do m = 1,88 !! s_Z が Ob_altitude のある区間に挟まれるとき, !! その区間の Obaltitude(m), Obaltitude(m+1) を結ぶ !! 直線で Temp を線形補間する 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(i) の時は, Ob_TempZ(i) = 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 ! 初期化 !!!--地表面気圧・温度から大気最下層の気圧を求める--!!! z_Press(RegZMin+1) = & & PressSfc - (Grav * PressSfc * DelZ) & & /((GasRUniv/ & & (MolWtMeanDry + (MolWtMeanWet - MolWtMeanDry) & & *SvapPress(6,TempSfc) & & *z_HumZ(RegZMin+1) / PressSfc)) & !! 地表面の相対湿度は大気最下層のものを代用 & *TempSfc) ! z_Press(RegZMin+1) = 997.750d2 ! 高度 125 m の気圧. (地表面気圧は 1010 hPa) ! za_MolFr(RegZMin+1,1:SpcNum) = a_MolFrIni(1:SpcNum) ! 地表面でのモル比, SpcWetMolFr から与えられる ! MolWtMeanDry = MolWtDry * (1.0d0 - sum(a_MolFrIni)) ! 地表面での乾燥大気分子量 ! MolWtMeanWet = dot_product(MolWtWet, a_MolFrIni) ! 地表面での湿潤大気分子量 za_MolFr(RegZMin+1,1) = SvapPress( 6,z_Temp(RegZMin+1)) * & ! 大気最下層のモル比 & (z_HumZ(RegZMin+1)/100)/z_Press(RegZMin+1) ! SpcWetMolFr からは与えない MolWtMeanDry = MolWtDry * 1.0d0 - za_MolFr(RegZMin+1,1)! 地表面での乾燥大気分子量 MolWtMeanWet = MolWtWet(1) * za_MolFr(RegZMin+1,1) ! 地表面での湿潤大気分子量 z_MolWtMean(RegZMin+1) = MolWtMeanDry + MolWtMeanWet 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 ! この飽和蒸気圧は Tetens(1930) を使用. ! ただし, Tetens(1930) の式中の温度は [℃] であり, また [hPa] で出てくるので ! 適宜変えている. do k = RegZMin+1,DimZMax-1 ! do s = 1,CondNum ! ! za_MolFr(k,IdxCG(s)) = SvapPress( SpcWetID(IdxCC(s)),z_Temp(k))* & ! & (z_HumZ(k)/100)/z_Press(k) ! end do za_MolFr(k,1) = SvapPress( 6,z_Temp(k))* (z_HumZ(k)/100)/z_Press(k) ! MolWtMeanDry = MolWtDry * (1.0d0 - sum(za_MolFr(k,1:SpcNum))) ! MolWtMeanWet = dot_product(MolWtWet, za_MolFr(k,1:SpcNum)) MolWtMeanDry = MolWtDry * (1.0d0 - za_MolFr(k,1)) MolWtMeanWet = MolWtWet(1) * za_MolFr(k,1) z_MolWtMean(k) = 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 !!!!!!!! !!!!!!!!MolWtMean を求める !!!!!!!!但し, 実際は MolWtMeanDry と MolWtMeanWet を使用する !!!!!!!!ため不要... !!!!!!!! ! do k = RegZMin+1,DimZMax-1 ! z_MolWtMean(k) = MolWtMeanDry + MolWtMeanWet ! end do !!!!!!!!! !!!!!!!!! end subroutine ECCM_Ima !!!------------------------------------------------------------------------------!!! subroutine Height(B,U) !!!!! プログラムが正しいか確認するためだけのもの !!!!!! !!!!! 特に意味はなし. !!!!!! implicit none integer,intent(out) :: B integer,intent(out) :: U integer :: k do k = RegZMin+1, DimZMax if (s_Z(k).ge.2500.and.s_Z(k).le.(2500 + DelZ)) then ! if (s_Z(k).ge.2500.and.s_Z(k+1).le.(5000 + 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 ! 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" 場を再現. ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! implicit none real(8),intent(out) :: z_RH(DimZMin:DimZMax) !Relative Humidity, Takemi, 2007 での名前 real(8),intent(out) :: Ztr !対流圏界面の高度. Takemi, 2007 の z_tr に相当 integer,intent(out) :: TR !対流圏界面の格子番号 integer,intent(out) :: LL !高度 1.5 km 直下の格子番号 (Low Level) integer :: k z_RH = 0.0d0 Ztr = 0.0d0 !!! まず対流圏界面高度を与える. Takemi, 2007 では 12 km だが, 必ずしも格子点が !!! ちょうど 12 km には乗らないので, 12 km か, 12 km の直上の格子点の位置を対流圏界面にする do k = RegZMin+1,DimZMax !! s_Z が ちょうど 12 km なら 12 km が対流圏界面 if (s_Z(k).eq.12000) then Ztr = s_Z(k) !! s_Z がちょうど 12 km にならなかったら 12 km より大きい s_Z の中の最小の s_Z を対流圏界面にする 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 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.5 km より上空の相対湿度を求めるものである. do k = RegZMin+1, DimZMax if (k.le.LL) then z_RH(k) = 0.0d0 ! if (s_Z(k).le.Ztr.and.s_Z(k).gt.1500) then else if (k.le.TR.and.k.gt.LL) then z_RH(k) = 1.0d0 - 0.750d0 * (s_Z(k) / Ztr)**1.25 !! (4/5) ではべき乗が認識されない !! なんで? ! else if (s_Z(k).gt.Ztr) then 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" 場を再現. ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! implicit none real(8),intent(out) :: z_RH(DimZMin:DimZMax) !Relative Humidity, Takemi, 2007 での名前 real(8),intent(out) :: Ztr !対流圏界面の高度. Takemi, 2007 の z_tr に相当 integer,intent(out) :: TR !対流圏界面の格子番号 integer,intent(out) :: LL !高度 1.5 km 直下の格子番号 (Low Level) integer :: k z_RH = 0.0d0 Ztr = 0.0d0 !!! まず対流圏界面高度を与える. Takemi, 2007 では 12 km だが, 必ずしも格子点が !!! ちょうど 12 km には乗らないので, 12 km か, 12 km の直上の格子点の位置を対流圏界面にする do k = RegZMin+1,DimZMax !! s_Z が ちょうど 12 km なら 12 km が対流圏界面 if (s_Z(k).eq.12000) then Ztr = s_Z(k) !! s_Z がちょうど 12 km にならなかったら 12 km より大きい s_Z の中の最小の s_Z を対流圏界面にする 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 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.5 km より上空の相対湿度を求めるものである. !! さらに, 乾燥大気のエントレインメントを考慮した熱帯の DRY の !! 湿度場では, 高度 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 % 以下にしない. 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" 場を再現. ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! implicit none real(8),intent(out) :: z_RH(DimZMin:DimZMax) !Relative Humidity, Takemi, 2007 での名前 real(8),intent(out) :: Ztr !対流圏界面の高度. Takemi, 2007 の z_tr に相当 integer,intent(out) :: TR !対流圏界面の格子番号 integer,intent(out) :: LL !高度 1.5 km 直下の格子番号 (Low Level) integer :: k z_RH = 0.0d0 Ztr = 0.0d0 !!! まず対流圏界面高度を与える. Takemi, 2007 では 12 km だが, 必ずしも格子点が !!! ちょうど 12 km には乗らないので, 12 km か, 12 km の直上の格子点の位置を対流圏界面にする do k = RegZMin+1,DimZMax !! s_Z が ちょうど 12 km なら 12 km が対流圏界面 if (s_Z(k).eq.12000) then Ztr = s_Z(k) !! s_Z がちょうど 12 km にならなかったら 12 km より大きい s_Z の中の最小の s_Z を対流圏界面にする 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 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.5 km より上空の相対湿度を求めるものである. !! さらに, 乾燥大気のエントレインメントを考慮した中緯度の DRY の !! 湿度場では, 高度 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 % 以下にしない. 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 ) 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 !湿潤気塊の分子量の作業変数 !! Takemi, 2007 を用いるために必要なものたち------- real(8) :: z_TakemiRH(DimZMin:DimZMax) 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) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!Takemi, 2007 による温位場!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 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 !!!!----大気最下層の気温・気圧・水蒸気混合比・分子量----!!!! !!!--地表面気圧・温度から大気最下層の気圧を求める--!!! 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) ! z_Press(RegZMin+1) = 987.750d2 ! 高度 125 m の気圧. (地表面気圧は 1000 hPa) z_Temp(RegZMin+1) = z_TakemiTheta(RegZMin+1) * (z_Press(RegZMin+1) / 100000.0d0)**(GasRDry / CpDry) 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)) !!!!--------------------------------------------!!!! do k = RegZMin+1, DimZMax !! 大気最下層から高度 1.5 km (k が L 以下) 以下の場 if (k.le.L-1) then za_MolFr(k+1,1) = 0.0180d0 * MolWtDry / MolWtWet(1) !! Takemi, 2007 では大気下層 1.5 km の混合比を ! !! 10 - 18 g/kg で変化させて与えている ! 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)) !!---低層で混合比を固定し続けると相対湿度が 100 % を超える可能性がある---!! !!--Takemi では暗黙の内に 100 % を超えたら前のステップの値にしている???--!! !!--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.5 km より上で, 対流圏界面以下 (L < k < T) の場 else if (k.le.T.and.k.gt.L) then z_TakemiRH(L) = z_LLhumidity(L) !この if 文の初めの k は L+1 から始まるが !z_TakemiRH(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 を求める. 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) real(8) :: z_humidity2(DimZMin:DimZMax) 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 ! デフォルト do k = RegZMin+1, 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) / 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, 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 * 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