!= Program Arare ! ! Authors:: SUGIYAMA Ko-ichiro, ODAKA Masatsugu ! Version:: $Id: arare.f90,v 1.47 2007/08/24 07:51:03 sugiyama Exp $ ! Tag Name:: $Name: arare4-20071012 $ ! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved. ! License:: See COPYRIGHT[link:../../COPYRIGHT] ! !== Overview ! ! 非静力学モデル deepconv/arare. ! !== Error Handling ! !== Known Bugs ! !== Note ! ! * 方程式系は準圧縮系. ! !== 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, 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 basicset, only : basicset_init, & & xz_DensBasicZ, xza_MixRtBasicZ, xz_PotTempBasicZ, & & xz_VelSoundBasicZ ! 積算値管理モジュール use StorePotTemp, only : StorePotTemp_init, StorePotTempClean, & & StorePotTempCond, z_Flux 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 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 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 , eddyviscosity ! 放射強制計算用モジュール use Radiation, only : Radiation_init, & & xz_RadHeatConst ! 地表フラックス計算用モジュール use HeatFlux_N1994, only : xz_HeatFluxBulk, xz_MixRtFluxBulk ! & pz_MomFluxBulk ! 湿潤飽和調節法計算用モジュール use MoistAdjust, only : MoistAdjustSvapPress, MoistAdjustNH4SH ! 雲物理パラメタリゼーション use WarmRainPrm, only : WarmRainPrm_Init, xz_Rain2GasHeat, xza_Rain2Gas, & & xza_Rain2GasNH4SH, xz_Rain2GasHeatNH4SH, & & xza_Cloud2Rain, xza_FallRain ! 湿潤気塊の浮力計算用モジュール 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_PotTempWork(:,:) real(8), allocatable :: xz_PotTempBl(:,:) real(8), allocatable :: xz_PotTempNl(:,:) real(8), allocatable :: xz_PotTempAl(:,:) 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(:,:,:) real(8) :: Time real(8) :: ReStartTime(2) real(8), allocatable :: DelTimeLFrog(:) real(8) :: DelTimeEular integer, allocatable :: NStepEular(:) integer :: NStepLFrog integer :: t, tau, t1, t2 integer :: kz !コマンドライン引数の解釈 ! 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 moistset_init() write(*,*) "OK" !積算値を保管するためのモジュールの初期化 ! NAMELIST から情報を得て, 基本場を設定する. call StorePotTemp_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_PotTempNl, xz_ExnerNl, pz_VelXNl, xr_VelZNl, & & xza_MixRtNl, xz_KmNl, xz_KhNl ) write(*,*) "OK" else !デフォルト設定の基本場, 擾乱場を作成する. call BasicEnv() call DisturbEnv(cfgfile, & & xz_PotTempBl, xz_ExnerBl, pz_VelXBl, xr_VelZBl, & & xza_MixRtBl, xz_KmBl, xz_KhBl ) ! 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 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_ExnerNl, & & pz_VelXNl, & & xr_VelZNl, & & xza_MixRtNl, & & xz_KmNl, & & xz_KhNl ) end if !時刻 0 の場合には出力 if ( Time == 0 ) then call HistoryFile_Output( & & Time, & & xz_PotTempNl, & & xz_ExnerNl, & & pz_VelXNl, & & xr_VelZNl, & & xza_MixRtNl, & & xz_KmNl, & & xz_KhNl ) end if !---------------------------------------------------------------------- ! 設定のチェック !---------------------------------------------------------------------- !CFL 条件のチェック call CFLCheckTimeShort( xz_VelSoundBasicZ ) !---------------------------------------------------------------------- ! 数値積分 !---------------------------------------------------------------------- do t1 = 1, NstepLFrog / NstepDisp do t2 = 1, NstepDisp t = (t1 - 1) * NstepDisp + t2 !時刻の設定 Time = Time + DelTimeLong !---------------------------------------------------------------- ! 渦粘性係数, 渦拡散係数を求める. !---------------------------------------------------------------- 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) & & ) ! call EddyViscosity(pz_VelXNl, xr_VelZNl, xz_PotTempNl, xz_KmNl, xz_KmAl) !値の上限下限の設定 ! * 値は正になることを保証する ! * 値の上限は 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 !---------------------------------------------------------------- ! 温位の移流計算. !---------------------------------------------------------------- !時間積分 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_PotTempBl, pz_VelXBl ) & !! & + 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 ) !------------------------------------------------------------- ! 長い時間ステップでの, 速度の移流, 拡散, 数値粘性, 浮力 !------------------------------------------------------------- pz_AccelVelXNl = & & + pz_AdvVelX(pz_VelXNl, xr_VelZNl) & & + pz_TurbVelX(xz_KmBl, pz_VelXBl, xr_VelZBl) & !! & + pz_MomFluxBulk(pz_VelXNl) & & + 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) !------------------------------------------------------------- ! 短い時間ステップの初期値作成 !------------------------------------------------------------- pz_VelXNs = pz_VelXBl xr_VelZNs = xr_VelZBl xz_ExnerNs = xz_ExnerBl !------------------------------------------------------------- ! 短い時間ステップの積分. オイラー法を利用 !------------------------------------------------------------- Eular: do tau = 1, NstepEular(t) !------------------------------------------------------------- ! 速度 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) !境界条件 call BoundaryXCyc_xz( xz_ExnerAs ) call BoundaryZSym_xz( xz_ExnerAs ) !------------------------------------------------------------- ! 速度 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 ) !------------------------------------------------------------- ! 短い時間ステップのループを回すための処置 !------------------------------------------------------------- xz_ExnerNs = xz_ExnerAs pz_VelXNs = pz_VelXAs xr_VelZNs = xr_VelZAs end do Eular !---------------------------------------------------------------- ! 最終的な短い時間ステップでの値を長い時間ステップでの値とみなす !---------------------------------------------------------------- xz_ExnerAl = xz_ExnerAs pz_VelXAl = pz_VelXAs xr_VelZAl = xr_VelZAs !------------------------------------------------------------- ! アッセリンの時間フィルタ. Nl の値をフィルタリング !------------------------------------------------------------- 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 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_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 end do !---------------------------------------------------------------- ! ファイル出力 !---------------------------------------------------------------- ! CFL 条件のチェック call MessageNotify( "M", "main", "Time = %f", d=(/Time/) ) call CFLCheckTimeLongVelX( pz_VelXNl ) call CFLCheckTimeLongVelZ( xr_VelZNl ) !ヒストリファイル出力 call HistoryFile_Output( & & Time, & & xz_PotTempNl, & & xz_ExnerNl, & & pz_VelXNl, & & xr_VelZNl, & & xza_MixRtNl, & & xz_KmNl, & & xz_KhNl & & ) ! do kz = DimZmin, DimZMax ! write(*,*) kz, z_Flux(kz) ! end do !積算値のクリア call StorePotTempClean call StoreMixRtClean call StoreBuoyClean call StoreStabClean end do !---------------------------------------------------------------- ! 出力ファイルのクローズ !---------------------------------------------------------------- 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 ) call ReStartFile_OutPut( & & Time, & & xz_PotTempNl, xz_ExnerNl, pz_VelXNl, xr_VelZNl, & & xza_MixRtNl, xz_KmNl, xz_KhNl ) 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_ExnerAl(DimXMin:DimXMax, DimZMin:DimZMax), & & xz_ExnerNs(DimXMin:DimXMax, DimZMin:DimZMax), & & xz_ExnerAs(DimXMin:DimXMax, DimZMin:DimZMax), & ! & xz_PotTempWork(DimXMin:DimXMax, DimZMin:DimZMax), & & xz_PotTempBl(DimXMin:DimXMax, DimZMin:DimZMax), & & xz_PotTempNl(DimXMin:DimXMax, DimZMin:DimZMax), & & xz_PotTempAl(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), & ! & xza_DelMixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum), & ! & DelTimeLFrog(NstepLong), NStepEular(NStepLong) ) pz_VelXNs = 0.0d0; pz_VelXAs = 0.0d0 xr_VelZNs = 0.0d0; xr_VelZAs = 0.0d0 xz_ExnerNs = 0.0d0; xz_ExnerAs = 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 DelTimeEular = 0.0d0 DelTimeLFrog = 0.0d0; NStepEular = 0.0d0 NstepLFrog = 0.0d0 end subroutine ArareAlloc end program arare