!--------------------------------------------------------------------- ! Copyright (C) GFD Dennou Club, 2005, 2006. All rights reserved. !--------------------------------------------------------------------- != Subroutine WarmRainPrm ! ! * Developer: SUGIYAMA Ko-ichiro ! * Version: $Id: $ ! * Tag Name: $Name: $ ! * Change History: ! !== Overview ! !暖かい雨のバルク法を用いた, 水蒸気と雨, 雲と雨の混合比の変換係数を求める. ! * 中島健介 (1994) で利用した定式をそのまま利用. ! !== Error Handling ! !== Bugs ! !== Note ! !== Future Plans ! ! module WarmRainPrm ! !暖かい雨のバルク法を用いた, 水蒸気と雨, 雲と雨の混合比の変換係数を求める. ! * 中島健介 (1994) で利用した定式をそのまま利用. ! !モジュール読み込み use gridset, only: DimXMin, &! x 方向の配列の下限 & DimXMax, &! x 方向の配列の上限 & DimZMin, &! z 方向の配列の下限 & DimZMax, &! z 方向の配列の上限 & SpcNum !化学種の数 use basicset, only: CpDry, &!乾燥成分の比熱 & MolWtWet, &! & SpcWetID, &! & SpcWetSymbol, &! & xz_DensBasicZ, &!基本場の密度 & xz_PotTempBasicZ, &!基本場の温位 & xz_ExnerBasicZ, &!基本場の無次元圧力 & xz_MixRtBasicZ !基本場の混合比 use average, only: xz_avr_xr use differentiate_center4, only: xr_dz_xz use ChemData, only: ChemData_OneSpcID use ChemCalc, only: xz_SvapPress, xz_LatentHeatPerMass, & & ReactHeatNH4SHPerMass use ConvUnit, only: xz_Vap2MixRt, xz_DelMixRtNH4SH, xz_Exner2Press !暗黙の型宣言禁止 implicit none !属性の指定 private !関数を public にする public WarmRainPrm_Init public WarmRainPrmSvapPress public WarmRainPrmNH4SH public xz_AutoCond public xz_Collect public xz_Evaporate public xz_FallRain real(8) :: MinError = 1.0d-12 !許容相対誤差 integer :: LoopNum = 0 integer :: GasNum(10) = 0 integer :: CloudNum(10) = 0 integer :: RainNum(10) = 0 integer :: NH3Num = 0 integer :: H2SNum = 0 integer :: NH4SHCloudNum = 0 integer :: NH4SHRainNum = 0 real(8), allocatable :: RainSW(:) save MinError save LoopNum, RainNum, CloudNum, GasNum save NH3Num, H2SNum, NH4SHCloudNum, NH4SHRainNum save RainSW contains subroutine WarmRainPrm_Init() !暗黙の型宣言禁止 implicit none !変数定義 integer :: s integer :: n1, n2, LoopNum2 !----------------------------------------------------------- ! 雨粒と雲粒と気体の ID の組を作る !----------------------------------------------------------- !初期化 allocate( RainSW(SpcNum) ) LoopNum = 0 LoopNum2 = 0 RainSW = 0.0d0 !化学種の中から雨粒を作るものを選び, その配列添え字と分子量を保管. SelectCloud: do s = 1, SpcNum !'Cloud' という文字列が含まれるものの個数を数える n1 = index(SpcWetSymbol(s), '-Cloud' ) if (n1 /= 0) then LoopNum = LoopNum + 1 GasNum(LoopNum) = minloc(SpcWetID, 1, SpcWetID == ChemData_OneSpcID(SpcWetSymbol(s)(1:n1-3) // '-g')) CloudNum(LoopNum) = s end if !'Rain' という文字列が含まれるものの個数を数える n2 = index(SpcWetSymbol(s), '-Rain' ) if (n2 /= 0) then LoopNum2 = LoopNum2 + 1 RainNum(LoopNum2) = s RainSW(s) = 1.0d0 end if ! NH4SH が存在する場合は LoopNum を 1 つ減らす if ( trim(SpcWetSymbol(s)) == 'NH4SH-s-Cloud' ) then LoopNum = LoopNum - 1 end if end do SelectCloud !----------------------------------------------------------- ! 硫化アンモニウム, およびアンモニアと硫化水素の ID を取得 ! 'Cloud' の方が 'Rain' よりも先に存在すると仮定している !----------------------------------------------------------- NH3Num = minloc(SpcWetID, 1, SpcWetID == ChemData_OneSpcID('NH3-g')) H2SNum = minloc(SpcWetID, 1, SpcWetID == ChemData_OneSpcID('H2S-g')) NH4SHCloudNum = minloc(SpcWetID, 1, SpcWetID == ChemData_OneSpcID('NH4SH-s')) NH4SHRainNum = NH4SHCloudNum + 1 !----------------------------------------------------------- ! 確認 !----------------------------------------------------------- if ( LoopNum == 0 ) then write(*,*) "WarmRainPrm: CloudNum = 0, please comment out of WarmRainPrm" stop end if write(*,*) "WarmRainPrm_Init, LoopNum: ", LoopNum write(*,*) "WarmRainPrm_Init, GasNum: ", GasNum write(*,*) "WarmRainPrm_Init, CloudNum: ", CloudNum write(*,*) "WarmRainPrm_Init, RainNum: ", RainNum write(*,*) "WarmRainPrm_Init, RainSW: ", RainSW write(*,*) "WarmRainPrm_Init, NH3Num: ", NH3Num write(*,*) "WarmRainPrm_Init, H2SNum: ", H2SNum write(*,*) "WarmRainPrm_Init, NH4SHNum: ", NH4SHCloudNum write(*,*) "WarmRainPrm_Init, NH4SHNum: ", NH4SHRainNum end subroutine WarmRainPrm_Init subroutine WarmRainPrmSvapPress(DelTime, xz_Exner, xz_PotTemp, xz_MixRt) !暗黙の型宣言禁止 implicit none !変数定義 real(8), intent(in) :: DelTime !時間刻み real(8), intent(in) :: xz_Exner(DimXMin:DimXMax, DimZMin:DimZMax) !無次元圧力の擾乱成分 real(8), intent(inout) :: xz_PotTemp(DimXMin:DimXMax, DimZMin:DimZMax) !温位の擾乱成分 real(8), intent(inout) :: xz_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) !混合比の擾乱成分 real(8) :: xz_PotTempAll(DimXMin:DimXMax, DimZMin:DimZMax) !温位の擾乱成分 + 平均成分 real(8) :: xz_ExnerAll(DimXMin:DimXMax, DimZMin:DimZMax) !無次元圧力の擾乱成分 + 平均成分 real(8) :: xz_TempAll(DimXMin:DimXMax, DimZMin:DimZMax) !温度の擾乱成分 + 平均成分 real(8) :: xz_MixRtAll(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) !混合比の擾乱成分 + 平均成分 real(8) :: xz_CNcr(DimXMin:DimXMax, DimZMin:DimZMax) ! real(8) :: xz_CLcr(DimXMin:DimXMax, DimZMin:DimZMax) ! real(8) :: xz_EVrv(DimXMin:DimXMax, DimZMin:DimZMax) ! real(8) :: xz_MixRtSat(DimXMin:DimXMax, DimZMin:DimZMax) !飽和混合比 real(8) :: xz_Gamma(DimXMin:DimXMax, DimZMin:DimZMax) !規格化された潜熱 integer :: s !温度, 圧力, 混合比の全量を求める !擾乱成分と平均成分の足し算 xz_PotTempAll = xz_PotTemp + xz_PotTempBasicZ xz_ExnerAll = xz_Exner + xz_ExnerBasicZ xz_TempAll = xz_ExnerAll * xz_PotTempAll xz_MixRtAll = xz_MixRt + xz_MixRtBasicZ do s = 1, LoopNum xz_CNcr = 0.0d0 xz_CLcr = 0.0d0 xz_EVrv = 0.0d0 write(*,*) 'pre: ', minval(xz_MixRt(:,:,GasNum(s))), maxval(xz_MixRt(:,:,GasNum(s))) write(*,*) 'pre: ', minval(xz_MixRt(:,:,CloudNum(s))), maxval(xz_MixRt(:,:,CloudNum(s))) write(*,*) 'pre: ', minval(xz_MixRt(:,:,RainNum(s))), maxval(xz_MixRt(:,:,RainNum(s))) !飽和蒸気圧 xz_MixRtSat = xz_Vap2MixRt(MolWtWet(CloudNum(s)), & & xz_SvapPress(SpcWetID(CloudNum(s)), xz_TempAll),& & xz_ExnerAll) write(*,*) 'sat: ', maxval(xz_SvapPress(SpcWetID(CloudNum(s)), xz_TempAll)) write(*,*) 'sat: ', maxval(xz_MixRtSat ) write(*,*) 'sat: ', maxval(xz_ExnerAll ) write(*,*) 'sat: ', maxval(xz_TempAll ) write(*,*) 'sat: ', SpcWetID(CloudNum(s)), MolWtWet(CloudNum(s)) !カテゴリ間の変換係数を求める xz_CNcr = xz_AutoCond(xz_MixRtAll(:,:,CloudNum(s))) xz_CLcr = xz_Collect(xz_MixRtAll(:,:,CloudNum(s)), xz_MixRtAll(:,:,RainNum(s))) xz_EVrv = xz_Evaporate(xz_MixRtSat(:,:), xz_MixRtAll(:,:,GasNum(s)), & & xz_MixRtAll(:,:,RainNum(s))) write(*,*) 'CN ', minval(xz_CNcr), maxval(xz_CNcr) write(*,*) 'CL ', minval(xz_CLcr), maxval(xz_CLcr) write(*,*) 'EV ', minval(xz_EVrv), maxval(xz_EVrv) !規格化された潜熱 xz_Gamma = xz_LatentHeatPerMass(SpcWetID(CloudNum(s)), xz_TempAll) & & / (xz_ExnerAll * CpDry) !温位 xz_PotTemp = & & xz_PotTemp & & - DelTime * xz_Gamma * xz_EVrv !水蒸気混合比 xz_MixRt(:,:,GasNum(s)) = & & xz_MixRt(:,:,GasNum(s)) & & + DelTime & & * ( & & + xz_EVrv & & ) !雲粒混合比 xz_MixRt(:,:,CloudNum(s)) = & & xz_MixRt(:,:,CloudNum(s)) & & + DelTime & & * ( & & - xz_CNcr & & - xz_CLcr & & ) !雨粒混合比 xz_MixRt(:,:,RainNum(s)) = & & xz_MixRt(:,:,RainNum(s)) & & + DelTime & & * ( & & + xz_CNcr & & + xz_CLcr & & - xz_EVrv & & ) write(*,*) 'nxt: ', minval(xz_MixRt(:,:,GasNum(s))), maxval(xz_MixRt(:,:,GasNum(s))) write(*,*) 'nxt: ', minval(xz_MixRt(:,:,CloudNum(s))), maxval(xz_MixRt(:,:,CloudNum(s))) write(*,*) 'nxt: ', minval(xz_MixRt(:,:,RainNum(s))), maxval(xz_MixRt(:,:,RainNum(s))) end do end subroutine WarmRainPrmSvapPress subroutine WarmRainPrmNH4SH(DelTime, xz_Exner, xz_PotTemp, xz_MixRt) ! !このルーチンは未完成 ! !暗黙の型宣言禁止 implicit none !変数定義 real(8), intent(in) :: DelTime !時間刻み real(8), intent(in) :: xz_Exner(DimXMin:DimXMax, DimZMin:DimZMax) !無次元圧力の擾乱成分 real(8), intent(inout) :: xz_PotTemp(DimXMin:DimXMax, DimZMin:DimZMax) !温位の擾乱成分 real(8), intent(inout) :: xz_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) !混合比の擾乱成分 real(8) :: xz_PotTempAll(DimXMin:DimXMax, DimZMin:DimZMax) !温位の擾乱成分 + 平均成分 real(8) :: xz_ExnerAll(DimXMin:DimXMax, DimZMin:DimZMax) !無次元圧力の擾乱成分 + 平均成分 real(8) :: xz_TempAll(DimXMin:DimXMax, DimZMin:DimZMax) !温度の擾乱成分 + 平均成分 real(8) :: xz_MixRtAll(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) !混合比の擾乱成分 + 平均成分 real(8) :: xz_DelMixRt(DimXMin:DimXMax, DimZMin:DimZMax) !混合比の擾乱成分 + 平均成分 real(8) :: xz_CNcr(DimXMin:DimXMax, DimZMin:DimZMax) ! real(8) :: xz_CLcr(DimXMin:DimXMax, DimZMin:DimZMax) ! real(8) :: xz_EVrv(DimXMin:DimXMax, DimZMin:DimZMax) ! real(8) :: xz_MixRtSat(DimXMin:DimXMax, DimZMin:DimZMax) !飽和混合比 real(8) :: xz_Gamma(DimXMin:DimXMax, DimZMin:DimZMax) !規格化された潜熱 !温度, 圧力, 混合比の全量を求める !擾乱成分と平均成分の足し算 xz_PotTempAll = xz_PotTemp + xz_PotTempBasicZ xz_ExnerAll = xz_Exner + xz_ExnerBasicZ xz_TempAll = xz_ExnerAll * xz_PotTempAll xz_MixRtAll = xz_MixRt + xz_MixRtBasicZ !NH4SH の生成量 xz_DelMixRt = xz_DelMixRtNH4SH(xz_TempAll, & & xz_Exner2Press(xz_ExnerAll), & & xz_MixRtAll(:,:,NH3Num), & & xz_MixRtAll(:,:,H2SNum), & & xz_MixRtAll(:,:,NH4SHCloudNum), & & MolWtWet(NH3Num), & & MolWtWet(H2SNum), & & MolWtWet(NH4SHCloudNum) ) !飽和蒸気圧混合比 ! 雨から気相への変換には, 量の多い NH3 の平衡時の混合比を用いる. xz_MixRtSat = xz_MixRtAll(:,:,NH3Num) & & + xz_DelMixRt * MolWtWet(NH3Num) / MolWtWet(NH4SHCloudNum) !カテゴリ間の変換係数を求める xz_CNcr = xz_AutoCond(xz_MixRtAll(:,:,NH4SHCloudNum)) xz_CLcr = xz_Collect(xz_MixRtAll(:,:,NH4SHCloudNum), xz_MixRtAll(:,:,NH4SHRainNum)) xz_EVrv = xz_Evaporate(xz_MixRtSat(:,:), xz_MixRtAll(:,:,NH3Num), & & xz_MixRtAll(:,:,NH4SHRainNum)) !規格化された反応熱 (NH4SH １kg に対する熱量) xz_Gamma = ReactHeatNH4SHPerMass / (xz_ExnerAll * CpDry) !温位 xz_PotTemp = & & xz_PotTemp & & - DelTime * xz_Gamma * xz_EVrv !水蒸気混合比 xz_MixRt(:,:,NH3Num) = & & xz_MixRt(:,:,NH3Num) & & + DelTime & & * ( & & + xz_EVrv * MolWtWet(NH3Num) / MolWtWet(NH4SHCloudNum) & & ) !水蒸気混合比 xz_MixRt(:,:,H2SNum) = & & xz_MixRt(:,:,H2SNum) & & + DelTime & & * ( & & + xz_EVrv * MolWtWet(H2SNum) / MolWtWet(NH4SHCloudNum) & & ) !雲粒混合比 xz_MixRt(:,:,NH4SHCloudNum) = & & xz_MixRt(:,:,NH4SHCloudNum) & & + DelTime & & * ( & & - xz_CNcr & & - xz_CLcr & & ) !雨粒混合比 xz_MixRt(:,:,NH4SHRainNum) = & & xz_MixRt(:,:,NH4SHRainNum) & & + DelTime & & * ( & & + xz_CNcr & & + xz_CLcr & & - xz_EVrv & & ) end subroutine WarmRainPrmNH4SH function xz_AutoCond(xz_MixRtC) ! ! 併合成長による雲水から雨水への変換 ! Berry (1968) のパラメタリゼーション ! !暗黙の型宣言禁止 implicit none !変数定義 real(8), intent(in) :: xz_MixRtC(DimXMin:DimXMax, DimZMin:DimZMax) !雲粒混合比 real(8) :: xz_AutoCond(DimXMin:DimXMax, DimZMin:DimZMax) !併合成長 real(8), parameter :: N0 = 5.0d7 real(8), parameter :: D0 = 3.66d-1 ! real(8), parameter :: K1 = 1.0d-3 !ケスラーの係数 ! real(8), parameter :: A = 1.0d-3 !雲水量の閾値 integer :: i, k do k = DimZMin, DimZmax do i = DimXMin, DimXMax xz_AutoCond(i,k) = & & xz_DensBasicZ(i,k) & & * (max(0.0d0, xz_MixRtC(i,k)) ** 3.0d0) * 1.0d6 & & / (60.0d0 * 2.0d0 * max(0.0d0, xz_MixRtC(i,k) ) & & + 60.0d0 * 2.66d-8 * N0 / (xz_DensBasicZ(i,k) * D0)) end do end do ! ! Kessler (1969) のパラメタリゼーション ! xz_AutoCond = K1 * (xz_MixRtC - A) end function xz_AutoCond function xz_Collect(xz_MixRtC, xz_MixRtR) ! ! 凝結による水蒸気から雲水への変換 ! Kessler (1969) のパラメタリゼーション ! !暗黙の型宣言禁止 implicit none !変数定義 real(8), intent(in) :: xz_MixRtC(DimXMin:DimXMax, DimZMin:DimZMax) !雲粒混合比 real(8), intent(in) :: xz_MixRtR(DimXMin:DimXMax, DimZMin:DimZMax) !雨粒混合比 real(8) :: xz_Collect(DimXMin:DimXMax, DimZMin:DimZMax) !衝突合体成長 integer :: i, k do k = DimZMin, DimZmax do i = DimXMin, DimXMax xz_Collect(i,k) = & & 2.2d0 * max(0.0d0, xz_MixRtC(i,k)) & & * ( max(0.0d0, xz_MixRtR(i,k)) * xz_DensBasicZ(i,k) )** 0.875d0 end do end do end function xz_Collect function xz_Evaporate(xz_MixRtSat, xz_MixRtV, xz_MixRtR) ! ! 蒸発による雨水から水蒸気への変換 ! !暗黙の型宣言禁止 implicit none !変数定義 real(8), intent(in) :: xz_MixRtSat(DimXMin:DimXMax, DimZMin:DimZMax) !蒸気混合比 real(8), intent(in) :: xz_MixRtV(DimXMin:DimXMax, DimZMin:DimZMax) !蒸気混合比(擾乱+基本場) real(8), intent(in) :: xz_MixRtR(DimXMin:DimXMax, DimZMin:DimZMax) !蒸気混合比 real(8) :: xz_Evaporate(DimXMin:DimXMax, DimZMin:DimZMax) !雨粒から水蒸気への変換 integer :: i, k do k = DimZMin, DimZmax do i = DimXMin, DimXMax xz_Evaporate(i,k) = & & 4.85d-2 * max(0.0d0, (xz_MixRtSat(i,k) - xz_MixRtV(i,k))) & & * ( max(0.0d0, xz_MixRtR(i,k)) * xz_DensBasicZ(i,k) ) ** 0.65d0 end do end do end function xz_Evaporate function xz_FallRain(Num, xz_MixRtR) ! !凝結による水蒸気から雲水への変換. 終端速度を計算してから評価. ! !暗黙の型宣言禁止 implicit none !変数定義 integer, intent(in) :: Num !配列添え字 real(8), intent(in) :: xz_MixRtR(DimXMin:DimXMax, DimZMin:DimZMax) !蒸気混合比 real(8) :: xz_FallRain(DimXMin:DimXMax, DimZMin:DimZMax) !雨粒の落下効果 real(8) :: xz_VelZRain(DimXMin:DimXMax, DimZMin:DimZMax) !雨粒落下速度 xz_VelZRain = 0.0d0 xz_FallRain = 0.0d0 where (xz_MixRtR > 1.0d-16) !雨粒終端速度 xz_VelZRain = 12.2d0 * (xz_MixRtR ** 0.125d0) !凝結による水蒸気から雲水への変換. !Dens の avr を取ってから割ると, ゼロ割が生じるので注意 xz_FallRain = & & xz_avr_xr( & & xr_dz_xz(xz_DensBasicZ * xz_VelZRain * xz_MixRtR) & & ) / xz_DensBasicZ & & * RainSW(Num) end where write(*,*) 'FallRain: ', minval(xz_FallRain), maxval(xz_FallRain) end function xz_FallRain end module WarmRainPrm