!= Program Arare ! ! Authors:: SUGIYAMA Ko-ichiro, ODAKA Masatsugu ! Version:: $Id: arare_mars1.f90,v 1.1 2011-03-30 13:53:07 sugiyama Exp $ ! Tag Name:: $Name: arare4-20120911 $ ! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved. ! License:: See COPYRIGHT[link:../../COPYRIGHT] ! !== Overview ! ! 非静力学モデル deepconv/arare. ! 火星湿潤対流用. ! ! !== Error Handling ! !== Known Bugs ! !== Note ! ! * 方程式系は準圧縮系. ! * 温位の移流部分をひとまとめにして計算(08/12/02) ! * 位置エネルギー, 弾性エネルギーを出力(08/12/02) ! * 雲密度の計算に乱流拡散, 数値粘性を導入(08/12/19) ! * 圧力方程式に散逸加熱・放射加熱の効果を考慮(09/10/12) ! * 質量収支チェック用プログラム. ! * 雲密度の計算には ! densitycloud_turb-masstest2.f90 (移流拡散だけ解く) ! を用いること. ! * データ出力には historyfileio_mmconv2_masstest3.f90 を用いること. ! * 主成分凝結計算用穴埋めプログラムを導入. ! * 解析量の初期値を代入するよう書き換え. ! * 雲密度の 2 次元 tendency の計算手順を変更. ! * 短い時間ステップ内で和を取り, その平均量を保管 ! * 保管した値を出力時間間隔内で足しこむ ! * 次の出力区間に移る前に値をクリア ! * 雲密度の移流を長い時間ステップで評価 ! * やむなくゼロに引き戻した雲密度を出力するよう変更 ! * 雲密度の計算に時間フィルターを適用 ! * 凝結率の気相質量を全量で評価 ! !== Future Plans ! ! 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_mmc, 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 StoreMom, only: StoreMom_init, StoreMomClean ! 湿潤ルーチン設定モジュール 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_mmc, only : ReStartFile_Open, ReStartFile_OutPut, & & ReStartFile_Close, ReStartFile_Get ! ヒストリファイル入出力モジュール use HistoryFileIO_mmc, 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 DynImpFunc_mmc, 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( ) call StoreMom_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_mmc() call DisturbEnv_mmc(cfgfile, & & xz_PotTempBl, xz_ExnerBl, pz_VelXBl, xr_VelZBl, & & 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_BuoyMolWt(xza_MixRtNl) & ! & + xr_BuoyDrag(xza_MixRtNl) & & + 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, & !(in) 現在の時間の雲密度 & DelTimeShort, & !(in) 時間間隔 & pz_VelXNl, xr_VelZNl, xz_DensCloudNl, & !(in) フラックス項の計算で用いる値 & xz_DensCloudBl, & !(in) 乱流拡散, 数値粘性で用いる値 & xz_KhBl, & !(in) 乱流拡散係数 & 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)凝結量 ! 境界条件 call BoundaryXCyc_xz( xz_DensCloudAs ) call BoundaryZSym_xz( xz_DensCloudAs ) ! write(*,*) "OK!" ! (2005/12/07 北守太一) ! * 凝結により発生した潜熱を計算 !=== 凝結熱の計算 call LatentHeat( & & xz_MassCondNs, & !(in) 凝結量 & xz_LatHeatPerMassNl, & !(in) 単位質量あたりの潜熱 & xz_ExnerNs, & !(in) エクスナー関数擾乱成分 & xz_PotTempNs, & !(in) 温位擾乱成分 & xz_QCond) !(out)凝結熱による温位変化率 ! ! 温位の移流計算. !(2008/06/20, 山下達也 : ! 主成分凝結を議論するため, 短い時間ステップのループ内 ! で計算するよう書き換えた. ) ! xz_PotTempAs = & & xz_PotTempNs & & + DelTimeEular & & * ( & ! & + 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_TendPotTempNl & & + xz_DispHeat( xz_KmBl ) & ! & + xz_RadHeatConst( xz_ExnerBasicZ ) & & + xz_RadHeatBalance( xz_ExnerBasicZ, xz_PotTempBasicZ, xz_ExnerNs, xz_PotTempNs ) & ! & + xz_HeatFluxDiff( xz_PotTempNl ) & & + xz_Qcond & ! & + xz_HeatFluxBulk( xz_PotTempNl, pz_VelXNl ) & ! & + xz_HeatFluxBulk( xz_PotTempNl ) & !! & + xz_NewtonCool( xz_PotTempBl ) & & ) !境界条件 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 & & ) !境界条件 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_QDis = xz_DispHeat( xz_KmBl ) 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 & & ) !境界条件 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, & !(in) & xz_PotTempAs, & !(in) & 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 call StoreMomClean ! 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