!= Module WarmRainPrm ! ! Authors:: SUGIYAMA Ko-ichiro, ODAKA Masatsugu ! Version:: $Id: warmrainprm.f90,v 1.23 2007/08/24 07:51:03 sugiyama Exp $ ! Tag Name:: $Name: arare4-20071012 $ ! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved. ! License:: See COPYRIGHT[link:../../COPYRIGHT] ! !== Overview ! !暖かい雨のバルク法を用いた, 水蒸気と雨, 雲と雨の混合比の変換係数を求める. ! * 中島健介 (1994) で利用した定式をそのまま利用. ! !== Error Handling ! !== Bugs ! !== Note ! !== Future Plans ! ! module WarmRainPrm ! !暖かい雨のバルク法を用いた, 水蒸気と雨, 雲と雨の混合比の変換係数を求める. ! * 中島健介 (1994) で利用した定式をそのまま利用. ! !モジュール読み込み use dc_message, only: MessageNotify use gridset, only: DimXMin, &! x 方向の配列の下限 & DimXMax, &! x 方向の配列の上限 & DimZMin, &! z 方向の配列の下限 & DimZMax, &! z 方向の配列の上限 & SpcNum !化学種の数 use basicset, only: PressBasis, &!温位の基準圧力 & CpDry, &!乾燥成分の比熱 & MolWtWet, &! & MolWtDry, &! & SpcWetID, &! & xz_DensBasicZ, &!基本場の密度 & xz_PotTempBasicZ, &!基本場の温位 & xz_ExnerBasicZ, &!基本場の無次元圧力 & xz_EffMolWtBasicZ, &!基本場の分子量の寄与 & xza_MixRtBasicZ, &!基本場の混合比 & GasRDry !乾燥成分の気体定数 use moistset, only: CondNum, &!凝結過程の数 & IdxCG, &!凝結過程(蒸気)の配列添え字 & IdxCC, &!凝結過程(雲)の配列添え字 & IdxCR, &!凝結過程(雲)の配列添え字 & CloudNum, &!雲の数 & RainNum, &!雨の数 & IdxC, &!雲の配列添え字 & IdxR, &!雨の配列添え字 & IdxNH3, &!NH3(蒸気)の配列添え字 & IdxH2S, &!H2S(蒸気)の配列添え字 & IdxNH4SHr !NH4SH(雨)の配列添え字 use average, only: xz_avr_xr use differentiate_center2, only: xr_dz_xz use ChemCalc, only: xz_SvapPress, xz_LatentHeat, ReactHeatNH4SH use MoistFunc,only: xz_DelMixRtNH4SH use StoreMixRt, only: StoreMixRtRain !暗黙の型宣言禁止 implicit none !属性の指定 private !関数を public にする public WarmRainPrm_Init public xz_Rain2GasHeat public xz_Rain2GasHeatNH4SH public xza_Rain2Gas public xza_Rain2GasNH4SH public xza_Cloud2Rain public xza_FallRain real(8) :: FactorJ = 1.0d0 !雲物理過程のパラメータ !木星では 3.0d0 !地球では 1.0d0 とする real(8) :: AutoConvTime = 1.0d3 !併合成長の時定数 [sec] real(8) :: MixRt_AutoConvCr = 1.0d-3 !併合成長を生じる臨界混合比 [kg/kg] save FactorJ, AutoConvTime, MixRt_AutoConvCr contains !!!=================================================================================!!! subroutine WarmRainPrm_Init( cfgfile ) !暗黙の型宣言禁止 implicit none !入力変数 character(*), intent(in) :: cfgfile !----------------------------------------------------------- ! NAMELIST から情報を取得 !----------------------------------------------------------- ! NAMELIST /warmrainprm/ & & FactorJ, AutoConvTime, MixRt_AutoConvCr open (10, FILE=cfgfile) read(10, NML=warmrainprm) close(10) call MessageNotify( "M", & & "WarmRainPrm_Init", "FactorJ = %f", d=(/FactorJ/) ) call MessageNotify( "M", & & "WarmRainPrm_Init", "AutoConvTime = %f", d=(/AutoConvTime/) ) call MessageNotify( "M", & & "WarmRainPrm_Init", "MixRt_AutoConvCr = %f", d=(/MixRt_AutoConvCr/) ) end subroutine WarmRainPrm_Init !!!=================================================================================!!! function xza_Rain2Gas(xz_Exner, xz_PotTemp, xza_MixRt, DelTime) ! ! 雨粒から蒸気への変換量を計算するためのルーチン ! ! 変換量および, 蒸気と雨粒の混合比は正の量なので, 計算の途中途中で ! 値が正になることを保証している. また, 元々存在する以上の雨粒が ! 蒸気に変換されないように, 元々の雨粒混合比を変換量の上限としている. ! ! 木星の場合は, FactorJ で変換量を加速する. ! !暗黙の型宣言禁止 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) :: DelTime !時間刻み real(8) :: xza_Rain2Gas(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) ! real(8) :: xz_TempAll(DimXMin:DimXMax, DimZMin:DimZMax) !温度の擾乱成分 + 平均成分 real(8) :: xz_PressAll(DimXMin:DimXMax, DimZMin:DimZMax) !全圧 real(8) :: xza_MixRtAll(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) !混合比の擾乱成分 + 平均成分 real(8) :: xz_NonSaturate(DimXMin:DimXMax, DimZMin:DimZMax) !未飽和度(飽和混合比と蒸気の混合比の差) integer :: s !温度, 圧力, 混合比の全量を求める !擾乱成分と平均成分の足し算 xz_TempAll = ( xz_PotTemp + xz_PotTempBasicZ ) * ( xz_Exner + xz_ExnerBasicZ ) xz_PressAll = PressBasis * ((xz_Exner + xz_ExnerBasicZ) ** (CpDry / GasRDry)) xza_Rain2Gas = 0.0d0 !混合比は正の値であることを保証 !移流拡散計算で負になることがあり得るので. xza_MixRtAll = max( 0.0d0, xza_MixRt + xza_MixRtBasicZ ) do s = 1, CondNum !飽和蒸気圧と混合比の差(飽和度)を計算. ! 雨から蒸気への変換量は飽和度に比例する. xz_NonSaturate = & & max( & & 0.0d0, & & xz_SvapPress(SpcWetID(IdxCC(s)), xz_TempAll) & & * MolWtWet(IdxCG(s)) & & / ( MolWtDry * xz_PressAll) & & - xza_MixRtAll(:,:,IdxCG(s)) & & ) !雨の変換量 ! 元々の雨粒の混合比以上に蒸発が生じないように上限値を設定 xza_Rain2Gas(:,:,IdxCR(s)) = & & - min( & & DelTime * 4.85d-2 * FactorJ * xz_NonSaturate & & * ( xza_MixRtAll(:,:,IdxCR(s)) * xz_DensBasicZ )** 0.65d0, & & xza_MixRtAll(:,:,IdxCR(s)) & & ) !蒸気の変換量 ! 雨粒の変換量とは符号が逆となる xza_Rain2Gas(:,:,IdxCG(s)) = - xza_Rain2Gas(:,:,IdxCR(s)) end do end function xza_Rain2Gas !!!=================================================================================!!! function xza_Rain2GasNH4SH(xz_Exner, xz_PotTemp, xza_MixRt, DelTime) ! ! 雨粒から蒸気への変換量を計算するためのルーチン ! ! 変換量および, 蒸気と雨粒の混合比は正の量なので, 計算の途中途中で ! 値が正になることを保証している. また, 元々存在する以上の雨粒が ! 蒸気に変換されないように, 元々の雨粒混合比を変換量の上限としている. ! !暗黙の型宣言禁止 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(in) :: DelTime !時間刻み real(8) :: xza_Rain2GasNH4SH(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) ! real(8) :: xz_TempAll(DimXMin:DimXMax, DimZMin:DimZMax) !温度の擾乱成分 + 平均成分 real(8) :: xz_PressAll(DimXMin:DimXMax, DimZMin:DimZMax) !圧力の擾乱成分 + 平均成分 real(8) :: xza_MixRtAll(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) !混合比の擾乱成分 + 平均成分 real(8) :: xz_NonSaturate(DimXMin:DimXMax, DimZMin:DimZMax) !未飽和度(飽和混合比と蒸気の混合比の差) !温度, 圧力, 混合比の全量を求める !擾乱成分と平均成分の足し算 xz_TempAll = ( xz_PotTemp + xz_PotTempBasicZ ) * ( xz_Exner + xz_ExnerBasicZ ) xz_PressAll = PressBasis * ((xz_Exner + xz_ExnerBasicZ) ** (CpDry / GasRDry)) xza_Rain2GasNH4SH = 0.0d0 !混合比は正の値であることを保証 !移流拡散計算で負になることがあり得るので. xza_MixRtAll = max( 0.0d0, xza_MixRt + xza_MixRtBasicZ ) !飽和蒸気圧と混合比の差(飽和度)を計算. ! 雨から蒸気への変換量は飽和度に比例する. ! 未飽和度を求めたいので, マイナスをかけ算している ! (DelMixRtNH4SH は, NH4SH が増加する方向, すなわち飽和度を正としている) xz_NonSaturate = & & max( & & 0.0d0, & & - xz_DelMixRtNH4SH( & & xz_TempAll, xz_PressAll, & & xza_MixRtAll(:,:,IdxNH3), xza_MixRtAll(:,:,IdxH2S), & & MolWtWet(IdxNH3), MolWtWet(IdxH2S) & & ) & & ) !雨の変換量 ! 元々の雨粒の混合比以上に蒸発が生じないように上限値を設定 xza_Rain2GasNH4SH(:,:,IdxNH4SHr) = & & - min( & & DelTime * 4.85d-2 * FactorJ * xz_NonSaturate & & * ( & & xza_MixRtAll(:,:,IdxNH4SHr) * xz_DensBasicZ & & ) ** 0.65d0, & & xza_MixRtAll(:,:,IdxNH4SHr) & & ) !蒸気の変換量 ! 雨粒の変換量とは符号が逆となる xza_Rain2GasNH4SH(:,:,IdxNH3) = & & - xza_Rain2GasNH4SH(:,:,IdxNH4SHr) * MolWtWet(IdxNH3) & & / MolWtWet(IdxNH4SHr) xza_Rain2GasNH4SH(:,:,IdxH2S) = & & - xza_Rain2GasNH4SH(:,:,IdxNH4SHr) * MolWtWet(IdxH2S) & & / MolWtWet(IdxNH4SHr) end function xza_Rain2GasNH4SH !!!=================================================================================!!! function xz_Rain2GasHeatNH4SH(xz_Exner, xza_DelMixRt) ! ! 雨粒から蒸気への変換量を計算するためのルーチン ! ! 変換量および, 蒸気と雨粒の混合比は正の量なので, 計算の途中途中で ! 値が正になることを保証している. また, 元々存在する以上の雨粒が ! 蒸気に変換されないように, 元々の雨粒混合比を変換量の上限としている. ! !暗黙の型宣言禁止 implicit none !変数定義 real(8), intent(in) :: xza_DelMixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) !混合比の変化量 real(8), intent(in) :: xz_Exner(DimXMin:DimXMax, DimZMin:DimZMax) !温度の擾乱成分 real(8) :: xz_Rain2GasHeatNH4SH(DimXMin:DimXMax, DimZMin:DimZMax) ! real(8) :: xz_ExnerAll(DimXMin:DimXMax, DimZMin:DimZMax) !圧力の全量を求める !擾乱成分と平均成分の足し算 xz_ExnerAll = xz_Exner + xz_ExnerBasicZ !雨から蒸気への相変化に伴う発熱 xz_Rain2GasHeatNH4SH = & & ReactHeatNH4SH * xza_DelMixRt(:,:,IdxNH4SHr) * xz_EffMolWtBasicZ & & / (xz_ExnerAll * CpDry) end function xz_Rain2GasHeatNH4SH !!!=================================================================================!!! function xz_Rain2GasHeat(xz_PotTemp, xz_Exner, xza_DelMixRt) ! ! 雨粒から蒸気への変換量を計算するためのルーチン ! ! 変換量および, 蒸気と雨粒の混合比は正の量なので, 計算の途中途中で ! 値が正になることを保証している. また, 元々存在する以上の雨粒が ! 蒸気に変換されないように, 元々の雨粒混合比を変換量の上限としている. ! !暗黙の型宣言禁止 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_DelMixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) !混合比の変化 real(8) :: xz_Rain2GasHeat(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xza_LatentHeat(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) real(8) :: xz_TempAll(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_ExnerAll(DimXMin:DimXMax, DimZMin:DimZMax) integer :: s !温度, 圧力, 混合比の全量を求める !擾乱成分と平均成分の足し算 xz_ExnerAll = xz_Exner + xz_ExnerBasicZ xz_TempAll = ( xz_PotTemp + xz_PotTempBasicZ ) * ( xz_Exner + xz_ExnerBasicZ ) xza_LatentHeat = 0.0d0 !雨から蒸気への相変化に伴う発熱 do s = 1, CondNum xza_LatentHeat(:,:,s) = & & xz_LatentHeat( SpcWetID(IdxCR(s)), xz_TempAll ) & & * xza_DelMixRt(:,:,IdxCR(s)) & & * xz_EffMolWtBasicZ & & / (xz_ExnerAll * CpDry) end do xz_Rain2GasHeat = sum( xza_LatentHeat, 3 ) end function xz_Rain2GasHeat !!!=================================================================================!!! function xza_Cloud2Rain( xza_MixRt, DelTime ) ! ! 雲粒から雨粒への変換量を計算するためのルーチン ! 併合成長は Berry (1968) のパラメタリゼーションを利用し, ! 衝突合体成長は Kessler (1969) のパラメタリゼーションを利用する. ! ! 変換量および, 雲粒と雨粒の混合比は正の量なので, 計算の途中途中で ! 値が正になることを保証している. また, 元々存在する以上の雲粒が ! 雨粒に変換されないように, 元々の雲粒混合比を変換量の上限としている. ! 正の値を保証するために, 引数として時間刻みが必要となる. ! (AutoConv, Collect は時間刻み幅での積分値を計算) ! ! このルーチンでは, 凝縮物質と反応生成物とを区別する必要が全くないので, ! ループを回す回数を LoopNum2 回としている. ! !暗黙の型宣言禁止 implicit none !変数定義 real(8), intent(in) :: xza_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) !混合比の擾乱成分 real(8) :: DelTime !時間刻み real(8) :: xza_Cloud2Rain(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) !雲から雨への変換量 real(8) :: xza_MixRtAll(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) !混合比の擾乱成分 + 平均成分 real(8) :: xz_AutoConv(DimXMin:DimXMax, DimZMin:DimZMax) !飽和混合比 real(8) :: xz_Collect(DimXMin:DimXMax, DimZMin:DimZMax) !規格化された潜熱 ! real(8), parameter :: N0 = 5.0d7 ! real(8), parameter :: D0 = 3.66d-1 integer :: s xza_Cloud2Rain = 0.0d0 !混合比は正の値を保証 !移流拡散で負になることもあり得るので. xza_MixRtAll = max( 0.0d0, xza_MixRt + xza_MixRtBasicZ ) do s = 1, CloudNum xz_AutoConv = 0.0d0 xz_Collect = 0.0d0 !併合成長 ! Kessler (1969) のパラメタリゼーション xz_AutoConv = & & DelTime / AutoConvTime & & * max( 0.0d0, ( xza_MixRtAll(:,:,IdxC(s)) - MixRt_AutoConvCr) ) ! ! Berry (1968) のパラメタリゼーション ! xz_AutoConv = & ! & DelTime & ! & * xz_DensBasicZ & ! & * ( xza_MixRtAll(:,:,IdxC(s)) ** 3.0d0 ) * 1.0d6 & ! & / ( 60.0d0 & ! & * ( & ! & 2.0d0 * xza_MixRtAll(:,:,IdxC(s)) & ! & + 2.66d-8 * N0 / ( xz_DensBasicZ * D0 ) & ! & ) & ! & ) !衝突合体成長 ! Kessler (1969) のパラメタリゼーション xz_Collect = & & DelTime & & * 2.2d0 * FactorJ * xza_MixRtAll(:,:,IdxC(s)) & & * ( & & xza_MixRtAll(:,:,IdxR(s)) * xz_DensBasicZ & & ) ** 0.875d0 !雲の変換量: 併合成長と合体衝突の和 ! 元々の変化量を上限値として設定する. 負の値となる. xza_Cloud2Rain(:,:,IdxC(s)) = & & - min( & & xza_MixRtAll(:,:,IdxC(s)), & & ( xz_AutoConv + xz_Collect ) & & ) !雨の変換量. 符号は雲の変換量とは反対. xza_Cloud2Rain(:,:,IdxR(s)) = - xza_Cloud2Rain(:,:,IdxC(s)) end do ! write(*,*) 'C2R: ', minval(xza_Cloud2Rain(:,:,1)), maxval(xza_Cloud2Rain(:,:,1)) ! write(*,*) 'C2R: ', minval(xza_Cloud2Rain(:,:,2)), maxval(xza_Cloud2Rain(:,:,2)) ! write(*,*) 'C2R: ', minval(xza_Cloud2Rain(:,:,3)), maxval(xza_Cloud2Rain(:,:,3)) end function xza_Cloud2Rain !!!=================================================================================!!! function xza_FallRain( xza_MixRt ) ! ! 雨粒の落下による移流を求める. ! !暗黙の型宣言禁止 implicit none !変数定義 real(8), intent(in) :: xza_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) !蒸気混合比(擾乱) real(8) :: xza_MixRtAll(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) !蒸気混合比(擾乱 + 平均場) real(8) :: xza_FallRain(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) !雨粒の落下効果 real(8) :: xz_VelZRain(DimXMin:DimXMax, DimZMin:DimZMax) !雨粒落下速度 integer :: s xza_MixRtAll = max( 0.0d0, xza_MixRt + xza_MixRtBasicZ ) xza_FallRain = 0.0d0 xz_VelZRain = 0.0d0 !落下による移流 ! Dens の avr を取ってから割ると, ゼロ割が生じるので注意 do s = 1, RainNum !雨粒終端速度 xz_VelZRain = 12.2d0 * FactorJ * ( xza_MixRtAll(:,:,IdxR(s)) ** 0.125d0 ) xza_FallRain(:,:,IdxR(s)) = & & xz_avr_xr( & & xr_dz_xz(xz_DensBasicZ & & * xz_VelZRain & & * xza_MixRtAll(:,:,IdxR(s)) & & ) & & ) / xz_DensBasicZ end do call StoreMixRtRain( xza_FallRain ) end function xza_FallRain end module WarmRainPrm