!--------------------------------------------------------------------- ! 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: SvapPress, LatentHeatPerMass use MoistFunc,only: Vap2MixRt, xz_DelMixRtNH4SH !暗黙の型宣言禁止 implicit none !属性の指定 private !関数を public にする public WarmRainPrm_Init public xz_LatentHeatRain2Gas public xz_Rain2Gas public xz_Cloud2Rain public xz_Cloud2RainNH4SH 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 !!!=================================================================================!!! function xz_Rain2Gas(xz_Exner, xz_PotTemp, xz_MixRt) ! ! 雨粒から蒸気への変換量を計算するためのルーチン ! ! 変換量および, 蒸気と雨粒の混合比は正の量なので, 計算の途中途中で ! 値が正になることを保証している. また, 元々存在する以上の雨粒が ! 蒸気に変換されないように, 元々の雨粒混合比を変換量の上限としている. ! !暗黙の型宣言禁止 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) :: xz_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) !混合比の擾乱成分 real(8) :: xz_Rain2Gas(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) ! real(8) :: xz_TempAll(DimXMin:DimXMax, DimZMin:DimZMax) !温度の擾乱成分 + 平均成分 real(8) :: xz_MixRtAll(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) !混合比の擾乱成分 + 平均成分 real(8) :: NonSaturate !未飽和度(飽和混合比と蒸気の混合比の差) integer :: s, i, k !温度, 圧力, 混合比の全量を求める !擾乱成分と平均成分の足し算 xz_TempAll = ( xz_PotTemp + xz_PotTempBasicZ ) * ( xz_Exner + xz_ExnerBasicZ ) xz_MixRtAll(:,:,:) = xz_MixRt(:,:,:) + xz_MixRtBasicZ(:,:,:) xz_Rain2Gas(:,:,:) = 0.0d0 do s = 1, LoopNum do k = DimZMin, DimZmax do i = DimXMin, DimXMax !飽和蒸気圧と混合比の差(未飽和度)を計算. ! 雨から蒸気への変換量は未飽和度に応じる. NonSaturate = & & Vap2MixRt( & & MolWtWet(CloudNum(s)), & & SvapPress(SpcWetID(CloudNum(s)), xz_TempAll(i,k)), & & (xz_Exner(i,k) + xz_ExnerBasicZ(i,k)) & & ) & & - xz_MixRtAll(i,k, GasNum(s)) !雨の変換量 ! 元々の雨粒の混合比以上に蒸発が生じないように上限値を設定 xz_Rain2Gas(i,k,RainNum(s)) = & & - min( & & 4.85d-2 * max(0.0d0, NonSaturate ) & & * ( & & max( 0.0d0, xz_MixRtAll(i,k,RainNum(s)) ) & & * xz_DensBasicZ(i,k) & & ) ** 0.65d0 , & & max( 0.0d0, xz_MixRtAll(i,k,RainNum(s)) ) & & ) !蒸気の変換量 ! 雨粒の変換量とは符号が逆となる xz_Rain2Gas(i,k,GasNum(s)) = - xz_Rain2Gas(i,k,RainNum(s)) end do end do end do ! write(*,*) 'R2V: ', minval(xz_Rain2Gas(:,:,1)), maxval(xz_Rain2Gas(:,:,1)) ! write(*,*) 'R2V: ', minval(xz_Rain2Gas(:,:,2)), maxval(xz_Rain2Gas(:,:,2)) ! write(*,*) 'R2V: ', minval(xz_Rain2Gas(:,:,3)), maxval(xz_Rain2Gas(:,:,3)) end function xz_Rain2Gas !!!=================================================================================!!! function xz_LatentHeatRain2Gas(xz_Exner, xz_PotTemp, xz_MixRt) ! ! 雨粒から蒸気への変換量を計算するためのルーチン ! ! 変換量および, 蒸気と雨粒の混合比は正の量なので, 計算の途中途中で ! 値が正になることを保証している. また, 元々存在する以上の雨粒が ! 蒸気に変換されないように, 元々の雨粒混合比を変換量の上限としている. ! !暗黙の型宣言禁止 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) :: xz_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) !混合比の擾乱成分 real(8) :: xz_LatentHeatRain2Gas(DimXMin:DimXMax, DimZMin:DimZMax) ! real(8) :: xz_Rain2Gas(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) ! real(8) :: xz_TempAll(DimXMin:DimXMax, DimZMin:DimZMax) !温度の擾乱成分 + 平均成分 real(8) :: xz_MixRtAll(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) !混合比の擾乱成分 + 平均成分 real(8) :: NonSaturate !未飽和度(飽和混合比と蒸気の混合比の差) integer :: s, i, k !温度, 圧力, 混合比の全量を求める !擾乱成分と平均成分の足し算 xz_TempAll = ( xz_PotTemp + xz_PotTempBasicZ ) * ( xz_Exner + xz_ExnerBasicZ ) xz_MixRtAll(:,:,:) = xz_MixRt(:,:,:) + xz_MixRtBasicZ(:,:,:) xz_Rain2Gas(:,:,:) = 0.0d0 xz_LatentHeatRain2Gas(:,:) = 0.0d0 do s = 1, LoopNum do k = DimZMin, DimZmax do i = DimXMin, DimXMax !飽和蒸気圧と混合比の差(未飽和度)を計算. ! 雨から蒸気への変換量は未飽和度に応じる. NonSaturate = & & Vap2MixRt( & & MolWtWet(CloudNum(s)), & & SvapPress(SpcWetID(CloudNum(s)), xz_TempAll(i,k)), & & (xz_Exner(i,k) + xz_ExnerBasicZ(i,k)) & & ) & & - xz_MixRtAll(i,k, GasNum(s)) !雨の変換量 ! 元々の雨粒の混合比以上に蒸発が生じないように上限値を設定 xz_Rain2Gas(i,k,RainNum(s)) = & & - min( & & 4.85d-2 * max(0.0d0, NonSaturate ) & & * ( & & max( 0.0d0, xz_MixRtAll(i,k,RainNum(s)) ) & & * xz_DensBasicZ(i,k) & & ) ** 0.65d0 , & & max( 0.0d0, xz_MixRtAll(i,k,RainNum(s)) ) & & ) !雨から蒸気への相変化に伴う発熱 xz_LatentHeatRain2Gas(i,k) = & & xz_LatentHeatRain2Gas(i,k) & & + LatentHeatPerMass( SpcWetID(RainNum(s)), xz_TempAll(i,k) ) & & * xz_Rain2Gas(i,k,RainNum(s)) end do end do end do end function xz_LatentHeatRain2Gas !!!=================================================================================!!! function xz_Cloud2Rain( xz_MixRt ) ! ! 雲粒から雨粒への変換量を計算するためのルーチン ! 併合成長は Kessler (1969) のパラメタリゼーションを利用し, ! 衝突合体成長は Kessler (1969) のパラメタリゼーションを利用する. ! ! 変換量および, 雲粒と雨粒の混合比は正の量なので, 計算の途中途中で ! 値が正になることを保証している. また, 元々存在する以上の雲粒が ! 雨粒に変換されないように, 元々の雲粒混合比を変換量の上限としている. ! !暗黙の型宣言禁止 implicit none !変数定義 real(8), intent(in) :: xz_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) !混合比の擾乱成分 real(8) :: xz_Cloud2Rain(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) !雲から雨への変換量 real(8) :: xz_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 :: i, k, s xz_MixRtAll = xz_MixRt + xz_MixRtBasicZ xz_Cloud2Rain = 0.0d0 do s = 1, LoopNum xz_AutoConv = 0.0d0 xz_Collect = 0.0d0 do k = DimZMin, DimZmax do i = DimXMin, DimXMax !併合成長 ! Kessler (1969) のパラメタリゼーション xz_AutoConv(i,k) = & & 1.0d-3 * max( 0.0d0, ( xz_MixRtAll(i,k,CloudNum(s)) - 1.0d-3) ) ! ! Berry (1968) のパラメタリゼーション ! xz_AutoConv(i,k) = & ! & xz_DensBasicZ(i,k) & ! & * ( & ! & max( 0.0d0, xz_MixRtAll(i,k,CloudNum(s)) ) ** 3.0d0 & ! & ) * 1.0d6 & ! & / 60.0d0 * ( & ! & 2.0d0 * max( 0.0d0, xz_MixRtAll(i,k,CloudNum(s)) ) & ! & + 60.0d0 * 2.66d-8 * N0 / ( xz_DensBasicZ(i,k) * D0 ) & ! & ) !衝突合体成長 ! Kessler (1969) のパラメタリゼーション xz_Collect(i,k) = & & 2.2d0 * max( 0.0d0, xz_MixRtAll(i,k, CloudNum(s)) ) & & * ( & & max( 0.0d0, xz_MixRtAll(i,k,RainNum(s)) ) * xz_DensBasicZ(i,k) & & ) ** 0.875d0 !雲の変換量: 併合成長と合体衝突の和 ! 元々の変化量を上限値として設定する. 負の値となる. xz_Cloud2Rain(i,k,CloudNum(s)) = & & - min( & & max( 0.0d0, xz_MixRtAll(i,k,CloudNum(s)) ), & & ( xz_AutoConv(i,k) + xz_Collect(i,k) ) & & ) !雨の変換量. 符号は雲の変換量とは反対. xz_Cloud2Rain(i,k,RainNum(s)) = - xz_Cloud2Rain(i,k,CloudNum(s)) end do end do end do ! write(*,*) 'C2R: ', minval(xz_Cloud2Rain(:,:,1)), maxval(xz_Cloud2Rain(:,:,1)) ! write(*,*) 'C2R: ', minval(xz_Cloud2Rain(:,:,2)), maxval(xz_Cloud2Rain(:,:,2)) ! write(*,*) 'C2R: ', minval(xz_Cloud2Rain(:,:,3)), maxval(xz_Cloud2Rain(:,:,3)) end function xz_Cloud2Rain !!!=================================================================================!!! function xz_Cloud2RainNH4SH( xz_MixRt ) ! ! 雲粒から雨粒への変換量を計算するためのルーチン ! 併合成長は Kessler (1969) のパラメタリゼーションを利用し, ! 衝突合体成長は Kessler (1969) のパラメタリゼーションを利用する. ! ! 変換量および, 雲粒と雨粒の混合比は正の量なので, 計算の途中途中で ! 値が正になることを保証している. また, 元々存在する以上の雲粒が ! 雨粒に変換されないように, 元々の雲粒混合比を変換量の上限としている. ! !暗黙の型宣言禁止 implicit none !変数定義 real(8), intent(in) :: xz_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) !混合比の擾乱成分 real(8) :: xz_Cloud2RainNH4SH(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) !雲から雨への変換量 real(8) :: xz_MixRtAll(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) !混合比の擾乱成分 + 平均成分 real(8) :: xz_AutoConv(DimXMin:DimXMax, DimZMin:DimZMax) !飽和混合比 real(8) :: xz_Collect(DimXMin:DimXMax, DimZMin:DimZMax) !規格化された潜熱 integer :: i, k xz_MixRtAll = xz_MixRt + xz_MixRtBasicZ xz_Cloud2RainNH4SH = 0.0d0 xz_AutoConv = 0.0d0 xz_Collect = 0.0d0 do k = DimZMin, DimZmax do i = DimXMin, DimXMax !併合成長 ! Kessler (1969) のパラメタリゼーション xz_AutoConv(i,k) = & & 1.0d-3 * max( 0.0d0, ( xz_MixRtAll(i,k,NH4SHCloudNum) - 1.0d-3) ) !衝突合体成長 ! Kessler (1969) のパラメタリゼーション xz_Collect(i,k) = & & 2.2d0 * max( 0.0d0, xz_MixRtAll(i,k, NH4SHCloudNum) ) & & * ( & & max( 0.0d0, xz_MixRtAll(i,k,NH4SHRainNum) ) * xz_DensBasicZ(i,k) & & ) ** 0.875d0 !雲の変換量: 併合成長と合体衝突の和 ! 元々の変化量を上限値として設定する. 負の値となる. xz_Cloud2RainNH4SH(i,k,NH4SHCloudNum) = & & - min( & & max( 0.0d0, xz_MixRtAll(i,k,NH4SHCloudNum) ),& & ( xz_AutoConv(i,k) + xz_Collect(i,k) ) & & ) !雨の変換量. 符号は雲の変換量とは反対. xz_Cloud2RainNH4SH(i,k,NH4SHRainNum) = - xz_Cloud2RainNH4SH(i,k,NH4SHCloudNum) end do end do ! write(*,*) 'C2R: ', minval(xz_Cloud2Rain(:,:,1)), maxval(xz_Cloud2Rain(:,:,1)) ! write(*,*) 'C2R: ', minval(xz_Cloud2Rain(:,:,2)), maxval(xz_Cloud2Rain(:,:,2)) ! write(*,*) 'C2R: ', minval(xz_Cloud2Rain(:,:,3)), maxval(xz_Cloud2Rain(:,:,3)) end function xz_Cloud2RainNH4SH !!!=================================================================================!!! function xz_FallRain( xz_MixRt, num ) ! ! 雨粒の落下による移流を求める. ! ! このルーチンの引数 xz_MixRt は 2 次元配列である. ! 引数に与えられた混合比に対し, 移流を計算する. ! もしも, xz_MixRt が雨粒の混合比を表すなら(RainSW = 1)なら ! その値を出力する. ! もしも, xz_MixRt が蒸気もしくは雲粒の混合比を表すなら(RainSW = 0) ! ならば, 落下による移流はゼロとする. ! !暗黙の型宣言禁止 implicit none !変数定義 integer, intent(in) :: num real(8), intent(in) :: xz_MixRt(DimXMin:DimXMax, DimZMin:DimZMax) !蒸気混合比(擾乱) real(8) :: xz_MixRtAll(DimXMin:DimXMax, DimZMin:DimZMax) !蒸気混合比(擾乱 + 平均場) real(8) :: xz_FallRain(DimXMin:DimXMax, DimZMin:DimZMax) !雨粒の落下効果 real(8) :: xz_VelZRain(DimXMin:DimXMax, DimZMin:DimZMax) !雨粒落下速度 xz_MixRtAll = xz_MixRt + xz_MixRtBasicZ(:,:,num) xz_FallRain = 0.0d0 xz_VelZRain = 0.0d0 where (xz_MixRt > 1.0d-16) !雨粒終端速度 xz_VelZRain = 12.2d0 * (xz_MixRtAll ** 0.125d0) !落下による移流 ! Dens の avr を取ってから割ると, ゼロ割が生じるので注意 xz_FallRain = & & - xz_avr_xr( & & xr_dz_xz(xz_DensBasicZ * xz_VelZRain * xz_MixRtAll) & & ) / xz_DensBasicZ & & * RainSW(num) end where ! write(*,*) 'MixRt: ', minval( xz_MixRt ), maxval( xz_MixRt ) ! write(*,*) 'Fall: ', minval( xz_FallRain ), maxval( xz_FallRain ) end function xz_FallRain end module WarmRainPrm