!= Module MoistAdjust ! ! Authors:: SUGIYAMA Ko-ichiro ! Version:: $Id: moistadjust.f90,v 1.16 2006/11/22 15:40:55 sugiyama Exp $ ! Tag Name:: $Name: arare4-20061224 $ ! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved. ! License:: See COPYRIGHT[link:../../COPYRIGHT] ! !== Overview ! !湿潤飽和調節法 ! !== Error Handling ! !== Bugs ! !== Note ! !== Future Plans ! ! Module MoistAdjust ! !湿潤飽和調節法を行うためのパッケージ型モジュール ! * 飽和蒸気圧を使う化学種は, 一元的に実行する ! * 化学反応については, それぞれの化学反応毎に実行する. !モジュール読み込み use dc_message, only: MessageNotify use gridset, only: SpcNum, &!化学種の数 & DimXMin, &! x 方向の配列の下限 & DimXMax, &! x 方向の配列の上限 & DimZMin, &! z 方向の配列の下限 & DimZMax ! z 方向の配列の上限 use basicset, only: MolWtWet, & & SpcWetID, & & SpcWetSymbol, & & xz_ExnerBasicZ, & & xz_PotTempBasicZ, & & xz_EffMolWtBasicZ,& & xza_MixRtBasicZ, & & PressBasis, &!温位の基準圧力 [Pa] & CpDry, &!乾燥成分の平均定圧比熱 [J/K kg] & MolWtDry, &!乾燥成分の平均分子量 [kg/mol] & GasRDry !乾燥成分の気体定数 [J/K kg] use average, only: xz_avr_pz, xz_avr_xr, pz_avr_xz, xr_avr_xz use ChemData, only: ChemData_OneSpcID use ChemCalc, only: xz_SvapPress, xz_LatentHeat, ReactHeatNH4SH use MoistFunc,only: xz_DMixRtSatDPotTemp, xz_DelMixRtNH4SH !暗黙の型宣言禁止 implicit none !属性の指定 private !関数の属性を public に変更 public MoistAdjust_Init !初期化ルーチン public MoistAdjustPrm !変数を出力. 解析ルーチンで利用するため public MoistAdjustSvapPress !飽和蒸気圧を用いた飽和湿潤調節(簡易版) public MoistAdjustNH4SH !化学反応の圧平衡定数を用いた飽和湿潤調節 !変数定義 integer :: LoopNum = 0 integer :: GasNum(10) = 0 integer :: CloudNum(10) = 0 integer :: NH3Num = 0 integer :: H2SNum = 0 integer :: NH4SHNum = 0 save LoopNum, CloudNum, GasNum save NH3Num, H2SNum, NH4SHNum contains !!!------------------------------------------------------------------!!! subroutine MoistAdjust_Init( ) ! !計算に利用する化学種の ID の対を取得 ! * 気相と凝縮相(雲粒)の ID の対. ! * NH4SH 生成反応を計算するための, NH4SH と NH3, H2S の ID. !暗黙の型宣言禁止 implicit none !変数定義 integer :: s integer :: n1 !----------------------------------------------------------- ! 雲粒と気体の ID の組を作る !----------------------------------------------------------- !初期化 LoopNum = 0 !化学種の中から雲粒を作るものを選び, その配列添え字と分子量を保管. SelectCloud: do s = 1, SpcNum ! NH4SH については無視 if ( trim(SpcWetSymbol(s)) == 'NH4SH-s-Cloud' ) then cycle SelectCloud end if !'Cloud' という文字列が含まれるものの個数を数える n1 = index(SpcWetSymbol(s), '-Cloud' ) if (n1 /= 0) then LoopNum = LoopNum + 1 CloudNum(LoopNum)= s GasNum(LoopNum) = minloc(SpcWetID, 1, SpcWetID == ChemData_OneSpcID(SpcWetSymbol(s)(1:n1-3) // '-g')) end if end do SelectCloud !----------------------------------------------------------- ! 硫化アンモニウム, およびアンモニアと硫化水素の ID を取得 !----------------------------------------------------------- NH3Num = minloc(SpcWetID, 1, SpcWetID == ChemData_OneSpcID('NH3-g')) H2SNum = minloc(SpcWetID, 1, SpcWetID == ChemData_OneSpcID('H2S-g')) NH4SHNum = minloc(SpcWetID, 1, SpcWetID == ChemData_OneSpcID('NH4SH-s')) !----------------------------------------------------------- ! 確認 !----------------------------------------------------------- if ( LoopNum == 0 ) then call MessageNotify( "W", & & "MoistAdjust: ", "CloudNum = 0, please comment out of MoistAdjust" ) ! stop end if call MessageNotify( "M", "MoistAdjust_Init", "LoopNum=%d", i=(/LoopNum/) ) call MessageNotify( "M", "MoistAdjust_Init", "CloudNum=%d", i=(/CloudNum/) ) call MessageNotify( "M", "MoistAdjust_Init", "GasNum=%d", i=(/GasNum/) ) call MessageNotify( "M", "MoistAdjust_Init", "NH3Num=%d", i=(/NH3Num/) ) call MessageNotify( "M", "MoistAdjust_Init", "H2SNum=%d", i=(/H2SNum/) ) call MessageNotify( "M", "MoistAdjust_Init", "NH4SHNum=%d", i=(/NH4SHNum/) ) end subroutine MoistAdjust_Init subroutine MoistAdjustPrm( M1, M2, M3, M4, M5, M6 ) implicit none integer, intent(out) :: M1, M2(10), M3(10), M4, M5, M6 M1 = LoopNum M2 = CloudNum M3 = GasNum M4 = NH3Num M5 = H2SNum M6 = NH4SHNum ! write(*,*) M1, M2, M3, M4, M5, M6 end subroutine MoistAdjustPrm !!!------------------------------------------------------------------!!! subroutine MoistAdjustSvapPress(xz_Exner, xz_PotTemp, xza_MixRt) ! ! 飽和蒸気圧を用いた湿潤飽和調整法の実行 ! この副プログラムでは, 予め決めた回数だけ反復改良を行う. ! !暗黙の型宣言禁止 implicit none !変数定義 real(8),intent(in) :: xz_Exner(DimXMin:DimXMax, DimZMin:DimZMax) !エクスナー関数 real(8),intent(inout):: xz_PotTemp(DimXMin:DimXMax, DimZMin:DimZMax) !温位 real(8),intent(inout):: xza_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) !混合比 integer, parameter :: ItrNum = 4 !反復改良の回数 real(8) :: xz_MixRtV_pre(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_MixRtV_nxt(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_MixRtC_pre(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_MixRtC_nxt(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_PotTemp_pre(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_PotTemp_nxt(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_ExnerAll(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_TempAll(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_DelMixRt(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_MixRtSat(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_Cond(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_Evap(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_Gamma(DimXMin:DimXMax, DimZMin:DimZMax) integer :: i, s ! 添え字 !--------------------------------------------------------------------- ! 初期化 !--------------------------------------------------------------------- xz_MixRtV_pre = 0.0d0 xz_MixRtV_nxt = 0.0d0 xz_MixRtC_pre = 0.0d0 xz_MixRtC_nxt = 0.0d0 xz_PotTemp_pre = 0.0d0 xz_PotTemp_nxt = 0.0d0 xz_ExnerAll = 0.0d0 xz_TempAll = 0.0d0 xz_DelMixRt = 0.0d0 xz_MixRtSat = 0.0d0 xz_Gamma = 0.0d0 !--------------------------------------------------------------------- ! 湿潤飽和調節法の実行 ! ループを回すのは, 雲についてだけ. !--------------------------------------------------------------------- LoopSvapPress: do s = 1, LoopNum !湿潤飽和法では圧力は変化させない. xz_ExnerAll = xz_Exner + xz_ExnerBasicZ !今までに得られた擾乱成分の値を暫定量とみなす. 添え字を追加 xz_MixRtV_pre = xza_MixRt(:,:,GasNum(s)) + xza_MixRtBasicZ(:,:,GasNum(s)) xz_MixRtC_pre = xza_MixRt(:,:,CloudNum(s)) + xza_MixRtBasicZ(:,:,CloudNum(s)) xz_PotTemp_pre = xz_PotTemp Adjusting: do i = 1, ItrNum !--------------------------------------------------------------- ! 飽和蒸気圧から飽和混合比を求める !--------------------------------------------------------------- !温度 xz_TempAll = ( xz_PotTemp_pre + xz_PotTempBasicZ ) * xz_ExnerAll !飽和蒸気圧から飽和混合比を計算(基本場からの差). xz_MixRtSat = & & xz_SvapPress(SpcWetID(CloudNum(s)), xz_TempAll) & & * MolWtWet(CloudNum(s)) & & / (MolWtDry * PressBasis * (xz_ExnerAll ** (CpDry / GasRDry)) ) ! write(*,*) "A", minval( xz_SvapPress(SpcWetID(CloudNum(s)), xz_TempAll) ), maxval( xz_SvapPress(SpcWetID(CloudNum(s)), xz_TempAll) ) ! write(*,*) "B", minval( xz_MixRtSat ), maxval( xz_MixRtSat ) !規格化された潜熱 xz_Gamma = xz_LatentHeat(SpcWetID(CloudNum(s)), xz_TempAll) & & * xz_EffMolWtBasicZ & & / (xz_ExnerAll * CpDry) ! write(*,*) "C", minval( xz_Gamma ), maxval( xz_Gamma ) !凝結量を求める. ! 凝結が生じる場合には, xz_MixRtV_pre - xz_MixRtSat は必ず正の値となる. ! 蒸発が生じる場合には, 蒸発量は - MixRtC を超えることはない. xz_DelMixRt = & & ( xz_MixRtV_pre - xz_MixRtSat ) & & / (1.0d0 + xz_Gamma * xz_DMixRtSatDPotTemp( & & SpcWetID(CloudNum(s)), MolWtWet(CloudNum(s)), xz_TempAll, xz_ExnerAll & & ) ) ! Write(*,*) "D1", & ! & minval( xz_DMixRtSatDPotTemp( & ! & SpcWetID(CloudNum(s)), MolWtWet(CloudNum(s)), xz_TempAll, xz_ExnerAll & ! & ) ), & ! & maxval( xz_DMixRtSatDPotTemp( & ! & SpcWetID(CloudNum(s)), MolWtWet(CloudNum(s)), xz_TempAll, xz_ExnerAll & ! & ) ) ! write(*,*) "D2", minval( xz_MixRtV_pre - xz_MixRtSat ), maxval( xz_MixRtV_pre - xz_MixRtSat ) ! write(*,*) "D3", minval( xz_DelMixRt ), maxval( xz_DelMixRt ) xz_Cond = max( 0.0d0, min( xz_MixRtV_pre, xz_DelMixRt ) ) xz_Evap = max( 0.0d0, min( xz_MixRtC_pre, - xz_DelMixRt ) ) ! write(*,*) "E", minval( xz_Cond ), maxval( xz_Cond ) ! write(*,*) "F", minval( xz_Evap ), maxval( xz_Evap ) !より真に近い値を計算 xz_PotTemp_nxt = xz_PotTemp_pre + xz_Gamma * ( xz_Cond - xz_Evap ) xz_MixRtV_nxt = xz_MixRtV_pre - xz_Cond + xz_Evap xz_MixRtC_nxt = xz_MixRtC_pre + xz_Cond - xz_Evap !繰り返しのための変数定義 xz_PotTemp_pre = xz_PotTemp_nxt xz_MixRtV_pre = xz_MixRtV_nxt xz_MixRtC_pre = xz_MixRtC_nxt end do Adjusting xz_PotTemp = xz_PotTemp_nxt xza_MixRt(:,:,GasNum(s)) = xz_MixRtV_nxt - xza_MixRtBasicZ(:,:,GasNum(s)) xza_MixRt(:,:,CloudNum(s)) = xz_MixRtC_nxt - xza_MixRtBasicZ(:,:,CloudNum(s)) end do LoopSvapPress end subroutine MoistAdjustSvapPress !!!--------------------------------------------------------------------------!!! subroutine MoistAdjustNH4SH(xz_Exner, xz_PotTemp, xza_MixRt ) ! ! NH3 + H2S --> NH4SH の生成反応の圧平衡定数 Kp を用いた飽和湿潤調節法 ! !暗黙の型宣言禁止 implicit none !変数定義 real(8),intent(in) :: xz_Exner(DimXMin:DimXMax, DimZMin:DimZMax) !エクスナー関数 real(8),intent(inout):: xz_PotTemp(DimXMin:DimXMax, DimZMin:DimZMax) !温位 real(8),intent(inout):: xza_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) !凝縮成分の混合比 real(8) :: xz_PotTemp_pre(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_PotTemp_nxt(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_MixRtNH3_pre(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_MixRtNH3_nxt(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_MixRtH2S_pre(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_MixRtH2S_nxt(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_MixRtNH4SH_pre(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_MixRtNH4SH_nxt(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_ExnerAll(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_PressAll(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_TempAll(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_Gamma(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_DelMixRt(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_Cond(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_Evap(DimXMin:DimXMax, DimZMin:DimZMax) integer :: i integer, parameter :: ItrNum = 2 !--------------------------------------------------------------------- ! 初期化 !--------------------------------------------------------------------- xz_PotTemp_nxt = 0.0d0 xz_MixRtNH3_nxt = 0.0d0 xz_MixRtH2S_nxt = 0.0d0 xz_MixRtNH4SH_nxt = 0.0d0 xz_Gamma = 0.0d0 xz_DelMixRt = 0.0d0 !湿潤飽和法では圧力は変化させない. xz_ExnerAll = xz_Exner + xz_ExnerBasicZ xz_PressAll = PressBasis * (xz_ExnerAll ** (CpDry / GasRDry)) !今までに得られた擾乱成分の値を暫定量とみなす. 添え字を追加 xz_MixRtNH3_pre(:,:) = xza_MixRt(:,:,NH3Num) + xza_MixRtBasicZ(:,:,NH3Num) xz_MixRtH2S_pre(:,:) = xza_MixRt(:,:,H2SNum) + xza_MixRtBasicZ(:,:,H2SNum) xz_MixRtNH4SH_pre(:,:) = xza_MixRt(:,:,NH4SHNum) + xza_MixRtBasicZ(:,:,NH4SHNum) xz_PotTemp_pre = xz_PotTemp AdjustNH4SH: do i = 1, ItrNum !--------------------------------------------------------------- ! 変数の初期化 !--------------------------------------------------------------- !温度 xz_TempAll = ( xz_PotTemp_pre + xz_PotTempBasicZ ) * xz_ExnerAll !規格化された反応熱 (NH4SH １kg に対する熱量) xz_Gamma = ReactHeatNH4SH * xz_EffMolWtBasicZ / ( xz_ExnerAll * CpDry ) !NH4SH の生成量 xz_DelMixRt = & & xz_DelMixRtNH4SH( & & xz_TempAll, xz_PressAll, xz_MixRtNH3_pre, xz_MixRtH2S_pre, & & MolWtWet(NH3Num), MolWtWet(H2SNum) & & ) xz_Cond = max( 0.0d0, xz_DelMixRt ) xz_Evap = max( 0.0d0, min( - xz_DelMixRt, xz_MixRtNH4SH_pre ) ) !--------------------------------------------------------------- ! より真に近い値を求める飽和蒸気圧から飽和混合比を求める !--------------------------------------------------------------- ! NH4SH の混合比を修正 xz_MixRtNH4SH_nxt = xz_MixRtNH4SH_pre + xz_Cond - xz_Evap ! DelPress を元に, NH3 と H2S の混合比を修正 xz_MixRtNH3_nxt = xz_MixRtNH3_pre - ( xz_Cond - xz_Evap ) & & * MolWtWet(NH3Num) / MolWtWet(NH4SHNum) xz_MixRtH2S_nxt = xz_MixRtH2S_pre - ( xz_Cond - xz_Evap ) & & * MolWtWet(H2SNum) / MolWtWet(NH4SHNum) !温位を修正 xz_PotTemp_nxt = xz_PotTemp_pre + xz_Gamma * ( xz_Cond - xz_Evap ) !ループを回すための変数変化 xz_PotTemp_pre = xz_PotTemp_nxt xz_MixRtNH3_pre = xz_MixRtNH3_nxt xz_MixRtH2S_pre = xz_MixRtH2S_nxt xz_MixRtNH4SH_pre = xz_MixRtNH4SH_nxt end do AdjustNH4SH xz_PotTemp = xz_PotTemp_nxt xza_MixRt(:,:,NH3Num) = xz_MixRtNH3_nxt - xza_MixRtBasicZ(:,:,NH3Num) xza_MixRt(:,:,H2SNum) = xz_MixRtH2S_nxt - xza_MixRtBasicZ(:,:,H2SNum) xza_MixRt(:,:,NH4SHNum) = xz_MixRtNH4SH_nxt - xza_MixRtBasicZ(:,:,NH4SHNum) end subroutine MoistAdjustNH4SH end Module MoistAdjust