非静力学モデル deepconv/arare.
program arare
!
!非静力学モデル deepconv/arare.
!
!----- モジュール読み込み ------
!----- 型宣言, 文字列処理 ----
use dc_types, only : STRING
use dc_string, only : StoA
!----- メッセージ出力 -----
use dc_message, only: MessageNotify
! コマンドライン引数解釈
use argset, only : argset_init
!----- 管理モジュール -----
! 化学量計算モジュール
use ChemCalc, only: ChemCalc_init
use chemdata, only: chemdata_init
! 入出力ファイル名管理モジュール
use fileset, only : fileset_init, InitFile
! デバッグ出力管理モジュール
use debugset, only : debugset_init
! 時間管理モジュール
use timeset, only : timeset_init, NstepLong, NstepShort, DelTimeLong, DelTimeShort, NstepDisp
! ! 格子点管理モジュール
! use gridset, only : gridset_init, &
! & DimXMin, DimXMax, DimZMin, DimZMax, SpcNum
! 格子点管理モジュール
use gridset, only : gridset_init, DimXMin, DimXMax, DimZMin, DimZMax, SpcNum, RegXMin, RegXMax, RegZMin, RegZMax, DelX, DelZ, Xmin, Xmax, s_Z
! ! 基本場設定モジュール
! use basicset, only : basicset_init, &
! & xz_DensBasicZ, xza_MixRtBasicZ, xz_PotTempBasicZ, &
! & xz_VelSoundBasicZ, xz_ExnerBasicZ
! 基本場設定モジュール
use basicset, only : basicset_init, xz_DensBasicZ, xza_MixRtBasicZ, xz_PotTempBasicZ, xz_VelSoundBasicZ, xz_ExnerBasicZ, CpDry, GasRDry, Grav, PressSfc
! 平均操作モジュール(スカラー格子点, フラックス格子点間の変換)
use average, only : xz_avr_pz, xz_avr_xr
! 積算値管理モジュール
use StorePotTemp, only : StorePotTemp_init, StorePotTempClean, StorePotTempCond
use StoreDensCloud, only : StoreDensCloud_init, StoreDensCloudClean, StoreDensCloudCond, StoreDensCloudFill, StoreDensCloudFillZero, StoreDensCloudMean, xz_DensCloudAdv, xz_DensCloudTurb, xz_DensCloudDiff, xz_DensCloudCond, xz_DensCloudFill, xz_DensCloudFillZero, xz_DensCloudAdv2, xz_DensCloudTurb2, xz_DensCloudDiff2, xz_DensCloudCond2, xz_DensCloudFill2, xz_DensCloudFillZero2, DensCloudAdvTendSum, DensCloudTurbTendSum, DensCloudDiffTendSum, DensCloudCondTendSum, DensCloudFillTendSum, DensCloudFillZeroTendSum
use StoreMixRt, only : StoreMixRt_init, StoreMixRtClean, StoreMixRtCond, StoreMixRtFill1, StoreMixRtFill2
use StoreBuoy, only : StoreBuoy_init, StoreBuoyClean
use StoreStab, only : StoreStab_init, StoreStabClean
! 湿潤ルーチン設定モジュール
use moistset, only: moistset_init
!----- 下請けモジュール -----
! 数値摩擦計算モジュール
use damping, only : damping_init, DampSponge_xz, DampSponge_pz, DampSponge_xr
! 時間積分フィルターモジュール
use timefilter, only : AsselinFilter_xz, AsselinFilter_xr, AsselinFilter_pz, AsselinFilter_xza
! 境界条件適用モジュール
use boundary, only : BoundaryXCyc_xz, BoundaryZSym_xz, BoundaryXCyc_xza, BoundaryZSym_xza, BoundaryXCyc_pz, BoundaryZSym_pz, BoundaryXCyc_xr, BoundaryZAntiSym_xr
! CFL 条件確認モジュール
use cflcheck, only : CFLCheckTimeShort, CFLCheckTimeLongVelX, CFLCheckTimeLongVelZ
! 負の湿潤量の補填計算モジュール
! use fillnegative, only : FillNegative_init, xza_FillNegative_xza
use fillnegative, only : xz_FillNegativeMMC_xz
! 主成分凝結対流用
!----- 入出力モジュール -----
! リスタートファイル入出力モジュール
use RestartFileIO, only : ReStartFile_Open, ReStartFile_OutPut, ReStartFile_Close, ReStartFile_Get
! ヒストリファイル入出力モジュール
use HistoryFileIO, only : HistoryFile_Open, HistoryFile_OutPut, HistoryFile_Close
!----- 力学過程 -----
! 力学過程計算用関数モジュール
use DynFunc, only : xz_AdvScalar, xz_AdvKm, xza_AdvScalar, pz_AdvVelX, xr_Buoy, xr_AdvVelZ, pz_GradPi
! 力学過程陰解法計算用関数モジュール
! use DynImpFunc, only : xz_Exner_init, xz_Exner, xr_GradPi
! 微量成分凝結用
use DynImpFuncMarscond2, only : xz_Exner_init, xz_Exner, xr_GradPi
! 主成分凝結用
!----- 物理過程 -----
! 数値拡散計算用モジュール
use NumDiffusion, only : NumDiffusion_Init, xz_NumDiffScalar, xz_NumDiffKm, xza_NumDiffScalar, pz_NumDiffVelX, xr_NumDiffVelZ
! 乱流拡散計算用モジュール
use Turbulence, only : Turbulence_Init, xz_TurbScalar, xza_TurbScalar, pz_TurbVelX, xr_TurbVelZ , xz_ShearKm , xz_DispKm, xz_DispHeat
! 放射強制計算用モジュール
use Radiation, only : Radiation_init, xz_RadHeatBalance !加熱率と冷却率が各時刻で
!釣り合うよう調整
! & xz_RadHeatConst !加熱率は一定
! 地表フラックス計算用モジュール
! use HeatFlux, only : xz_HeatFluxBulk, xz_MixRtFluxBulk
! use HeatFlux, only : xz_HeatFluxDiff, xza_MixRtFluxDiff
use HeatFlux_N1994, only : xz_HeatFluxBulk, xz_MixRtFluxBulk
! 湿潤飽和調節法計算用モジュール
use MoistAdjust, only : MoistAdjustSvapPress, MoistAdjustNH4SH
! 雲物理パラメタリゼーション
use WarmRainPrm, only : WarmRainPrm_Init, xz_Rain2GasHeat, xza_Rain2Gas, xza_Rain2GasNH4SH, xz_Rain2GasHeatNH4SH, xza_Cloud2Rain, xza_FallRain
! (2008/06/27 山下達也追加)
! 雲物理パラメータの初期化サブルーチン
use cloudset, only : cloudset_init
! 湿潤気塊の浮力計算用モジュール
use MoistBuoyancy,only : MoistBuoy_Init, xz_BuoyMoistKm, xr_BuoyMolWt, xr_BuoyDrag
! 安定度の計算
use ECCM, only: ECCM_Stab
!暗黙の型宣言禁止
implicit none
!内部変数
character(80) :: cfgfile
real(8), allocatable :: pz_VelXBl(:,:)
real(8), allocatable :: pz_VelXNl(:,:)
real(8), allocatable :: pz_VelXAl(:,:)
real(8), allocatable :: pz_VelXNs(:,:)
real(8), allocatable :: pz_VelXAs(:,:)
real(8), allocatable :: xr_VelZBl(:,:)
real(8), allocatable :: xr_VelZNl(:,:)
real(8), allocatable :: xr_VelZAl(:,:)
real(8), allocatable :: xr_VelZNs(:,:)
real(8), allocatable :: xr_VelZAs(:,:)
real(8), allocatable :: xz_ExnerBl(:,:)
real(8), allocatable :: xz_ExnerNl(:,:)
real(8), allocatable :: xz_ExnerAl(:,:)
real(8), allocatable :: xz_ExnerNs(:,:)
real(8), allocatable :: xz_ExnerAs(:,:)
real(8), allocatable :: xz_ExnerSum(:,:)
real(8), allocatable :: xz_PotTempWork(:,:)
real(8), allocatable :: xz_PotTempBl(:,:)
real(8), allocatable :: xz_PotTempNl(:,:)
real(8), allocatable :: xz_PotTempAl(:,:)
real(8), allocatable :: xz_PotTempNs(:,:)
real(8), allocatable :: xz_PotTempAs(:,:)
real(8), allocatable :: xz_PotTempSum(:,:)
real(8), allocatable :: xz_TempSum(:,:)
real(8), allocatable :: xz_KmBl(:,:)
real(8), allocatable :: xz_KmNl(:,:)
real(8), allocatable :: xz_KmAl(:,:)
real(8), allocatable :: xz_KhBl(:,:)
real(8), allocatable :: xz_KhNl(:,:)
real(8), allocatable :: xz_KhAl(:,:)
real(8), allocatable :: xza_MixRtWork(:,:,:)
real(8), allocatable :: xza_MixRtBl(:,:,:)
real(8), allocatable :: xza_MixRtNl(:,:,:)
real(8), allocatable :: xza_MixRtAl(:,:,:)
real(8), allocatable :: pz_AccelVelXNl(:,:)
real(8), allocatable :: xr_AccelVelZNl(:,:)
real(8), allocatable :: xza_DelMixRt(:,:,:)
! 長い時間ステップで評価すべき温位の項を格納する配列
! 2008/06/20 山下達也 追加
real(8), allocatable :: xz_TendPotTempNl(:,:)
real(8) :: Time
real(8) :: ReStartTime(2)
real(8), allocatable :: DelTimeLFrog(:)
real(8) :: DelTimeEular
integer, allocatable :: NStepEular(:)
integer :: NStepLFrog
integer :: t, tau, t1, t2, k
! 以下の変数は山下が追加(2008/05/07)
! 単位質量あたりの潜熱
real(8), allocatable :: xz_LatHeatPerMassNl(:,:)
! 過飽和度
real(8), allocatable :: xz_SatRatioBl(:,:)
real(8), allocatable :: xz_SatRatioNl(:,:)
real(8), allocatable :: xz_SatRatioAl(:,:)
real(8), allocatable :: xz_SatRatioNs(:,:)
real(8), allocatable :: xz_SatRatioAs(:,:)
! 雲密度
real(8), allocatable :: xz_DensCloudBl(:,:)
real(8), allocatable :: xz_DensCloudNl(:,:)
real(8), allocatable :: xz_DensCloudAl(:,:)
real(8), allocatable :: xz_DensCloudNs(:,:)
real(8), allocatable :: xz_DensCloudAs(:,:)
real(8), allocatable :: xz_DensCloudWorkAs(:,:) ! 作業変数
real(8) :: DensCloudMin ! 負の雲密度穴埋めの為の作業変数
! 凝結量
real(8), allocatable :: xz_MassCondNs(:,:)
real(8), allocatable :: xz_MassCondNl(:,:)
! 凝結熱による温位変化率
real(8), allocatable :: xz_QCond(:,:)
! 放射による温位変化率
real(8), allocatable :: xz_QRad(:,:)
! 熱散逸による温位変化率
real(8), allocatable :: xz_QDis(:,:)
! 作業変数
real(8), allocatable :: worknum(:,:)
real(8) :: workcloud
! 各格子点の質量密度(擾乱成分の寄与)
real(8), allocatable :: xz_MassDens(:,:)
! 領域全体の気相の質量
real(8) :: MassTotal
! 温位 2 次移流項の全質量の時間変化に対する寄与
real(8), allocatable :: xz_MassTend(:,:)
real(8) :: MassTend
! 領域全体の運動エネルギー
real(8) :: KineticEnergyTotal
! 各格子点の運動エネルギー
real(8), allocatable :: xz_KineticEnergy(:,:)
! 領域全体の雲の質量
real(8) :: CloudMassTotal
! 領域全体の雲の質量(fillnegative の前, write 文でチェック)
real(8) :: CloudMassTotal01
! 領域全体の雲の質量(fillnegative の後, write 文でチェック)
real(8) :: CloudMassTotal02
! 各格子点のポテンシャルエネルギー
real(8), allocatable :: xz_PotentialEnergy(:,:)
! 領域全体のポテンシャルエネルギー
real(8) :: PotentialEnergyTotal
! 各格子点の弾性エネルギー(1 次量)
real(8), allocatable :: xz_ElasticEnergyFO(:,:)
! 領域全体の弾性エネルギー(1 次量)
real(8) :: ElasticEnergyFOTotal
! 各格子点の弾性エネルギー(2 次量)
real(8), allocatable :: xz_ElasticEnergySO(:,:)
! 領域全体の弾性エネルギー(2 次量)
real(8) :: ElasticEnergySOTotal
! z 座標(2 次元配列)
real(8), allocatable :: xz_Z(:,:)
!コマンドライン引数の解釈
! NAMELIST ファイル名の読み込み
call argset_init(cfgfile)
!デバッグ設定
call debugset_init(cfgfile)
!物質特性の初期化
call chemdata_init()
!時刻に関する設定の初期化
! NAMELIST から必要な情報を読み取り, 時間関連の変数の設定を行う.
call timeset_init(cfgfile)
!格子点情報の初期化
! NAMELIST から情報を得て, 格子点を計算する
call gridset_init(cfgfile)
!化学計算ルーチンの初期化
call chemcalc_init()
!基本場の情報の初期化
! NAMELIST から情報を得て, 基本場を設定する.
call basicset_init(cfgfile)
!I/O ファイル名の初期化
! NAMELIST ファイル名を指定し, deepconv/arare の
! 出力ファイル名を NAMELIST から得る
call fileset_init(cfgfile)
!雲物理パラメータの初期化サブルーチン
call cloudset_init(cfgfile)
!湿潤ルーチンの共有変数の初期化
call moistset_init()
!write(*,*) "OK"
!積算値を保管するためのモジュールの初期化
! NAMELIST から情報を得て, 基本場を設定する.
call StorePotTemp_init( )
call StoreDensCloud_init( )
call StoreMixRt_init( )
call StoreBuoy_init( )
call StoreStab_init( )
!write(*,*) "OK"
!内部変数の初期化. とりあえずゼロを入れて値を確定させておく.
call ArareAlloc
!write(*,*) "OK"
!予報変数の初期化
! ReStartFile が設定されている場合にはファイルを読み込み,
! 設定されていない場合には, デフォルトの基本場と擾乱場を作る.
write(*,*) InitFile
if (trim(InitFile) /= '') then
!基本場, 擾乱場の初期値を netCDF ファイルから取得する.
call ReStartFile_Get( ReStartTime, xz_PotTempBl, xz_ExnerBl, pz_VelXBl, xr_VelZBl, xza_MixRtBl, xz_KmBl, xz_KhBl, xz_DensCloudBl, xz_SatRatioBl, xz_PotTempNl, xz_ExnerNl, pz_VelXNl, xr_VelZNl, xza_MixRtNl, xz_KmNl, xz_KhNl, xz_DensCloudNl, xz_SatRatioNl )
!write(*,*) "OK"
else
!デフォルト設定の基本場, 擾乱場を作成する.
call BasicEnv()
call DisturbEnv(cfgfile, xz_PotTempBl, xz_ExnerBl, pz_VelXBl, xr_VelZBl, xza_MixRtBl, xz_KmBl, xz_KhBl, xz_DensCloudBl, xz_SatRatioBl )
! 1 ループ目だけは, bl と nl の値を同じにしておく.
! 1 ステップ目はオイラー法で解く必要があるが,
! 1 ステップ目とそれ以外のステップを別々にコーディングしたくない為
xz_PotTempNl = xz_PotTempBl
xz_ExnerNl = xz_ExnerBl
pz_VelXNl = pz_VelXBl
xr_VelZNl = xr_VelZBl
xza_MixRtNl = xza_MixRtBl
xz_KmNl = xz_KmBl
xz_KhNl = xz_KhBl
xz_DensCloudNl = xz_DensCloudBl
xz_SatRatioNl = xz_SatRatioBl
xz_PotTempSum = xz_PotTempBl + xz_PotTempBasicZ
xz_ExnerSum = xz_ExnerBl + xz_ExnerBasicZ
xz_TempSum = xz_PotTempSum * xz_ExnerSum
xz_KineticEnergy = 0.5d0 * PressSfc * DelX * DelZ / GasRDry * (xz_ExnerBasicZ)**((CpDry - GasRDry)/GasRDry) / (xz_PotTempBasicZ) * ( xz_avr_pz(pz_VelXNl)**2.0d0 + xz_avr_xr(xr_VelZNl)**2.0d0 )
KineticEnergyTotal = sum( xz_KineticEnergy(RegXMin:RegXMax,RegZMin:RegZMax) )
xz_MassDens = (xz_ExnerBasicZ + xz_ExnerNl)**( (CpDry - GasRDry)/GasRDry ) / (xz_PotTempBasicZ + xz_PotTempNl) - (xz_ExnerBasicZ )**( (CpDry - GasRDry)/GasRDry ) / xz_PotTempBasicZ
MassTotal = (Xmax - Xmin) * PressSfc /Grav * ( xz_ExnerBasicZ(RegXMin,RegZMin)**(CpDry/GasRDry) - xz_ExnerBasicZ(RegXMin,RegZMax)**(CpDry/GasRDry) ) + PressSfc * DelX * DelZ / GasRDry * sum( xz_MassDens(RegXMin:RegXMax,RegZMin:RegZMax))
CloudMassTotal = DelX * DelZ * sum( xz_DensCloudNl(RegXMin:RegXMax,RegZMin:RegZMax) )
do k = RegZMin, RegZMax
xz_Z(:,k) = s_Z(k)
end do
xz_PotentialEnergy = - xz_DensBasicZ * xz_PotTempNl * xz_Z / xz_PotTempBasicZ
PotentialEnergyTotal = DelX * DelZ * sum( xz_PotentialEnergy(RegXMin:RegXMax,RegZMin:RegZMax) )
xz_MassTend = - xz_DensBasicZ * xz_TendPotTempNl / xz_PotTempBasicZ
MassTend = DelX * DelZ * sum( xz_MassTend(RegXMin:RegXMax,RegZMin:RegZMax) )
!write(*,*) "OK"
end if
!write(*,*) "OK"
!----------------------------------------------------------------------
! パッケージ型モジュールの初期化
! デフォルトの値から変更する必要のあるルーチンのみ初期化
!----------------------------------------------------------------------
call Damping_Init( cfgfile ) !波の減衰係数の初期化
call NumDiffusion_Init() !数値拡散項の初期化
call Turbulence_Init() !乱流計算の初期化
call WarmRainPrm_Init( cfgfile ) !暖かい雨のパラメタリゼーションの初期化
! call FillNegative_Init( xza_MixRtBasicZ, xz_DensBasicZ)
! !移流による負の量の処理
call Radiation_Init( cfgfile ) !放射強制の初期化
call MoistBuoy_Init() !分子量に対する浮力計算ルーチンの初期化
call xz_Exner_Init() !陰解法の初期化
!----------------------------------------------------------------------
! 時刻とループ回数の初期化
!----------------------------------------------------------------------
NstepLFrog = NstepLong
NstepEular = NstepShort
DelTimeLFrog = DelTimeLong * 2.0d0
DelTimeEular = DelTimeShort
if ( trim(InitFile) /= '') then
Time = ReStartTime(2) !リスタート開始時刻
else
Time = 0.0d0 !計算開始時刻
NstepEular(1) = NstepShort /2 ! 1 ループ目だけ
DelTimeLFrog(1) = DelTimeLong ! 1 ループ目だけ
end if
!----------------------------------------------------------------
! ファイル入出力
!----------------------------------------------------------------
!ファイルオープン
call HistoryFile_Open( )
if ( trim(InitFile) /= '') then
call HistoryFile_Output( ReStartTime(2), xz_PotTempNl, xz_PotTempSum, xz_TempSum, xz_ExnerNl, pz_VelXNl, xr_VelZNl, xza_MixRtNl, xz_KmNl, xz_KhNl, xz_DensCloudNl, xz_SatRatioNl, MassTotal, MassTend, KineticEnergyTotal, CloudMassTotal, PotentialEnergyTotal, ElasticEnergyFOTotal, ElasticEnergySOTotal )
end if
!時刻 0 の場合には出力
if ( Time == 0 ) then
call HistoryFile_Output( Time, xz_PotTempNl, xz_PotTempSum, xz_TempSum, xz_ExnerNl, pz_VelXNl, xr_VelZNl, xza_MixRtNl, xz_KmNl, xz_KhNl, xz_DensCloudNl, xz_SatRatioNl, MassTotal, MassTend, KineticEnergyTotal, CloudMassTotal, PotentialEnergyTotal, ElasticEnergyFOTotal, ElasticEnergySOTotal )
end if
!----------------------------------------------------------------------
! 設定のチェック
!----------------------------------------------------------------------
!CFL 条件のチェック
call CFLCheckTimeShort( xz_VelSoundBasicZ )
! write(*,*) "OK"
!----------------------------------------------------------------------
! 数値積分
!----------------------------------------------------------------------
call MessageNotify( "M", "main", "Time integration" )
do t1 = 1, NstepLFrog / NstepDisp
do t2 = 1, NstepDisp
t = (t1 - 1) * NstepDisp + t2
!時刻の設定
Time = Time + DelTimeLong
! write(*,*) "OK"
!----------------------------------------------------------------
! 渦粘性係数, 渦拡散係数を求める.
!----------------------------------------------------------------
xz_KmAl = xz_KmBl + DelTimeLFrog(t) * ( + xz_AdvKm(xz_KmNl, pz_VelXNl, xr_VelZNl) + xz_BuoyMoistKm(xz_PotTempBl, xz_ExnerBl, xza_MixRtBl) + xz_ShearKm(xz_KmBl, pz_VelXBl, xr_VelZBl) + xz_NumDiffKm(xz_KmBl) + xz_DispKm(xz_KmBl) )
! write(*,*) "OK"
!値の上限下限の設定
! * 値は正になることを保証する
! * 値の上限は 800 とする. 中島健介(1994, 学位論文)参照
xz_KmAl = max( 0.0d0, min( xz_KmAl, 800.0d0 ) )
! xz_KmAl = max( 0.0d0, min( xz_KmAl, 300.0d0 ) )
!境界条件
call BoundaryXCyc_xz( xz_KmAl )
call BoundaryZSym_xz( xz_KmAl )
!渦拡散定数を決める
xz_KhAl = 3.0d0 * xz_KmAl
! write(*,*) "OK"
! !----------------------------------------------------------------
! ! 温位の移流計算.
! !----------------------------------------------------------------
! !時間積分
! xz_PotTempAl = &
! & xz_PotTempBl &
! & + DelTimeLFrog(t) &
! & * ( &
! & + xz_AdvScalar( xz_PotTempNl, pz_VelXNl, xr_VelZNl) &
! & + xz_AdvScalar( xz_PotTempBasicZ, pz_VelXNl, xr_VelZNl) &
! & + xz_TurbScalar(xz_PotTempBl, xz_KhBl) &
! & + xz_TurbScalar(xz_PotTempBasicZ, xz_KhBl) &
! & + xz_NumDiffScalar(xz_PotTempBl) &
! & + xz_DispHeat( xz_KmBl ) &
! & + xz_RadHeatConst( xz_ExnerBl ) &
! & + xz_HeatFluxDiff( xz_PotTempNl ) &
!!! & + xz_HeatFluxBulk( xz_PotTempNl ) &
!!! & + xz_NewtonCool( xz_PotTempBl ) &
! & )
!
! !境界条件
! call BoundaryXCyc_xz( xz_PotTempAl )
! call BoundaryZSym_xz( xz_PotTempAl )
! !----------------------------------------------------------------
! ! 凝縮成分の混合比の移流計算.
! !----------------------------------------------------------------
! xza_MixRtAl = &
! & xza_MixRtBl &
! & + DelTimeLFrog(t) &
! & * ( &
! & + xza_AdvScalar(xza_MixRtNl, pz_VelXNl, xr_VelZNl)&
! & + xza_AdvScalar(xza_MixRtBasicZ, pz_VelXNl, xr_VelZNl)&
! & + xza_TurbScalar(xza_MixRtBl, xz_KhBl) &
! & + xza_TurbScalar(xza_MixRtBasicZ,xz_KhBl) &
! & + xza_NumDiffScalar(xza_MixRtBl) &
! & + xza_FallRain(xza_MixRtBl) &
! & + xza_MixRtFluxDiff(xza_MixRtNl) &
!!! & + xza_MixRtFluxBulk(xza_MixRtNl) &
! & )
!
! !移流によって負になった部分を埋める
! xza_MixRtWork = xza_MixRtAl
! xza_MixRtAl = xza_FillNegative_xza( xza_MixRtWork )
!
! !埋めた/削った量を保管
! call StoreMixRtFill1( (xza_MixRtAl - xza_MixRtWork) / DelTimeLFrog(t) )
!
! !境界条件
! call BoundaryXCyc_xza( xza_MixRtAl )
! call BoundaryZSym_xza( xza_MixRtAl )
!
!
! !-------------------------------------------------------------
! ! 暖かい雨のパラメタリゼーション.
! ! 雲<-->雨 の変換を行う
! !-------------------------------------------------------------
! !雲から雨への変換分を追加する.
! ! xza_Cloud2Rain 関数は, 引数として時間刻みを取ることで,
! ! 時間積分値を出力する.
!
! !これまでの値を作業配列に保管
! xza_MixRtWork = xza_MixRtAl
!
! !雨への変化量を計算
! xza_MixRtAl = xza_MixRtWork &
! & + xza_Cloud2Rain( xza_MixRtAl, DelTimeLFrog(t) )
!
! !雲から雨への変換量を保管
! call StoreMixRtCond( (xza_MixRtAl - xza_MixRtWork) / DelTimeLFrog(t) )
!
! !境界条件
! call BoundaryXCyc_xza( xza_MixRtAl )
! call BoundaryZSym_xza( xza_MixRtAl )
!
!
! !-------------------------------------------------------------
! ! 湿潤飽和調節法
! ! 水蒸気<-->雲の変換を行う.
! !-------------------------------------------------------------
! !これまでの値を作業配列に保管
! xz_PotTempWork = xz_PotTempAl
! xza_MixRtWork = xza_MixRtAl
!
! !湿潤調節法を適用
! call MoistAdjustSvapPress( xz_ExnerNl, xz_PotTempAl, xza_MixRtAl )
! call MoistAdjustNH4SH( xz_ExnerNl, xz_PotTempAl, xza_MixRtAl )
!
! !湿潤調節法による温位と混合比の変化量を保管
! call StorePotTempCond( (xz_PotTempAl - xz_PotTempWork) / DelTimeLFrog(t) )
! call StoreMixRtCond( (xza_MixRtAl - xza_MixRtWork) / DelTimeLFrog(t) )
!
! !境界条件
! call BoundaryXCyc_xza( xza_MixRtAl )
! call BoundaryZSym_xza( xza_MixRtAl )
!
!
! !-------------------------------------------------------------
! ! 暖かい雨のパラメタリゼーション.
! ! 蒸気<-->雨 の変換を行う
! !-------------------------------------------------------------
! !雨から蒸気への変換に伴う混合比の変化を計算
! ! xza_Rain2Gas 関数は, 引数として時間刻みを取ることで,
! ! 時間積分値を出力する.
!
! !これまでの値を作業配列に保管
! xz_PotTempWork = xz_PotTempAl
! xza_MixRtWork = xza_MixRtAl
!
! !雨から蒸気への混合比変化を求める
! ! 温位の計算において, 混合比変化が必要となるため,
! ! 混合比変化を 1 つの配列として用意する.
! xza_DelMixRt = 0.0d0
! xza_DelMixRt = &
! & ( &
! & + xza_Rain2Gas( &
! & xz_ExnerNl, xz_PotTempAl, xza_MixRtAl, DelTimeLFrog(t) &
! & ) &
! & + xza_Rain2GasNH4SH( &
! & xz_ExnerNl, xz_PotTempAl, xza_MixRtAl, DelTimeLFrog(t) &
! & ) &
! & )
!
! !温位の計算. 雨から蒸気への変換に伴う潜熱・反応熱を追加.
! xz_PotTempAl = &
! & xz_PotTempWork &
! & + ( &
! & + xz_Rain2GasHeat( xz_PotTempAl, xz_ExnerNl, xza_DelMixRt ) &
! & + xz_Rain2GasHeatNH4SH( xz_ExnerNl, xza_DelMixRt ) &
! & )
!
! !混合比の計算. 雨から蒸気への変換分を追加
! xza_MixRtAl = xza_MixRtWork + xza_DelMixRt
!
! !雨から蒸気の値の保管
! call StorePotTempCond( (xz_PotTempAl - xz_PotTempWork) / DelTimeLFrog(t) )
! call StoreMixRtCond( xza_DelMixRt / DelTimeLFrog(t) )
!
! !境界条件
! call BoundaryXCyc_xz( xz_PotTempAl )
! call BoundaryZSym_xz( xz_PotTempAl )
! call BoundaryXCyc_xza( xza_MixRtAl )
! call BoundaryZSym_xza( xza_MixRtAl )
!
!-------------------------------------------------------------
! 長い時間ステップでの, 速度の移流, 拡散, 数値粘性, 浮力
! (2008/06/20, 山下達也 : 温位の移流, 乱流, 拡散, 数値粘性の項
! を格納する作業配列を追加)
!-------------------------------------------------------------
xz_PotTempSum = xz_PotTempNl + xz_PotTempBasicZ
xz_TendPotTempNl = + xz_AdvScalar( xz_PotTempSum, pz_VelXNl, xr_VelZNl) + xz_TurbScalar(xz_PotTempBl, xz_KhBl) + xz_TurbScalar(xz_PotTempBasicZ, xz_KhBl) + xz_NumDiffScalar(xz_PotTempBl)
! xz_TendPotTempNl = &
! & + xz_AdvScalar( xz_PotTempNl, pz_VelXNl, xr_VelZNl) &
! & + xz_AdvScalar( xz_PotTempBasicZ, pz_VelXNl, xr_VelZNl) &
! & + xz_TurbScalar(xz_PotTempBl, xz_KhBl) &
! & + xz_TurbScalar(xz_PotTempBasicZ, xz_KhBl) &
! & + xz_NumDiffScalar(xz_PotTempBl)
! do k = DimZMin, DimZMax
! write(*,*) "TendPotTempNl",xz_TendPotTempNl(1,k)
! end do
! write(*,*) "PotTempBl(50,1)",xz_PotTempBl(50,1)
! write(*,*) "TendPotTempNl(50,1)",xz_TendPotTempNl(50,1)
! worknum = xz_NumDiffScalar(xz_PotTempBl)
! write(*,*) "xz_NumDiffScalar(50,1)",worknum(50,1)
! worknum = - xz_AdvScalar( xz_PotTempNl, pz_VelXNl, xr_VelZNl) &
! & - xz_AdvScalar( xz_PotTempBasicZ, pz_VelXNl, xr_VelZNl)
! write(*,*) "xz_AdvScalar(50,1)",worknum(50,1)
! worknum = xz_TurbScalar(xz_PotTempBl, xz_KhBl) &
! & + xz_TurbScalar(xz_PotTempBasicZ, xz_KhBl)
! write(*,*) "xz_TurbScalar(50,1)",worknum(50,1)
! write(*,*) "PotTempBl(50,2)",xz_PotTempBl(50,2)
! write(*,*) "xz_NumDiffScalar(50,2)",worknum(50,2)
! write(*,*) "PotTempBl(50,3)",xz_PotTempBl(50,3)
! write(*,*) "xz_NumDiffScalar(50,3)",worknum(50,3)
pz_AccelVelXNl = + pz_AdvVelX(pz_VelXNl, xr_VelZNl) + pz_TurbVelX(xz_KmBl, pz_VelXBl, xr_VelZBl) + pz_NumDiffVelX(pz_VelXBl)
xr_AccelVelZNl = + xr_AdvVelZ(xr_VelZNl, pz_VelXNl) + xr_Buoy(xz_PotTempNl) + xr_TurbVelZ(xz_KmBl, pz_VelXBl, xr_VelZBl) + xr_NumDiffVelZ(xr_VelZBl)
! write(*,*) "OK"
! xz_TendPotTempNl = &
! & + xz_AdvScalar( xz_PotTempNl, pz_VelXNl, xr_VelZNl) &
! & + xz_AdvScalar( xz_PotTempBasicZ, pz_VelXNl, xr_VelZNl) &
! & + xz_TurbScalar(xz_PotTempBl, xz_KhBl) &
! & + xz_TurbScalar(xz_PotTempBasicZ, xz_KhBl) &
! & + xz_NumDiffScalar(xz_PotTempBl)
!-------------------------------------------------------------
! 短い時間ステップの初期値作成
!-------------------------------------------------------------
pz_VelXNs = pz_VelXBl
xr_VelZNs = xr_VelZBl
xz_ExnerNs = xz_ExnerBl
xz_PotTempNs = xz_PotTempBl
xz_DensCloudNs = xz_DensCloudBl
xz_SatRatioNs = xz_SatRatioBl
!-------------------------------------------------------------
! 短い時間ステップの積分. オイラー法を利用
!-------------------------------------------------------------
Eular: do tau = 1, NstepEular(t)
!-------------------------------------------------------------
! 2008/04/14 山下追加
! 雲密度の計算, 凝結熱の計算, 凝結量の計算
!-------------------------------------------------------------
!=== 雲密度の移流計算
call DensityCloud( xz_DensCloudNs, DelTimeShort, pz_VelXNs, xr_VelZNs, xz_DensCloudNs, xz_DensCloudBl, xz_KhBl, xz_DensCloudWorkAs ) !(out)雲密度(負の密度調整前)
! 境界条件
call BoundaryXCyc_xz( xz_DensCloudWorkAs )
call BoundaryZSym_xz( xz_DensCloudWorkAs )
! ! 領域全体の雲の質量(fillnegative の保存性をチェック)
! CloudMassTotal01 = DelX * DelZ &
! & * sum( xz_DensCloudWorkAs(RegXMin:RegXMax,RegZMin:RegZMax) )
! write(*,*) "total cloud mass(before fillnegative)", CloudMassTotal01
! DensCloudMin = minval( xz_DensCloudWorkAs )
! write(*,*) "minimum value of DensCloud", DensCloudMin
xz_DensCloudAs = xz_FillNegativeMMC_xz( xz_DensCloudWorkAs )
! 埋めた/削った量を保管
call StoreDensCloudFill( (xz_DensCloudAs - xz_DensCloudWorkAs) / DelTimeShort )
xz_DensCloudWorkAs = xz_DensCloudAs
xz_DensCloudAs = max( 0.0d0, xz_DensCloudAs )
! やむなくゼロに戻した量を保管
call StoreDensCloudFillZero( (xz_DensCloudAs - xz_DensCloudWorkAs) / DelTimeShort )
! ! 移流によって負になった部分を埋める
! ! 負の部分が無くなるまで繰り返す
! do while( DensCloudMin < 0.0d0 )
! xz_DensCloudAs = xz_FillNegativeMMC_xz( xz_DensCloudWorkAs )
! DensCloudMin = minval( xz_DensCloudAs )
! write(*,*) "minimum value of DensCloud(in the loop)", DensCloudMin
! xz_DensCloudWorkAs = xz_DensCloudAs
! end do
! 境界条件
call BoundaryXCyc_xz( xz_DensCloudAs )
call BoundaryZSym_xz( xz_DensCloudAs )
! ! 領域全体の雲の質量(fillnegative の保存性をチェック)
! CloudMassTotal02 = DelX * DelZ &
! & * sum( xz_DensCloudAs(RegXMin:RegXMax,RegZMin:RegZMax) )
! write(*,*) "total cloud mass(after fillnegative)", CloudMassTotal02
! !=== 単位凝結量あたりの潜熱の計算
! call LatentHeatPerMass( xz_LatHeatPerMassNl )
!(out) 単位質量あたりの潜熱
! !=== 雲密度の凝結部分の計算
! call MassCondense( &
! & DelTimeShort, & !(in) 時間間隔
! & xz_LatHeatPerMassNl, & !(in) 単位質量あたりの潜熱
! & xz_SatRatioNs, & !(in) 過飽和度
! & xz_PotTempNs, & !(in) 温位
! & xz_ExnerNs, & !(in) 無次元圧力
! & xz_DensCloudAs, & !(inout) 雲密度
! & xz_MassCondNs) !(out)凝結量
xz_MassCondNs = 0.0d0
call StoreDensCloudCond( xz_MassCondNs )
! 境界条件
call BoundaryXCyc_xz( xz_DensCloudAs )
call BoundaryZSym_xz( xz_DensCloudAs )
! write(*,*) "OK!"
! (2005/12/07 北守太一)
! * 凝結により発生した潜熱を計算
!=== 凝結熱の計算
call LatentHeat( xz_MassCondNs, xz_LatHeatPerMassNl, xz_QCond) !(out)凝結熱による温位変化率
!
! 温位の移流計算.
!(2008/06/20, 山下達也 :
! 主成分凝結を議論するため, 短い時間ステップのループ内
! で計算するよう書き換えた. )
!
xz_PotTempAs = xz_PotTempNs + DelTimeEular * ( + xz_TendPotTempNl + xz_DispHeat( xz_KmBl ) + xz_Qcond )
!境界条件
call BoundaryXCyc_xz( xz_PotTempAs )
call BoundaryZSym_xz( xz_PotTempAs )
!-------------------------------------------------------------
! 速度 u の計算
!-------------------------------------------------------------
! pz_VelXAs = &
! & pz_VelXNs &
! & + DelTimeEular &
! & * ( &
! & - pz_GradPi(xz_ExnerNs, pz_VelXNs, xr_VelZNs) &
! & + pz_AccelVelXNl &
! & )
pz_VelXAs = 0.0d0
!境界条件
call BoundaryXCyc_pz( pz_VelXAs )
call BoundaryZSym_pz( pz_VelXAs )
!-------------------------------------------------------------
! エクスナー関数の計算
!-------------------------------------------------------------
! xz_ExnerAs = xz_Exner( &
! & xr_AccelVelZNl, &
! & pz_VelXNs, &
! & pz_VelXAs, &
! & xr_VelZNs, &
! & xz_ExnerNs)
! xz_ExnerAs = xz_Exner( &
! & xr_AccelVelZNl, &
! & pz_VelXNs, &
! & pz_VelXAs, &
! & xr_VelZNs, &
! & xz_ExnerNs, &
!! & xz_ExnerBasicZ, &
!! & xz_MassCondNl, &
! & xz_MassCondNs, &
! & xz_LatHeatPerMassNl)
! xz_QRad = xz_RadHeatBalance( xz_ExnerBasicZ, xz_PotTempBasicZ, xz_ExnerNs, xz_PotTempNs )
xz_QRad = 0.0d0
xz_QDis = xz_DispHeat( xz_KmBl )
xz_ExnerAs = 0.0d0
! xz_ExnerAs = xz_Exner( &
! & xr_AccelVelZNl, &
! & pz_VelXNs, &
! & pz_VelXAs, &
! & xr_VelZNs, &
! & xz_ExnerNs, &
! & xz_MassCondNs, &
! & xz_LatHeatPerMassNl, &
! & xz_QRad, &
! & xz_QDis)
!境界条件
call BoundaryXCyc_xz( xz_ExnerAs )
call BoundaryZSym_xz( xz_ExnerAs )
! write(*,*) "pz_VelXAs", pz_VelXAs(50,1)
! write(*,*) "xz_ExnerAs", xz_ExnerAs(50,1)
! write(*,*) "OK"
!-------------------------------------------------------------
! 速度 w の計算
!-------------------------------------------------------------
! xr_VelZAs = &
! & xr_VelZNs &
! & + DelTimeEular &
! & * ( &
! & - xr_GradPi(xz_ExnerAs,xz_ExnerNs,pz_VelXNs,xr_VelZNs) &
! & + xr_AccelVelZNl &
! & )
xr_VelZAs = 10.0d0
!境界条件
call BoundaryXCyc_xr( xr_VelZAs )
call BoundaryZAntiSym_xr( xr_VelZAs )
! worknum = - xr_GradPi(xz_ExnerAs,xz_ExnerNs,pz_VelXNs,xr_VelZNs)
! write(*,*) "- xr_GradPi", worknum(50,1)
! worknum = xr_AccelVelZNl
! write(*,*) "xr_AccelVelZNl", worknum(50,1)
! write(*,*) "xr_VelZAs", xr_VelZAs(50,1)
! write(*,*) "OK"
! (2008/06/16 山下達也)
!=== 過飽和度の計算
call SaturationRatio( xz_ExnerAs, xz_PotTempAs, xz_SatRatioAs) !(out)
!-------------------------------------------------------------
! 短い時間ステップのループを回すための処置
!-------------------------------------------------------------
xz_ExnerNs = xz_ExnerAs
xz_PotTempNs = xz_PotTempAs
pz_VelXNs = pz_VelXAs
xr_VelZNs = xr_VelZAs
xz_DensCloudNs = xz_DensCloudAs
xz_SatRatioNs = xz_SatRatioAs
end do Eular
! write(*,*) "xz_DensCloudCond(50,7) [sum]",xz_DensCloudCond(50,7)
! 短いループでの tendency を平均
call StoreDensCloudMean( )
! workcloud = xz_DensCloudCond(50,7) * dble( NstepEular(t) )
! write(*,*) "xz_DensCloudCond(50,7) [mean * step]", workcloud
! write(*,*) "xz_DensCloudCond(50,7) [mean]", xz_DensCloudCond(50,7)
! 平均された tendency の足しこみ
call DensCloudAdvTendSum( xz_DensCloudAdv )
call DensCloudTurbTendSum( xz_DensCloudTurb )
call DensCloudDiffTendSum( xz_DensCloudDiff )
call DensCloudCondTendSum( xz_DensCloudCond )
call DensCloudFillTendSum( xz_DensCloudFill )
call DensCloudFillZeroTendSum( xz_DensCloudFillZero )
! write(*,*) "xz_DensCloudCond2(50,7)",xz_DensCloudCond2(50,7)
! write(*,*) "xz_DensCloudCond(50,7) [after cleaning]",xz_DensCloudCond(50,7)
!----------------------------------------------------------------
! 最終的な短い時間ステップでの値を長い時間ステップでの値とみなす
!----------------------------------------------------------------
xz_ExnerAl = xz_ExnerAs
xz_PotTempAl = xz_PotTempAs
pz_VelXAl = pz_VelXAs
xr_VelZAl = xr_VelZAs
xz_DensCloudAl = xz_DensCloudAs
xz_SatRatioAl = xz_SatRatioAs
!-------------------------------------------------------------
! アッセリンの時間フィルタ. Nl の値をフィルタリング
! 1 ステップ目は Euler 法で計算するため, フィルタをかけない.
! (2008/08/25 山下達也)
!-------------------------------------------------------------
if ( t /= 1 ) then
call AsselinFilter_xz(xz_ExnerAl, xz_ExnerNl, xz_ExnerBl)
call AsselinFilter_pz(pz_VelXAl, pz_VelXNl, pz_VelXBl )
call AsselinFilter_xr(xr_VelZAl, xr_VelZNl, xr_VelZBl )
call AsselinFilter_xz(xz_PotTempAl, xz_PotTempNl, xz_PotTempBl)
call AsselinFilter_xz(xz_KmAl, xz_KmNl, xz_KmBl)
call AsselinFilter_xz(xz_DensCloudAl, xz_DensCloudNl, xz_DensCloudBl)
! call AsselinFilter_xza(xza_MixRtAl, xza_MixRtNl, xza_MixRtBl)
else
write(*,*) "OK"
end if
! call AsselinFilter_xz(xz_ExnerAl, xz_ExnerNl, xz_ExnerBl)
! call AsselinFilter_pz(pz_VelXAl, pz_VelXNl, pz_VelXBl )
! call AsselinFilter_xr(xr_VelZAl, xr_VelZNl, xr_VelZBl )
! call AsselinFilter_xz(xz_PotTempAl, xz_PotTempNl, xz_PotTempBl)
! call AsselinFilter_xz(xz_KmAl, xz_KmNl, xz_KmBl)
!! call AsselinFilter_xza(xza_MixRtAl, xza_MixRtNl, xza_MixRtBl)
!! call AsselinFilter_xz(xz_DensCloudAl, xz_DensCloudNl, xz_DensCloudBl)
!-------------------------------------------------------------
! スポンジ層
!-------------------------------------------------------------
call DampSponge_pz( pz_VelXAl, pz_VelXBl, DelTimeLFrog(t) )
call DampSponge_xr( xr_VelZAl, xr_VelZBl, DelTimeLFrog(t) )
! xz_ExnerAl = xz_DampSponge( xz_ExnerAl, xz_ExnerBl, DelTimeLFrog(t) )
! pz_VelXAl = pz_DampSponge( pz_VelXAl, pz_VelXBl, DelTimeLFrog(t) )
! xr_VelZAl = xr_DampSponge( xr_VelZAl, xr_VelZBl, DelTimeLFrog(t) )
! xz_PotTempAl = xz_DampSponge( xz_PotTempAl, xz_PotTempBl, DelTimeLFrog(t) )
! xz_KmAl = xz_DampSponge( xz_KmAl, xz_KmBl, DelTimeLFrog(t) )
!--------------------------------------------------------------
! 混合比がゼロ以下にならないための処置
!---------------------------------------------------------------
! xza_MixRtWork = xza_MixRtAl
! xza_MixRtAl = max( - xza_MixRtBasicZ, xza_MixRtWork )
! call StoreMixRtFill2( (xza_MixRtAl - xza_MixRtWork) / DelTimeLFrog(t) )
!--------------------------------------------------------------
! 解析値
!---------------------------------------------------------------
call ECCM_Stab( xz_PotTempAl, xz_ExnerAl, xza_MixRtAl )
!----------------------------------------------------------------
! 長い時間ステップのループを回すための処置
!----------------------------------------------------------------
xz_PotTempBl = xz_PotTempNl
xza_MixRtBl = xza_MixRtNl
xz_ExnerBl = xz_ExnerNl
pz_VelXBl = pz_VelXNl
xr_VelZBl = xr_VelZNl
xz_KmBl = xz_KmNl
xz_KhBl = xz_KhNl
xz_DensCloudBl = xz_DensCloudNl
xz_SatRatioBl = xz_SatRatioNl
xz_PotTempNl = xz_PotTempAl
xza_MixRtNl = xza_MixRtAl
xz_ExnerNl = xz_ExnerAl
pz_VelXNl = pz_VelXAl
xr_VelZNl = xr_VelZAl
xz_KmNl = xz_KmAl
xz_KhNl = xz_KhAl
xz_DensCloudNl = xz_DensCloudAl
xz_SatRatioNl = xz_SatRatioAl
! xz_PotTempSum = xz_PotTempNl + xz_PotTempBasicZ
xz_ExnerSum = xz_ExnerNl + xz_ExnerBasicZ
xz_TempSum = xz_PotTempSum * xz_ExnerSum
! 2008/08/16(山下 達也)
! 基本場と擾乱場の和を出力する為に追加
! 積分量の計算
! 各格子点の質量密度(擾乱成分の寄与)
xz_MassDens = (xz_ExnerBasicZ + xz_ExnerNl)**( (CpDry - GasRDry)/GasRDry ) / (xz_PotTempBasicZ + xz_PotTempNl) - (xz_ExnerBasicZ )**( (CpDry - GasRDry)/GasRDry ) / xz_PotTempBasicZ
! 領域全体の気相の質量
MassTotal = (Xmax - Xmin) * PressSfc /Grav * ( xz_ExnerBasicZ(RegXMin,RegZMin)**(CpDry/GasRDry) - xz_ExnerBasicZ(RegXMin,RegZMax)**(CpDry/GasRDry) ) + PressSfc * DelX * DelZ / GasRDry * sum( xz_MassDens(RegXMin:RegXMax,RegZMin:RegZMax))
! 温位の 2 次移流項を残すことによって生じてしまう
! 全質量への寄与
xz_MassTend = - xz_DensBasicZ * xz_TendPotTempNl / xz_PotTempBasicZ
MassTend = DelX * DelZ * sum( xz_MassTend(RegXMin:RegXMax,RegZMin:RegZMax) )
! ! 各格子点の運動エネルギー(密度を和の密度で評価)
! xz_KineticEnergy = 0.5d0 * PressSfc * DelX * DelZ / GasRDry &
! & * (xz_ExnerBasicZ + xz_ExnerNl )**((CpDry - GasRDry)/GasRDry) &
! & / (xz_PotTempBasicZ + xz_PotTempNl) &
! & * ( xz_avr_pz(pz_VelXNl)**2.0d0 + xz_avr_xr(xr_VelZNl)**2.0d0 &
! & )
! 各格子点の運動エネルギー(密度を基本場の密度で評価)
xz_KineticEnergy = 0.5d0 * PressSfc * DelX * DelZ / GasRDry * (xz_ExnerBasicZ)**((CpDry - GasRDry)/GasRDry) / (xz_PotTempBasicZ) * ( xz_avr_pz(pz_VelXNl)**2.0d0 + xz_avr_xr(xr_VelZNl)**2.0d0 )
! 領域全体の運動エネルギー
KineticEnergyTotal = sum( xz_KineticEnergy(RegXMin:RegXMax,RegZMin:RegZMax) )
! 領域全体の雲の質量
CloudMassTotal = DelX * DelZ * sum( xz_DensCloudNl(RegXMin:RegXMax,RegZMin:RegZMax) )
! 鉛直座標を 2 次元配列化
do k = RegZMin, RegZMax
xz_Z(:,k) = s_Z(k)
end do
! 各格子点の位置エネルギー
xz_PotentialEnergy = - xz_DensBasicZ * xz_PotTempNl * xz_Z / xz_PotTempBasicZ
! 領域全体の位置エネルギー
PotentialEnergyTotal = DelX * DelZ * sum( xz_PotentialEnergy(RegXMin:RegXMax,RegZMin:RegZMax) )
! 各格子点の弾性エネルギー(1 次量)
xz_ElasticEnergyFO = xz_DensBasicZ * (CpDry * xz_PotTempBasicZ / xz_VelSoundBasicZ )**2.0d0 * xz_ExnerBasicZ * xz_ExnerNl
! 領域全体の弾性エネルギー(1 次量)
ElasticEnergyFOTotal = DelX * DelZ * sum( xz_ElasticEnergyFO(RegXMin:RegXMax,RegZMin:RegZMax) )
! 各格子点の弾性エネルギー(2 次量)
xz_ElasticEnergySO = 0.5d0 * xz_DensBasicZ * (CpDry * xz_PotTempBasicZ * xz_ExnerNl / xz_VelSoundBasicZ )**2.0d0
! 領域全体の弾性エネルギー(2 次量)
ElasticEnergySOTotal = DelX * DelZ * sum( xz_ElasticEnergySO(RegXMin:RegXMax,RegZMin:RegZMax) )
end do
!----------------------------------------------------------------
! ファイル出力
!----------------------------------------------------------------
! CFL 条件のチェック
call MessageNotify( "M", "main", "Time = %f", d=(/Time/) )
call CFLCheckTimeLongVelX( pz_VelXNl )
call CFLCheckTimeLongVelZ( xr_VelZNl )
! write(*,*) "xz_DensCloudCond2(50,7) [before historyfile output]",xz_DensCloudCond2(50,7)
!ヒストリファイル出力
call HistoryFile_Output( Time, xz_PotTempNl, xz_PotTempSum, xz_TempSum, xz_ExnerNl, pz_VelXNl, xr_VelZNl, xza_MixRtNl, xz_KmNl, xz_KhNl, xz_DensCloudNl, xz_SatRatioNl, MassTotal, MassTend, KineticEnergyTotal, CloudMassTotal, PotentialEnergyTotal, ElasticEnergyFOTotal, ElasticEnergySOTotal )
!積算値のクリア
call StorePotTempClean
call StoreDensCloudClean
call StoreMixRtClean
call StoreBuoyClean
call StoreStabClean
! write(*,*) "xz_DensCloudCond2(50,7) [after cleaning]",xz_DensCloudCond2(50,7)
end do
! write(*,*) "OK"
!----------------------------------------------------------------
! 出力ファイルのクローズ
!----------------------------------------------------------------
call HistoryFile_Close
!----------------------------------------------------------------
! リスタートファイルの作成
!----------------------------------------------------------------
call ReStartFile_Open( )
call ReStartFile_OutPut( Time - DelTimeLong, xz_PotTempBl, xz_ExnerBl, pz_VelXBl, xr_VelZBl, xza_MixRtBl, xz_KmBl, xz_KhBl, xz_DensCloudBl, xz_SatRatioBl )
call ReStartFile_OutPut( Time, xz_PotTempNl, xz_ExnerNl, pz_VelXNl, xr_VelZNl, xza_MixRtNl, xz_KmNl, xz_KhNl, xz_DensCloudNl, xz_SatRatioNl )
call ReStartFile_Close
contains
subroutine ArareAlloc
!
!初期化として, 配列を定義し, 値としてゼロを代入する.
!
!暗黙の型宣言禁止
implicit none
!配列割り当て
allocate( pz_VelXBl(DimXMin:DimXMax, DimZMin:DimZMax), pz_VelXNl(DimXMin:DimXMax, DimZMin:DimZMax), pz_VelXAl(DimXMin:DimXMax, DimZMin:DimZMax), pz_VelXNs(DimXMin:DimXMax, DimZMin:DimZMax), pz_VelXAs(DimXMin:DimXMax, DimZMin:DimZMax), xr_VelZBl(DimXMin:DimXMax, DimZMin:DimZMax), xr_VelZNl(DimXMin:DimXMax, DimZMin:DimZMax), xr_VelZAl(DimXMin:DimXMax, DimZMin:DimZMax), xr_VelZNs(DimXMin:DimXMax, DimZMin:DimZMax), xr_VelZAs(DimXMin:DimXMax, DimZMin:DimZMax), xz_ExnerBl(DimXMin:DimXMax, DimZMin:DimZMax), xz_ExnerNl(DimXMin:DimXMax, DimZMin:DimZMax), xz_ExnerNs(DimXMin:DimXMax, DimZMin:DimZMax), xz_ExnerAl(DimXMin:DimXMax, DimZMin:DimZMax), xz_ExnerAs(DimXMin:DimXMax, DimZMin:DimZMax), xz_ExnerSum(DimXMin:DimXMax, DimZMin:DimZMax), xz_PotTempWork(DimXMin:DimXMax, DimZMin:DimZMax), xz_PotTempSum(DimXMin:DimXMax, DimZMin:DimZMax), xz_PotTempBl(DimXMin:DimXMax, DimZMin:DimZMax), xz_PotTempNl(DimXMin:DimXMax, DimZMin:DimZMax), xz_PotTempAl(DimXMin:DimXMax, DimZMin:DimZMax), xz_PotTempNs(DimXMin:DimXMax, DimZMin:DimZMax), xz_PotTempAs(DimXMin:DimXMax, DimZMin:DimZMax), xz_TempSum(DimXMin:DimXMax, DimZMin:DimZMax), xz_KmBl(DimXMin:DimXMax, DimZMin:DimZMax), xz_KmNl(DimXMin:DimXMax, DimZMin:DimZMax), xz_KmAl(DimXMin:DimXMax, DimZMin:DimZMax), xz_KhBl(DimXMin:DimXMax, DimZMin:DimZMax), xz_KhNl(DimXMin:DimXMax, DimZMin:DimZMax), xz_KhAl(DimXMin:DimXMax, DimZMin:DimZMax), xza_MixRtWork(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum), xza_MixRtBl(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum), xza_MixRtNl(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum), xza_MixRtAl(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum), pz_AccelVelXNl(DimXMin:DimXMax, DimZMin:DimZMax), xr_AccelVelZNl(DimXMin:DimXMax, DimZMin:DimZMax), xz_TendPotTempNl(DimXMin:DimXMax, DimZMin:DimZMax), xza_DelMixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum), DelTimeLFrog(NstepLong), NStepEular(NStepLong), xz_LatHeatPerMassNl(DimXMin:DimXMax, DimZMin:DimZMax), xz_MassCondNs(DimXMin:DimXMax, DimZMin:DimZMax), xz_MassCondNl(DimXMin:DimXMax, DimZMin:DimZMax), xz_DensCloudBl(DimXMin:DimXMax, DimZMin:DimZMax), xz_DensCloudNl(DimXMin:DimXMax, DimZMin:DimZMax), xz_DensCloudAl(DimXMin:DimXMax, DimZMin:DimZMax), xz_DensCloudNs(DimXMin:DimXMax, DimZMin:DimZMax), xz_DensCloudAs(DimXMin:DimXMax, DimZMin:DimZMax), xz_DensCloudWorkAs(DimXMin:DimXMax, DimZMin:DimZMax), xz_SatRatioBl(DimXMin:DimXMax, DimZMin:DimZMax), xz_SatRatioNl(DimXMin:DimXMax, DimZMin:DimZMax), xz_SatRatioAl(DimXMin:DimXMax, DimZMin:DimZMax), xz_SatRatioNs(DimXMin:DimXMax, DimZMin:DimZMax), xz_SatRatioAs(DimXMin:DimXMax, DimZMin:DimZMax), xz_QCond(DimXMin:DimXMax, DimZMin:DimZMax), xz_QRad(DimXMin:DimXMax, DimZMin:DimZMax), xz_QDis(DimXMin:DimXMax, DimZMin:DimZMax), worknum(DimXMin:DimXMax, DimZMin:DimZMax), xz_MassDens(DimXMin:DimXMax, DimZMin:DimZMax), xz_MassTend(DimXMin:DimXMax, DimZMin:DimZMax), xz_KineticEnergy(DimXMin:DimXMax, DimZMin:DimZMax), xz_PotentialEnergy(DimXMin:DimXMax, DimZMin:DimZMax), xz_ElasticEnergyFO(DimXMin:DimXMax, DimZMin:DimZMax), xz_ElasticEnergySO(DimXMin:DimXMax, DimZMin:DimZMax), xz_Z(DimXMin:DimXMax, DimZMin:DimZMax) )
pz_VelXNs = 0.0d0
pz_VelXAs = 0.0d0
xr_VelZNs = 0.0d0
xr_VelZAs = 0.0d0
xz_ExnerNs = 0.0d0
xz_ExnerAs = 0.0d0
xz_PotTempNs = 0.0d0
xz_PotTempAs = 0.0d0
pz_VelXBl = 0.0d0
pz_VelXNl = 0.0d0
pz_VelXAl = 0.0d0
xr_VelZBl = 0.0d0
xr_VelZNl = 0.0d0
xr_VelZAl = 0.0d0
xz_ExnerBl = 0.0d0
xz_ExnerNl = 0.0d0
xz_ExnerAl = 0.0d0
xz_KmBl = 0.0d0
xz_KmNl = 0.0d0
xz_KmAl = 0.0d0
xz_KhBl = 0.0d0
xz_KhNl = 0.0d0
xz_KhAl = 0.0d0
xz_PotTempBl = 0.0d0
xz_PotTempNl = 0.0d0
xz_PotTempAl = 0.0d0
xza_MixRtBl = 0.0d0
xza_MixRtNl = 0.0d0
xza_MixRtAl = 0.0d0
pz_AccelVelXNl = 0.0d0
xr_AccelVelZNl = 0.0d0
xza_DelMixRt = 0.0d0
xz_LatHeatPerMassNl = 0.0d0
xz_MassCondNs = 0.0d0
xz_MassCondNl = 0.0d0
xz_DensCloudBl = 0.0d0
xz_DensCloudNl = 0.0d0
xz_DensCloudAl = 0.0d0
xz_DensCloudNs = 0.0d0
xz_DensCloudAs = 0.0d0
xz_DensCloudWorkAs = 0.0d0
xz_SatRatioBl = 0.0d0
xz_SatRatioNl = 0.0d0
xz_SatRatioAl = 0.0d0
xz_SatRatioNs = 0.0d0
xz_SatRatioAs = 0.0d0
xz_QCond = 0.0d0
xz_QRad = 0.0d0
xz_QDis = 0.0d0
DelTimeEular = 0.0d0
DelTimeLFrog = 0.0d0
NStepEular = 0.0d0
NstepLFrog = 0.0d0
worknum = 0.0d0
xz_MassDens = 0.0d0
xz_MassTend = 0.0d0
xz_KineticEnergy = 0.0d0
xz_PotentialEnergy = 0.0d0
xz_ElasticEnergyFO = 0.0d0
xz_ElasticEnergySO = 0.0d0
xz_Z = 0.0d0
end subroutine ArareAlloc
end program arare