!= Module WarmRainPrm_3d ! ! Authors:: SUGIYAMA Ko-ichiro, ODAKA Masatsugu ! Version:: $Id: warmrainprm_3d.f90,v 1.1 2008-06-19 16:53:25 odakker Exp $ ! Tag Name:: $Name: arare4-20100306 $ ! 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_3d ! !暖かい雨のバルク法を用いた, 水蒸気と雨, 雲と雨の混合比の変換係数を求める. ! * 中島健介 (1994) で利用した定式をそのまま利用. ! !モジュール読み込み use dc_types, only : DP use dc_iounit, only : FileOpen use dc_message, only : MessageNotify use gridset_3d, only : DimXMin, &!x 方向の配列の下限 & DimXMax, &!x 方向の配列の上限 & DimYMin, &!y 方向の配列の上限 & DimYMax, &!y 方向の配列の上限 & DimZMin, &!z 方向の配列の下限 & DimZMax, &!z 方向の配列の上限 & SpcNum !化学種の数 use basicset_3d,only : PressBasis, &!温位の基準圧力 & CpDry, &!乾燥成分の比熱 & MolWtWet, &! & MolWtDry, &! & SpcWetID, &! & xyz_DensBasicZ, &!基本場の密度 & xyz_PotTempBasicZ, &!基本場の温位 & xyz_ExnerBasicZ, &!基本場の無次元圧力 & xyz_EffMolWtBasicZ,&!基本場の分子量の寄与 & xyza_MixRtBasicZ, &!基本場の混合比 & GasRDry !乾燥成分の気体定数 use moistset, only : CondNum, &!凝結過程の数 & IdxCG, &!凝結過程(蒸気)の配列添え字 & IdxCC, &!凝結過程(雲)の配列添え字 & IdxCR, &!凝結過程(雲)の配列添え字 & CloudNum, &!雲の数 & RainNum, &!雨の数 & IdxC, &!雲の配列添え字 & IdxR, &!雨の配列添え字 & IdxNH3, &!NH3(蒸気)の配列添え字 & IdxH2S, &!H2S(蒸気)の配列添え字 & IdxNH4SHr !NH4SH(雨)の配列添え字 use xyz_base_module, only : xyz_avr_xyr use xyz_deriv_module,only : xyr_dz_xyz use ChemCalc_3d, only : xyz_SvapPress, xyz_LatentHeat, ReactHeatNH4SH use MoistFunc_3d,only : xyz_DelMixRtNH4SH use StoreMixRt, only : StoreMixRtRain !暗黙の型宣言禁止 implicit none !属性の指定 private !関数を public にする public WarmRainPrm_Init public xyz_Rain2GasHeat public xyz_Rain2GasHeatNH4SH public xyza_Rain2Gas public xyza_Rain2GasNH4SH public xyza_Cloud2Rain public xyza_FallRain real(DP) :: FactorJ = 1.0d0 !雲物理過程のパラメータ !木星では 3.0d0 !地球では 1.0d0 とする real(DP) :: AutoConvTime = 1.0d3 !併合成長の時定数 [sec] real(DP) :: MixRt_AutoConvCr = 1.0d-3 !併合成長を生じる臨界混合比 [kg/kg] save FactorJ, AutoConvTime, MixRt_AutoConvCr contains !!!=================================================================================!!! subroutine WarmRainPrm_Init( cfgfile ) !暗黙の型宣言禁止 implicit none !入力変数 character(*), intent(in) :: cfgfile !内部変数 integer :: unit !装置番号 !----------------------------------------------------------- ! NAMELIST から情報を取得 !----------------------------------------------------------- ! NAMELIST から情報を取得 NAMELIST /warmrainprm/ & & FactorJ, AutoConvTime, MixRt_AutoConvCr call FileOpen(unit, file=cfgfile, mode='r') read(unit, NML=warmrainprm) close(unit) 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 xyza_Rain2Gas(xyz_Exner, xyz_PotTemp, xyza_MixRt, DelTime) ! ! 雨粒から蒸気への変換量を計算するためのルーチン ! ! 変換量および, 蒸気と雨粒の混合比は正の量なので, 計算の途中途中で ! 値が正になることを保証している. また, 元々存在する以上の雨粒が ! 蒸気に変換されないように, 元々の雨粒混合比を変換量の上限としている. ! ! 木星の場合は, FactorJ で変換量を加速する. ! !暗黙の型宣言禁止 implicit none !変数定義 real(DP), intent(in) :: xyz_PotTemp(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !温位の擾乱成分 real(DP), intent(in) :: xyz_Exner(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !温度の擾乱成分 real(DP), intent(in) :: xyza_MixRt(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum) !混合比の擾乱成分 real(DP) :: DelTime !時間刻み real(DP) :: xyza_Rain2Gas(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum) ! real(DP) :: xyz_TempAll(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !温度の擾乱成分 + 平均成分 real(DP) :: xyz_PressAll(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !全圧 real(DP) :: xyza_MixRtAll(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum) !混合比の擾乱成分 + 平均成分 real(DP) :: xyz_NonSaturate(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !未飽和度(飽和混合比と蒸気の混合比の差) integer :: s !温度, 圧力, 混合比の全量を求める !擾乱成分と平均成分の足し算 xyz_TempAll = ( xyz_PotTemp + xyz_PotTempBasicZ ) * ( xyz_Exner + xyz_ExnerBasicZ ) xyz_PressAll = PressBasis * ((xyz_Exner + xyz_ExnerBasicZ) ** (CpDry / GasRDry)) xyza_Rain2Gas = 0.0d0 !混合比は正の値であることを保証 !移流拡散計算で負になることがあり得るので. xyza_MixRtAll = max( 0.0d0, xyza_MixRt + xyza_MixRtBasicZ ) do s = 1, CondNum !飽和蒸気圧と混合比の差(飽和度)を計算. ! 雨から蒸気への変換量は飽和度に比例する. xyz_NonSaturate = & & max( & & 0.0d0, & & xyz_SvapPress(SpcWetID(IdxCC(s)), xyz_TempAll) & & * MolWtWet(IdxCG(s)) & & / ( MolWtDry * xyz_PressAll) & & - xyza_MixRtAll(:,:,:,IdxCG(s)) & & ) !雨の変換量 ! 元々の雨粒の混合比以上に蒸発が生じないように上限値を設定 xyza_Rain2Gas(:,:,:,IdxCR(s)) = & & - min( & & DelTime * 4.85d-2 * FactorJ * xyz_NonSaturate & & * ( xyza_MixRtAll(:,:,:,IdxCR(s)) * xyz_DensBasicZ )** 0.65d0,& & xyza_MixRtAll(:,:,:,IdxCR(s)) & & ) !蒸気の変換量 ! 雨粒の変換量とは符号が逆となる xyza_Rain2Gas(:,:,:,IdxCG(s)) = - xyza_Rain2Gas(:,:,:,IdxCR(s)) end do end function xyza_Rain2Gas !!!=================================================================================!!! function xyza_Rain2GasNH4SH(xyz_Exner, xyz_PotTemp, xyza_MixRt, DelTime) ! ! 雨粒から蒸気への変換量を計算するためのルーチン ! ! 変換量および, 蒸気と雨粒の混合比は正の量なので, 計算の途中途中で ! 値が正になることを保証している. また, 元々存在する以上の雨粒が ! 蒸気に変換されないように, 元々の雨粒混合比を変換量の上限としている. ! !暗黙の型宣言禁止 implicit none !変数定義 real(DP), intent(in) :: xyz_PotTemp(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !温位の擾乱成分 real(DP), intent(in) :: xyz_Exner(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !温度の擾乱成分 real(DP), intent(in) :: xyza_MixRt(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum) !混合比の擾乱成分 real(DP), intent(in) :: DelTime !時間刻み real(DP) :: xyza_Rain2GasNH4SH(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum) ! real(DP) :: xyz_TempAll(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !温度の擾乱成分 + 平均成分 real(DP) :: xyz_PressAll(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !圧力の擾乱成分 + 平均成分 real(DP) :: xyza_MixRtAll(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum) !混合比の擾乱成分 + 平均成分 real(DP) :: xyz_NonSaturate(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !未飽和度(飽和混合比と蒸気の混合比の差) !温度, 圧力, 混合比の全量を求める !擾乱成分と平均成分の足し算 xyz_TempAll = ( xyz_PotTemp + xyz_PotTempBasicZ ) * ( xyz_Exner + xyz_ExnerBasicZ ) xyz_PressAll= PressBasis * ((xyz_Exner + xyz_ExnerBasicZ) ** (CpDry / GasRDry)) xyza_Rain2GasNH4SH = 0.0d0 !混合比は正の値であることを保証 !移流拡散計算で負になることがあり得るので. xyza_MixRtAll = max( 0.0d0, xyza_MixRt + xyza_MixRtBasicZ ) !飽和蒸気圧と混合比の差(飽和度)を計算. ! 雨から蒸気への変換量は飽和度に比例する. ! 未飽和度を求めたいので, マイナスをかけ算している ! (DelMixRtNH4SH は, NH4SH が増加する方向, すなわち飽和度を正としている) xyz_NonSaturate = & & max( & & 0.0d0, & & - xyz_DelMixRtNH4SH( & & xyz_TempAll, xyz_PressAll, & & xyza_MixRtAll(:,:,:,IdxNH3), & & xyza_MixRtAll(:,:,:,IdxH2S), & & MolWtWet(IdxNH3), & & MolWtWet(IdxH2S) & & ) & & ) !雨の変換量 ! 元々の雨粒の混合比以上に蒸発が生じないように上限値を設定 xyza_Rain2GasNH4SH(:,:,:,IdxNH4SHr) = & & - min( & & DelTime * 4.85d-2 * FactorJ * xyz_NonSaturate & & * ( & & xyza_MixRtAll(:,:,:,IdxNH4SHr) * xyz_DensBasicZ & & ) ** 0.65d0, & & xyza_MixRtAll(:,:,:,IdxNH4SHr) & & ) !蒸気の変換量 ! 雨粒の変換量とは符号が逆となる xyza_Rain2GasNH4SH(:,:,:,IdxNH3) = & & - xyza_Rain2GasNH4SH(:,:,:,IdxNH4SHr) * MolWtWet(IdxNH3) & & / MolWtWet(IdxNH4SHr) xyza_Rain2GasNH4SH(:,:,:,IdxH2S) = & & - xyza_Rain2GasNH4SH(:,:,:,IdxNH4SHr) * MolWtWet(IdxH2S) & & / MolWtWet(IdxNH4SHr) end function xyza_Rain2GasNH4SH !!!=================================================================================!!! function xyz_Rain2GasHeatNH4SH(xyz_Exner, xyza_DelMixRt) ! ! 雨粒から蒸気への変換量を計算するためのルーチン ! ! 変換量および, 蒸気と雨粒の混合比は正の量なので, 計算の途中途中で ! 値が正になることを保証している. また, 元々存在する以上の雨粒が ! 蒸気に変換されないように, 元々の雨粒混合比を変換量の上限としている. ! !暗黙の型宣言禁止 implicit none !変数定義 real(DP), intent(in) :: xyza_DelMixRt(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum) !混合比の変化量 real(DP), intent(in) :: xyz_Exner(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !温度の擾乱成分 real(DP) :: xyz_Rain2GasHeatNH4SH(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) ! real(DP) :: xyz_ExnerAll(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !圧力の全量を求める !擾乱成分と平均成分の足し算 xyz_ExnerAll = xyz_Exner + xyz_ExnerBasicZ !雨から蒸気への相変化に伴う発熱 xyz_Rain2GasHeatNH4SH = & & ReactHeatNH4SH * xyza_DelMixRt(:,:,:,IdxNH4SHr) * xyz_EffMolWtBasicZ & & / (xyz_ExnerAll * CpDry) end function xyz_Rain2GasHeatNH4SH !!!=================================================================================!!! function xyz_Rain2GasHeat(xyz_PotTemp, xyz_Exner, xyza_DelMixRt) ! ! 雨粒から蒸気への変換量を計算するためのルーチン ! ! 変換量および, 蒸気と雨粒の混合比は正の量なので, 計算の途中途中で ! 値が正になることを保証している. また, 元々存在する以上の雨粒が ! 蒸気に変換されないように, 元々の雨粒混合比を変換量の上限としている. ! !暗黙の型宣言禁止 implicit none !変数定義 real(DP), intent(in) :: xyz_PotTemp(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !温位の擾乱成分 real(DP), intent(in) :: xyz_Exner(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !温度の擾乱成分 real(DP), intent(in) :: xyza_DelMixRt(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum) !混合比の変化 real(DP) :: xyz_Rain2GasHeat(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) real(DP) :: xyza_LatentHeat(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum) real(DP) :: xyz_TempAll(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) real(DP) :: xyz_ExnerAll(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) integer :: s !温度, 圧力, 混合比の全量を求める !擾乱成分と平均成分の足し算 xyz_ExnerAll = xyz_Exner + xyz_ExnerBasicZ xyz_TempAll = ( xyz_PotTemp + xyz_PotTempBasicZ ) * ( xyz_Exner + xyz_ExnerBasicZ ) xyza_LatentHeat = 0.0d0 !雨から蒸気への相変化に伴う発熱 do s = 1, CondNum xyza_LatentHeat(:,:,:,s) = & & xyz_LatentHeat( SpcWetID(IdxCR(s)), xyz_TempAll ) & & * xyza_DelMixRt(:,:,:,IdxCR(s)) & & * xyz_EffMolWtBasicZ & & / (xyz_ExnerAll * CpDry) end do xyz_Rain2GasHeat = sum( xyza_LatentHeat, 3 ) end function xyz_Rain2GasHeat !!!=================================================================================!!! function xyza_Cloud2Rain( xyza_MixRt, DelTime ) ! ! 雲粒から雨粒への変換量を計算するためのルーチン ! 併合成長は Berry (1968) のパラメタリゼーションを利用し, ! 衝突合体成長は Kessler (1969) のパラメタリゼーションを利用する. ! ! 変換量および, 雲粒と雨粒の混合比は正の量なので, 計算の途中途中で ! 値が正になることを保証している. また, 元々存在する以上の雲粒が ! 雨粒に変換されないように, 元々の雲粒混合比を変換量の上限としている. ! 正の値を保証するために, 引数として時間刻みが必要となる. ! (AutoConv, Collect は時間刻み幅での積分値を計算) ! ! このルーチンでは, 凝縮物質と反応生成物とを区別する必要が全くないので, ! ループを回す回数を LoopNum2 回としている. ! !暗黙の型宣言禁止 implicit none !変数定義 real(DP), intent(in) :: xyza_MixRt(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum) !混合比の擾乱成分 real(DP) :: DelTime !時間刻み real(DP) :: xyza_Cloud2Rain(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum) !雲から雨への変換量 real(DP) :: xyza_MixRtAll(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum) !混合比の擾乱成分 + 平均成分 real(DP) :: xyz_AutoConv(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !飽和混合比 real(DP) :: xyz_Collect(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !規格化された潜熱 ! real(DP), parameter :: N0 = 5.0d7 ! real(DP), parameter :: D0 = 3.66d-1 integer :: s xyza_Cloud2Rain = 0.0d0 !混合比は正の値を保証 !移流拡散で負になることもあり得るので. xyza_MixRtAll = max( 0.0d0, xyza_MixRt + xyza_MixRtBasicZ ) do s = 1, CloudNum xyz_AutoConv = 0.0d0 xyz_Collect = 0.0d0 !併合成長 ! Kessler (1969) のパラメタリゼーション xyz_AutoConv = & & DelTime / AutoConvTime & & * max( 0.0d0, ( xyza_MixRtAll(:,:,:,IdxC(s)) - MixRt_AutoConvCr) ) ! ! Berry (1968) のパラメタリゼーション ! xyz_AutoConv = & ! & DelTime & ! & * xyz_DensBasicZ & ! & * ( xyza_MixRtAll(:,:,:,IdxC(s)) ** 3.0d0 ) * 1.0d6 & ! & / ( 60.0d0 & ! & * ( & ! & 2.0d0 * xyza_MixRtAll(:,,:,:,IdxC(s)) & ! & + 2.66d-8 * N0 / ( xyz_DensBasicZ * D0 ) & ! & ) & ! & ) !衝突合体成長 ! Kessler (1969) のパラメタリゼーション xyz_Collect = & & DelTime & & * 2.2d0 * FactorJ * xyza_MixRtAll(:,:,:,IdxC(s)) & & * ( & & xyza_MixRtAll(:,:,:,IdxR(s)) * xyz_DensBasicZ & & ) ** 0.875d0 !雲の変換量: 併合成長と合体衝突の和 ! 元々の変化量を上限値として設定する. 負の値となる. xyza_Cloud2Rain(:,:,:,IdxC(s)) = & & - min( & & xyza_MixRtAll(:,:,:,IdxC(s)), & & ( xyz_AutoConv + xyz_Collect ) & & ) !雨の変換量. 符号は雲の変換量とは反対. xyza_Cloud2Rain(:,:,:,IdxR(s)) = - xyza_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 xyza_Cloud2Rain !!!=================================================================================!!! function xyza_FallRain( xyza_MixRt ) ! ! 雨粒の落下による移流を求める. ! !暗黙の型宣言禁止 implicit none !変数定義 real(DP), intent(in) :: xyza_MixRt(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum) !蒸気混合比(擾乱) real(DP) :: xyza_MixRtAll(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum) !蒸気混合比(擾乱 + 平均場) real(DP) :: xyza_FallRain(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum) !雨粒の落下効果 real(DP) :: xyz_VelZRain(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !雨粒落下速度 integer :: s xyza_MixRtAll = max( 0.0d0, xyza_MixRt + xyza_MixRtBasicZ ) xyza_FallRain = 0.0d0 xyz_VelZRain = 0.0d0 !落下による移流 ! Dens の avr を取ってから割ると, ゼロ割が生じるので注意 do s = 1, RainNum !雨粒終端速度 xyz_VelZRain = 12.2d0 * FactorJ * ( xyza_MixRtAll(:,:,:,IdxR(s)) ** 0.125d0 ) xyza_FallRain(:,:,:,IdxR(s)) = & & xyz_avr_xyr( & & xyr_dz_xyz(xyz_DensBasicZ & & * xyz_VelZRain & & * xyza_MixRtAll(:,:,:,IdxR(s)) & & ) & & ) / xyz_DensBasicZ end do call StoreMixRtRain( xyza_FallRain ) end function xyza_FallRain end module WarmRainPrm_3d