!= Module MoistBuoyancy ! ! Authors:: SUGIYAMA Ko-ichiro ! Version:: $Id: moistbuoyancy.f90,v 1.15 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 MoistBuoyancy ! ! ! !モジュール読み込み use dc_message, only: MessageNotify use gridset, only: SpcNum, &! 化学種の数 & DimXMin, &! x 方向の配列の下限 & DimXMax, &! x 方向の配列の上限 & DimZMin, &! z 方向の配列の下限 & DimZMax, &! z 方向の配列の上限 & DelX, &! x 方向の格子点間隔 & DelZ ! z 方向の格子点間隔 use basicset, only: MolWtDry, &!乾燥成分の分子量 & Grav, &!重力加速度 & SpcWetID, &!凝縮成分のID & SpcWetSymbol, &!凝縮成分の名前 & MolWtWet, &!凝縮成分の分子量 & CpDry, &!乾燥成分の平均定圧比熱 [J/K kg] & xz_PotTempBasicZ, &!基本場の温位 & xz_EffMolWtBasicZ, &!基本場の分子量効果 & xz_ExnerBasicZ, &!基本場の無次元圧力 & xza_MixRtBasicZ !基本場の混合比 use ChemCalc, only: xz_LatentHeat !潜熱 use ChemData, only: ChemData_OneSpcID use average, only: xr_avr_xz, xz_avr_xr use differentiate_center2, only: xr_dz_xz !暗黙の型宣言禁止 implicit none !属性の指定 private !関数を public に設定 public MoistBuoy_Init public xz_BuoyMoistKm public xr_BuoyMolWt public xr_BuoyDrag !変数の定義 integer :: LoopNum = 0 integer :: GasCount = 0 integer :: GasNum(10) = 0 integer :: CloudNum(10) = 0 real(8) :: Cm = 2.0d-1 real(8) :: MixLen = 0.0d0 real(8), allocatable :: xz_MixRtBasicZPerMolWt(:,:) !基本場の混合比 / 分子量 の総和 real(8), allocatable :: xr_MixRtBasicZPerMolWt(:,:) !基本場の混合比 / 分子量 の総和 real(8), allocatable :: xz_MixRtBasicZ(:,:) !基本場の混合比の総和 real(8), allocatable :: xr_MixRtBasicZ(:,:) !基本場の混合比の総和 save GasNum, CloudNum, LoopNum save Cm, MixLen save xz_MixRtBasicZPerMolWt save xr_MixRtBasicZPerMolWt save xz_MixRtBasicZ save xr_MixRtBasicZ contains !!!------------------------------------------------------------------!!! subroutine MoistBuoy_Init( ) use ChemData, only: ChemData_OneSpcID !暗黙の型宣言禁止 implicit none !変数定義 integer :: s integer :: n1, n2 real(8), allocatable :: xza_MixRtBasicZPerMolWt(:,:,:) !基本場の混合比 / 分子量 !----------------------------------------------------------- ! 混合距離 !----------------------------------------------------------- MixLen = sqrt(DelX * DelZ) !----------------------------------------------------------- ! 雲粒と気体の ID の組を作る !----------------------------------------------------------- !化学種の中から雲粒を作るものを選び, その配列添え字と分子量を保管. LoopNum = 0 GasCount = 0 SelectCloud: do s = 1, SpcNum !'Cloud' という文字列が含まれるものの個数を数える n1 = index(SpcWetSymbol(s), '-Cloud' ) if (n1 /= 0) then LoopNum = LoopNum + 1 CloudNum(LoopNum)= s end if n2 = index(SpcWetSymbol(s), '-g' ) if (n2 /= 0) then GasCount = GasCount + 1 GasNum(GasCount)= s end if end do SelectCloud !----------------------------------------------------------- ! 確認 !----------------------------------------------------------- call MessageNotify( "M", & & "MoistBuoy_Init", "LoopNum = %d", i=(/LoopNum/) ) call MessageNotify( "M", & & "MoistBuoy_Init", "GasCount = %d", i=(/GasCount/)) call MessageNotify( "M", & & "MoistBuoy_Init", "GasNum = %d", i=(/GasNum/) ) call MessageNotify( "M", & & "MoistBuoy_Init", "CloudNum= % ad", i=(/CloudNum/)) if ( LoopNum /= GasCount ) then call MessageNotify( "E", & & "MoistBuoy_Init", "GasCount /= LoopNum" ) stop end if !----------------------------------------------------------- ! 作業配列の初期化. !----------------------------------------------------------- allocate( & & xza_MixRtBasicZPerMolWt(DimXMin:DimXMax, DimZMin:DimZMax, LoopNum), & & xr_MixRtBasicZPerMolWt(DimXMin:DimXMax, DimZMin:DimZMax), & & xz_MixRtBasicZPerMolWt(DimXMin:DimXMax, DimZMin:DimZMax), & & xr_MixRtBasicZ(DimXMin:DimXMax, DimZMin:DimZMax), & & xz_MixRtBasicZ(DimXMin:DimXMax, DimZMin:DimZMax) & & ) xza_MixRtBasicZPerMolWt = 0.0d0 xr_MixRtBasicZPerMolWt = 0.0d0 xz_MixRtBasicZPerMolWt = 0.0d0 xr_MixRtBasicZ = 0.0d0 xz_MixRtBasicZ = 0.0d0 call MessageNotify( "M", & & "MoistBuoy_Init", "MolWtWet = %f", d=(/pack(MolWtWet(:), .true.) /) ) do s = 1, LoopNum xza_MixRtBasicZPerMolWt(:,:,GasNum(s)) = & & xza_MixRtBasicZ(:,:,GasNum(s)) / MolWtWet(GasNum(s)) end do xz_MixRtBasicZPerMolWt = sum(xza_MixRtBasicZPerMolWt, 3) xr_MixRtBasicZPerMolWt = xr_avr_xz( sum(xza_MixRtBasicZPerMolWt, 3) ) xz_MixRtBasicZ = sum(xza_MixRtBasicZ, 3) xr_MixRtBasicZ = xr_avr_xz( sum(xza_MixRtBasicZ, 3) ) end subroutine MoistBuoy_Init !!!------------------------------------------------------------------------!!! function xz_BuoyMoistKm(xz_PotTemp, xz_Exner, xza_MixRt) ! !本来は NH4SH の凝結熱を考慮しないとならないのだが, ここでは !とりあえず無視する. ! !暗黙の型宣言禁止 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) :: xza_MixRtAll(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) !凝縮成分の混合比 real(8) :: xz_BuoyMoistKm(DimXMin:DimXMax, DimZMin:DimZMax) ! real(8) :: xza_LatentHeat(DimXMin:DimXMax, DimZMin:DimZMax, LoopNum) !潜熱 real(8) :: xz_TempAll(DimXMin:DimXMax, DimZMin:DimZMax) !温度 real(8) :: xz_EffHeat(DimXMin:DimXMax, DimZMin:DimZMax) ! real(8) :: xz_EffPotTemp(DimXMin:DimXMax, DimZMin:DimZMax) real(8) :: xz_EffMolWt(DimXMin:DimXMax, DimZMin:DimZMax) ! real(8) :: xza_MixRtPerMolWt(DimXMin:DimXMax, DimZMin:DimZMax, LoopNum) !混合比/分子量 integer :: s !温度, 圧力, 混合比の全量を求める !擾乱成分と平均成分の足し算 xz_TempAll = ( xz_PotTemp + xz_PotTempBasicZ ) * ( xz_Exner + xz_ExnerBasicZ ) xza_MixRtAll = xza_MixRtBasicZ + xza_MixRt xza_LatentHeat = 0.0d0 !作業配列の初期化. 気体のみ利用 do s = 1, LoopNum xza_MixRtPerMolWt(:,:,GasNum(s)) = & & xza_MixRt(:,:,GasNum(s)) / MolWtWet(GasNum(s)) end do !温度の効果 xz_EffPotTemp = xz_PotTemp / xz_PotTempBasicZ !分子量効果 + 引きづりの効果 xz_EffMolWt = & & + sum(xza_MixRtPerMolWt, 3) & & / ( 1.0d0 / MolWtDry + xz_MixRtBasicZPerMolWt ) & & - sum(xza_MixRt, 3) / ( 1.0d0 + xz_MixRtBasicZ ) !蒸気が蒸発する場合の潜熱を計算 ! 分子量の部分はいつでも効くが潜熱は飽和していないと効かないので, ! 雲の混合比がゼロの時には, 潜熱の寄与はゼロとなるように調節している do s = 1, LoopNum xza_LatentHeat(:,:,s) = & & xz_LatentHeat( SpcWetID(CloudNum(s)), xz_TempAll ) & & * xza_MixRtAll(:,:,GasNum(s)) & & * ( 5.0d-1 + sign( 5.0d-1, (xza_MixRtAll(:,:,CloudNum(s)) - 1.0d-4) ) ) end do xz_EffHeat = sum( xza_LatentHeat, 3 ) * xz_EffMolWtBasicZ & & / ( CpDry * xz_ExnerBasicZ ) !乱流拡散係数の時間発展式の浮力項を決める xz_BuoyMoistKm = & & - 3.0d0 * Grav * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) & & * xz_avr_xr( & & xr_dz_xz( & & xz_EffHeat & & + xz_PotTempBasicZ / xz_EffMolWtBasicZ & & * ( 1.0d0 + xz_EffPotTemp + xz_EffMolWt ) & & ) & & ) & & / ( 2.0d0 * xz_PotTempBasicZ / xz_EffMolWtBasicZ) end function xz_BuoyMoistKm !!!------------------------------------------------------------------------!!! function xr_BuoyMolWt(xza_MixRt) ! ! 鉛直方向の運動方程式に現れる浮力項のうち, ! 分子量の効果だけを求める ! !暗黙の型宣言禁止 implicit none !変数定義 real(8), intent(in) :: xza_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) !凝縮成分の混合比 real(8) :: xr_BuoyMolWt(DimXMin:DimXMax, DimZMin:DimZMax) !浮力項(分子量効果) real(8) :: xza_MixRtPerMolWt(DimXMin:DimXMax, DimZMin:DimZMax, LoopNum) !混合比/分子量 integer :: s !初期化 xr_BuoyMolWt = 0.0d0 !作業配列の初期化. 気体のみ利用 do s = 1, LoopNum xza_MixRtPerMolWt(:,:,GasNum(s)) = & & xza_MixRt(:,:,GasNum(s)) / MolWtWet(GasNum(s)) end do !浮力項の計算 xr_BuoyMolWt = & & + Grav * xr_avr_xz( sum(xza_MixRtPerMolWt, 3) ) & & / ( 1.0d0 / MolWtDry + xr_MixRtBasicZPerMolWt ) & & - Grav * xr_avr_xz( sum(xza_MixRt(:,:,1:LoopNum), 3) ) & & / ( 1.0d0 + xr_MixRtBasicZ ) end function xr_BuoyMolWt !!!------------------------------------------------------------------------!!! function xr_BuoyDrag(xza_MixRt) ! ! 鉛直方向の運動方程式に現れる浮力項のうち, ! 引きづりのの効果だけを求める ! !暗黙の型宣言禁止 implicit none !変数定義 real(8), intent(in) :: xza_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) !凝縮成分の混合比 real(8) :: xr_BuoyDrag(DimXMin:DimXMax, DimZMin:DimZMax) !浮力項(分子量効果) !初期化 xr_BuoyDrag = 0.0d0 !浮力項の計算 xr_BuoyDrag = & & - Grav * xr_avr_xz( sum(xza_MixRt(:,:,LoopNum+1:SpcNum), 3) ) & & / ( 1.0d0 + xr_MixRtBasicZ ) end function xr_BuoyDrag end module MoistBuoyancy