!= Module MoistBuoyancy_3d ! ! Authors:: SUGIYAMA Ko-ichiro, ODAKA Masatsugu ! Version:: $Id: moistbuoyancy_3d.f90,v 1.2 2011-03-28 01:48:54 sugiyama Exp $ ! Tag Name:: $Name: arare4-20120911 $ ! Copyright:: Copyright (C) GFD Dennou Club, 2008. All rights reserved. ! License:: See COPYRIGHT[link:../../COPYRIGHT] ! !== Overview ! ! !== Error Handling ! !== Bugs ! !== Note ! ! !== Future Plans ! ! module MoistBuoyancy_3d ! ! ! !モジュール読み込み use dc_types, only : DP use dc_message, only: MessageNotify use gridset_3d, only: SpcNum, &! 化学種の数 & DimXMin, &! x 方向の配列の下限 & DimXMax, &! x 方向の配列の上限 & DimYMin, &! y 方向の配列の下限 & DimYMax, &! y 方向の配列の上限 & DimZMin, &! z 方向の配列の下限 & DimZMax, &! z 方向の配列の上限 & x_dx, &! x 方向の格子点間隔 & y_dy, &! y 方向の格子点間隔 & z_dz ! z 方向の格子点間隔 use basicset_3d, only: MolWtDry, &!乾燥成分の分子量 & Grav, &!重力加速度 & SpcWetID, &!凝縮成分のID & MolWtWet, &!凝縮成分の分子量 & CpDry, &!乾燥成分の平均定圧比熱 [J/K kg] & xyz_PotTempBasicZ, &!基本場の温位 & xyz_EffMolWtBasicZ, &!基本場の分子量効果 & xyz_ExnerBasicZ, &!基本場の無次元圧力 & xyza_MixRtBasicZ !基本場の混合比 use moistset_3d, only: CondNum, &!凝結過程の数 & IdxCG, &!凝結過程(蒸気)の配列添え字 & IdxCC, &!凝結過程(雲)の配列添え字 & GasNum, &!気体の数 & IdxG, &!気体の配列添え字 & IdxNH3, &!NH3(蒸気)の配列添え字 & IdxH2S !H2S(蒸気)の配列添え字 use ChemCalc_3d, only: xyz_LatentHeat, &!潜熱 & ReactHeatNH4SH !NH4SH の反応熱 use xyz_base_module, only : xyz_avr_xyr, xyr_avr_xyz use xyz_deriv_module,only : xyr_dz_xyz use StoreBuoy_3d, only: StoreBuoyMolWt, StoreBuoyDrag !暗黙の型宣言禁止 implicit none !属性の指定 private !関数を public に設定 public MoistBuoy_Init public xyz_BuoyMoistKm public xyr_BuoyMolWt public xyr_BuoyDrag !変数の定義 real(DP) :: Cm = 2.0d-1 real(DP) :: MixLen = 0.0d0 real(DP), allocatable :: xyz_MixRtBasicZPerMolWt(:,:,:) !基本場の混合比 / 分子量 の総和 real(DP), allocatable :: xyr_MixRtBasicZPerMolWt(:,:,:) !基本場の混合比 / 分子量 の総和 real(DP), allocatable :: xyz_MixRtBasicZ(:,:,:) !基本場の混合比の総和 real(DP), allocatable :: xyr_MixRtBasicZ(:,:,:) !基本場の混合比の総和 save Cm, MixLen save xyz_MixRtBasicZPerMolWt save xyr_MixRtBasicZPerMolWt save xyz_MixRtBasicZ save xyr_MixRtBasicZ contains !!!------------------------------------------------------------------!!! subroutine MoistBuoy_Init( ) !暗黙の型宣言禁止 implicit none !変数定義 integer :: s real(DP), allocatable :: xyza_MixRtBasicZPerMolWt(:,:,:,:) !基本場の混合比 / 分子量 real(DP) :: MinDelX real(DP) :: MinDelY real(DP) :: MinDelZ !----------------------------------------------------------- ! 混合距離 !----------------------------------------------------------- !混合距離 MinDelX = minval(x_dx) MinDelY = minval(y_dy) MinDelZ = minval(z_dz) ! MixLen = ( MinDelX * MinDelY * MinDelZ ) ** (1.0d0/3.0d0) MixLen = min( MinDelZ, min( MinDelX, MinDelY ) ) !----------------------------------------------------------- ! 作業配列の初期化. !----------------------------------------------------------- allocate( & & xyza_MixRtBasicZPerMolWt(DimXMin:DimXMax,DimYMin:DimYMax, DimZMin:DimZMax, GasNum), & & xyr_MixRtBasicZPerMolWt(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax), & & xyz_MixRtBasicZPerMolWt(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax), & & xyr_MixRtBasicZ(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax), & & xyz_MixRtBasicZ(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) & & ) xyza_MixRtBasicZPerMolWt = 0.0d0 xyr_MixRtBasicZPerMolWt = 0.0d0 xyz_MixRtBasicZPerMolWt = 0.0d0 xyr_MixRtBasicZ = 0.0d0 xyz_MixRtBasicZ = 0.0d0 call MessageNotify( "M", & & "MoistBuoy_Init", "MolWtWet = %f", d=(/pack(MolWtWet(:), .true.) /) ) do s = 1, GasNum xyza_MixRtBasicZPerMolWt(:,:,:,s) = & & xyza_MixRtBasicZ(:,:,:,IdxG(s)) / MolWtWet(IdxG(s)) end do xyz_MixRtBasicZPerMolWt = sum(xyza_MixRtBasicZPerMolWt, 4) xyr_MixRtBasicZPerMolWt = xyr_avr_xyz( sum(xyza_MixRtBasicZPerMolWt, 4) ) xyz_MixRtBasicZ = sum(xyza_MixRtBasicZ, 4) xyr_MixRtBasicZ = xyr_avr_xyz( sum(xyza_MixRtBasicZ, 4) ) end subroutine MoistBuoy_Init !!!------------------------------------------------------------------------!!! function xyz_BuoyMoistKm(xyz_PotTemp, xyz_Exner, xyza_MixRt) ! !暗黙の型宣言禁止 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) :: xyza_MixRtAll(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum) !凝縮成分の混合比 real(DP) :: xyza_MixRtAll2(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum) !凝縮成分の混合比 real(DP) :: xyz_BuoyMoistKm(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) ! real(DP) :: xyza_LatentHeat(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, GasNum) !潜熱 real(DP) :: xyz_TempAll(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !温度 real(DP) :: xyz_EffHeat(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) ! real(DP) :: xyz_EffPotTemp(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !浮力に対する温度差の寄与 real(DP) :: xyz_EffMolWt(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !浮力に対する分子量差の寄与 real(DP) :: xyza_MixRtPerMolWt(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, GasNum) !混合比/分子量 integer :: s !温度, 圧力, 混合比の全量を求める !擾乱成分と平均成分の足し算 xyz_TempAll = ( xyz_PotTemp + xyz_PotTempBasicZ ) * ( xyz_Exner + xyz_ExnerBasicZ ) xyza_MixRtAll = xyza_MixRtBasicZ + xyza_MixRt xyza_LatentHeat = 0.0d0 !作業配列の初期化. 気体のみ利用 do s = 1, GasNum xyza_MixRtPerMolWt(:,:,:,s) = xyza_MixRt(:,:,:,IdxG(s)) / MolWtWet(IdxG(s)) end do !温度の効果 xyz_EffPotTemp = xyz_PotTemp / xyz_PotTempBasicZ !分子量効果 + 引きづりの効果 xyz_EffMolWt = & & + sum(xyza_MixRtPerMolWt, 4) & & / ( 1.0d0 / MolWtDry + xyz_MixRtBasicZPerMolWt ) & & - sum(xyza_MixRt, 4) / ( 1.0d0 + xyz_MixRtBasicZ ) !蒸気が蒸発する場合の潜熱を計算 ! 分子量の部分はいつでも効くが潜熱は飽和していないと効かないので, ! 雲の混合比がゼロの時には, 潜熱の寄与はゼロとなるように調節している xyza_MixRtAll2 = xyza_MixRtAll ! xzya_MixRtAll2(:,:,:,IdxNH3) = & ! & xyza_MixRtAll(:,:,:,IdxNH3) - xyza_MixRtAll(:,:,:,IdxH2S) do s = 1, CondNum xyza_LatentHeat(:,:,:,s) = & & xyz_LatentHeat( SpcWetID(IdxCC(s)), xyz_TempAll ) & & * xyza_MixRtAll2(:,:,:,IdxCG(s)) & & * ( 5.0d-1 + sign( 5.0d-1, (xyza_MixRtAll2(:,:,:,IdxCC(s)) - 1.0d-4) ) ) end do xyz_EffHeat = ( sum( xyza_LatentHeat, 4 ) * xyz_EffMolWtBasicZ & ! & + ReactHeatNH4SH * & ! & (xyza_MixRtAll(:,:,:,IdxH2S) + xyza_MixRtAll(:,:,:,IdxH2S)) & & ) / ( CpDry * xyz_ExnerBasicZ ) !乱流拡散係数の時間発展式の浮力項を決める xyz_BuoyMoistKm = & & - 3.0d0 * Grav * ( Cm ** 2.0d0 ) * ( MixLen ** 2.0d0 ) & & * xyz_avr_xyr( & & xyr_dz_xyz( & & xyz_EffHeat & & + xyz_PotTempBasicZ / xyz_EffMolWtBasicZ & & * ( 1.0d0 + xyz_EffPotTemp + xyz_EffMolWt ) & & ) & & ) & & / ( 2.0d0 * xyz_PotTempBasicZ / xyz_EffMolWtBasicZ) end function xyz_BuoyMoistKm !!!------------------------------------------------------------------------!!! function xyr_BuoyMolWt(xyza_MixRt) ! ! 鉛直方向の運動方程式に現れる浮力項のうち, ! 分子量の効果だけを求める ! !暗黙の型宣言禁止 implicit none !変数定義 real(DP), intent(in) :: xyza_MixRt(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum) !凝縮成分の混合比 real(DP) :: xyr_BuoyMolWt(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !浮力項(分子量効果) real(DP) :: xyza_MixRtPerMolWt(DimXMin:DimXMax,DimYMin:DimYMax, DimZMin:DimZMax, GasNum) !混合比/分子量 integer :: s !初期化 xyr_BuoyMolWt = 0.0d0 !作業配列の初期化. 気体のみ利用 do s = 1, GasNum xyza_MixRtPerMolWt(:,:,:,s) = xyza_MixRt(:,:,:,IdxG(s)) / MolWtWet(IdxG(s)) end do !浮力項の計算 xyr_BuoyMolWt = & & + Grav * xyr_avr_xyz( sum(xyza_MixRtPerMolWt, 4) ) & & / ( 1.0d0 / MolWtDry + xyr_MixRtBasicZPerMolWt ) & & - Grav * xyr_avr_xyz( sum(xyza_MixRt(:,:,:,1:GasNum), 4) )& & / ( 1.0d0 + xyr_MixRtBasicZ ) call StoreBuoyMolWt(xyz_avr_xyr(xyr_BuoyMolWt)) end function xyr_BuoyMolWt !!!------------------------------------------------------------------------!!! function xyr_BuoyDrag(xyza_MixRt) ! ! 鉛直方向の運動方程式に現れる浮力項のうち, ! 引きづりのの効果だけを求める ! !暗黙の型宣言禁止 implicit none !変数定義 real(DP), intent(in) :: xyza_MixRt(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax, SpcNum) !凝縮成分の混合比 real(DP) :: xyr_BuoyDrag(DimXMin:DimXMax,DimYMin:DimYMax,DimZMin:DimZMax) !浮力項(分子量効果) !初期化 xyr_BuoyDrag = 0.0d0 !浮力項の計算 xyr_BuoyDrag = & & - Grav * xyr_avr_xyz( sum(xyza_MixRt(:,:,:,GasNum+1:SpcNum), 4) ) & & / ( 1.0d0 + xyr_MixRtBasicZ ) call StoreBuoyDrag(xyz_avr_xyr(xyr_BuoyDrag)) end function xyr_BuoyDrag end module MoistBuoyancy_3d