!= Program Arare ! ! Authors:: SUGIYAMA Ko-ichiro, ODAKA Masatsugu ! Version:: $Id: arare.f90,v 1.42 2006/11/27 02:18:45 sugiyama Exp $ ! Tag Name:: $Name: arare4-20061224 $ ! Copyright:: Copyright (C) GFD Dennou Club, 2006. All rights reserved. ! License:: See COPYRIGHT[link:../../COPYRIGHT] ! !== Overview ! ! 非静力学モデル deepconv/arare. ! !== Error Handling ! !== Known Bugs ! !== Note ! ! * 方程式系は準圧縮系. ! * 湿潤過程を除いた乾燥大気バージョン (2007/01/06) ! * Skamarock and Klemp 1989 の設定 ! * 拡散係数は 75 m^2/src 一定 ! !== Future Plans ! ! program arare ! !非静力学モデル deepconv/arare. ! !----- モジュール読み込み ------ !----- 型宣言, 文字列処理 ---- use dc_types, only : STRING use dc_string, only : StoA !----- メッセージ出力 ----- use dc_message, only: MessageNotify !----- コマンドライン引数解析 ----- use dc_args, only : ARGS, Open, Debug, Help, Strict, Close, Option !----- 管理モジュール ----- ! 化学量計算モジュール use chemcalc ! 入出力ファイル名管理モジュール use fileset, only : fileset_init, & & InitFile ! デバッグ出力管理モジュール use debugset, only : debugset_init ! 時間管理モジュール use timeset, only : timeset_init, & & NstepLong, NstepShort, DelTimeLong, DelTimeShort, & & TimeDisp ! 格子点管理モジュール use gridset, only : gridset_init, & & DimXMin, DimXMax, DimZMin, DimZMax, SpcNum ! 基本場設定モジュール use basicset, only : basicset_init, & & xz_DensBasicZ, xza_MixRtBasicZ, xz_PotTempBasicZ, & & xz_VelSoundBasicZ ! 積算値管理モジュール use storeset, only : storeset_init, storeclean, & & StoreStab, StoreMixRt, StoreCond use storeset2, only : store2set_init, store2clean, & & Store2Cond, Store2Flux, Store2Fill1, Store2Fill2 !----- 下請けモジュール ----- ! 数値摩擦計算モジュール use damping, only : damping_init, & & xz_DampSponge, pz_DampSponge, xr_DampSponge ! 時間積分フィルターモジュール use timefilter, only : AsselinFilter, xza_AsselinFilter_xza_xza_xza ! 境界条件適用モジュール use boundary, only : xz_BoundaryXCyc_xz, xz_BoundaryZSym_xz, & & xza_BoundaryXCyc_xza, xza_BoundaryZSym_xza, & & pz_BoundaryXCyc_pz, pz_BoundaryZSym_pz, & & xr_BoundaryXCyc_xr, 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_AdvScalar2, 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_NumDiffScalar2, & & 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, xz_BuoyKm ! 放射強制計算用モジュール use Radiation, only : Radiation_init, & & xz_RadHeatConst, xz_NewtonCool ! 地表フラックス計算用モジュール ! use HeatFlux, only : xz_HeatFluxBulk, xz_MixRtFluxBulk use HeatFlux, only : xz_HeatFluxDiff, xza_MixRtFluxDiff ! 湿潤飽和調節法計算用モジュール use MoistAdjust, only : MoistAdjust_Init, MoistAdjustSvapPress, MoistAdjustNH4SH ! 雲物理パラメタリゼーション use WarmRainPrm, only : WarmRainPrm_Init, xz_Rain2GasHeat, xza_Rain2Gas, & & xza_Rain2GasNH4SH, xz_Rain2GasHeatNH4SH, & & xza_Cloud2Rain, xz_FallRain, xza_FallRain ! 湿潤気塊の浮力計算用モジュール use MoistBuoyancy,only : MoistBuoy_Init, xz_BuoyMoistKm, xr_BuoyMolWt, & & xr_BuoyDrag ! 断熱上昇気塊の温度減率計算用モジュール use ECCM, only : eccm_init, 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), allocatable :: xz_Stab(:,:) real(8), allocatable :: xz_StabTemp(:,:) real(8), allocatable :: xz_StabMolWt(:,:) real(8) :: Time real(8) :: ReStartTime(2) real(8), allocatable :: DelTimeLFrog(:) real(8) :: DelTimeEular integer, allocatable :: NStepEular(:) integer :: NStepLFrog integer :: t, tau !----- コマンドライン行き数解析用変数 ----- type(ARGS) :: arg ! コマンドライン引数情報 logical :: OPT_namelist ! コマンドライン引数用論理変数 character(STRING) :: VAL_namelist ! コマンドライン引数の値 !---------------------------------------------------------------------- ! 変数参照型モジュールの初期化 !---------------------------------------------------------------------- !NAMELIST ファイルの取得 ! gt4f90io ライブラリの dc_args モジュールを利用. ! 指定可能なオプションは以下の通り. ! ! -N=VAL, --namelist=VAL ! specify Namelist file (default is 'arare.conf'). ! ! -D=VAL, --debug=VAL ! call dc_trace#SetDebug (display a lot of messages for debug). ! VAL is unit number (default is standard output) ! ! -h=VAL, -H=VAL, --help=VAL ! display this help and exit. VAL is unit number (default is ! standard output) ! call Open(arg) call Option(arg, StoA('-N', '--namelist'), & & OPT_namelist, VAL_namelist, & & help="specify Namelist file (default is 'arare.conf')." ) ! "-N/--namelist" オプションの設定 call Debug(arg) ! デバッグオプションの自動設定 call Help(arg) ! ヘルプオプションの自動設定 call Strict(arg) ! 無効なオプション指定時に警告を表示 !"-N/-namelist" オプションの解釈 ! 与えられていない場合はデフォルト値 (arare.conf) を ! NAMLIST ファイル名とする. ! if (OPT_namelist) then call MessageNotify( "M", "main", & & "Namelist file is '%c'", c1=trim(VAL_namelist) ) cfgfile=trim(VAL_namelist) else call MessageNotify( "W", "main", & & "Namelist file is not specified." ) call MessageNotify( "M", "main", & & "Use default Namelist file (arare.conf)." ) cfgfile="arare.conf" end if call Close(arg) !時刻に関する設定の初期化 ! NAMELIST から必要な情報を読み取り, 時間関連の変数の設定を行う. call timeset_init(cfgfile) !格子点情報の初期化 ! NAMELIST から情報を得て, 格子点を計算する call gridset_init(cfgfile) !基本場の情報の初期化 ! NAMELIST から情報を得て, 基本場を設定する. call basicset_init(cfgfile) !I/O ファイル名の初期化 ! NAMELIST ファイル名を指定し, deepconv/arare の ! 出力ファイル名を NAMELIST から得る call fileset_init(cfgfile) !積算値を保管するためのモジュールの初期化 ! NAMELIST から情報を得て, 基本場を設定する. call storeset_init( ) call store2set_init( ) !内部変数の初期化. とりあえずゼロを入れて値を確定させておく. call ArareAlloc !予報変数の初期化 ! ReStartFile が設定されている場合にはファイルを読み込み, ! 設定されていない場合には, デフォルトの基本場と擾乱場を作る. call ECCM_Init() 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 ) 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 end if !---------------------------------------------------------------------- ! パッケージ型モジュールの初期化 ! デフォルトの値から変更する必要のあるルーチンのみ初期化 !---------------------------------------------------------------------- call Damping_Init( cfgfile ) !波の減衰係数の初期化 call NumDiffusion_Init() !数値拡散項の初期化 call Turbulence_Init() !乱流計算の初期化 call MoistAdjust_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 ) !---------------------------------------------------------------------- ! 数値積分 !---------------------------------------------------------------------- LFrog: do t = 1, NstepLFrog !時刻の設定 Time = Time + DelTimeLong ! write(*,*) Time !---------------------------------------------------------------- ! 渦粘性係数, 渦拡散係数を求める. !---------------------------------------------------------------- xz_KmAl = & & xz_KmBl & & + DelTimeLFrog(t) & & * ( & & + xz_AdvScalar2(xz_KmNl, pz_VelXNl, xr_VelZNl) & ! & + xz_BuoyMoistKm(xz_PotTempBl, xz_ExnerBl, xza_MixRtBl) & & + xz_BuoyKm(xz_PotTempBl) & & + xz_ShearKm(xz_KmBl, pz_VelXBl, xr_VelZBl) & & + xz_NumDiffScalar2(xz_KmBl) & & + xz_DispKm(xz_KmBl) & & ) !値の上限下限の設定 ! * 値は正になることを保証する ! * 値の上限は 800 とする. 中島健介(1994, 学位論文)参照 ! xz_KmAl = max( 0.0d0, min( xz_KmAl, 800.0d0 ) ) xz_KmAl = max( 0.0d0, min( xz_KmAl, 75.0d0 ) ) !境界条件 xz_KmAl = xz_BoundaryXCyc_xz( xz_KmAl ) xz_KmAl = xz_BoundaryZSym_xz( xz_KmAl ) !渦拡散定数を決める ! xz_KhAl = 3.0d0 * xz_KmAl xz_KhAl = 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_PotTempNl ) & !! & + xz_NewtonCool( xz_PotTempBl ) & & ) !境界条件 xz_PotTempAl = xz_BoundaryXCyc_xz( xz_PotTempAl ) xz_PotTempAl = xz_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_MixRtAl ) !埋めた/削った量を保管 call Store2Fill1( (xza_MixRtAl - xza_MixRtWork) / DelTimeLFrog(t) ) !境界条件 xza_MixRtAl = xza_BoundaryXCyc_xza( xza_MixRtAl ) xza_MixRtAl = xza_BoundaryZSym_xza( xza_MixRtAl ) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! 乾燥大気バージョンのため以下の部分はコメントアウト !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! !------------------------------------------------------------- ! ! 暖かい雨のパラメタリゼーション. ! ! 雲<-->雨 の変換を行う ! !------------------------------------------------------------- ! !雲から雨への変換分を追加する. ! ! xza_Cloud2Rain 関数は, 引数として時間刻みを取ることで, ! ! 時間積分値を出力する. ! ! !これまでの値を作業配列に保管 ! xza_MixRtWork = xza_MixRtAl ! ! !雨への変化量を計算 ! xza_MixRtAl = xza_MixRtWork & ! & + xza_Cloud2Rain( xza_MixRtAl, DelTimeLFrog(t) ) ! ! !雲から雨への変換量を保管 ! call Store2Cond( (xza_MixRtAl - xza_MixRtWork) / DelTimeLFrog(t) ) ! ! !境界条件 ! xza_MixRtAl = xza_BoundaryXCyc_xza( xza_MixRtAl ) ! xza_MixRtAl = xza_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 StoreCond( (xz_PotTempAl - xz_PotTempWork) / DelTimeLFrog(t) ) ! call Store2Cond( (xza_MixRtAl - xza_MixRtWork) / DelTimeLFrog(t) ) ! ! !境界条件 ! xza_MixRtAl = xza_BoundaryXCyc_xza( xza_MixRtAl ) ! xza_MixRtAl = xza_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 StoreCond( (xz_PotTempAl - xz_PotTempWork) / DelTimeLFrog(t) ) ! call Store2Cond( xza_DelMixRt / DelTimeLFrog(t) ) ! ! !境界条件 ! xz_PotTempAl = xz_BoundaryXCyc_xz( xz_PotTempAl ) ! xz_PotTempAl = xz_BoundaryZSym_xz( xz_PotTempAl ) ! xza_MixRtAl = xza_BoundaryXCyc_xza( xza_MixRtAl ) ! xza_MixRtAl = xza_BoundaryZSym_xza( xza_MixRtAl ) ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !------------------------------------------------------------- ! 長い時間ステップでの, 速度の移流, 拡散, 数値粘性, 浮力 !------------------------------------------------------------- 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) !------------------------------------------------------------- ! 短い時間ステップの初期値作成 !------------------------------------------------------------- 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 & & ) !境界条件 pz_VelXAs = pz_BoundaryXCyc_pz( pz_VelXAs ) pz_VelXAs = pz_BoundaryZSym_pz( pz_VelXAs ) !------------------------------------------------------------- ! エクスナー関数の計算 !------------------------------------------------------------- xz_ExnerAs = xz_Exner( & & xr_AccelVelZNl, & & pz_VelXNs, & & pz_VelXAs, & & xr_VelZNs, & & xz_ExnerNs) !境界条件 xz_ExnerAs = xz_BoundaryXCyc_xz( xz_ExnerAs ) xz_ExnerAs = xz_BoundaryZSym_xz( xz_ExnerAs ) !------------------------------------------------------------- ! 速度 w の計算 !------------------------------------------------------------- xr_VelZAs = & & xr_VelZNs & & + DelTimeEular & & * ( & & - xr_GradPi(xz_ExnerAs,xz_ExnerNs,pz_VelXNs,xr_VelZNs) & & + xr_AccelVelZNl & & ) !境界条件 xr_VelZAs = xr_BoundaryXCyc_xr( xr_VelZAs ) xr_VelZAs = xr_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 の値をフィルタリング !------------------------------------------------------------- xz_ExnerNl = AsselinFilter(xz_ExnerAl, xz_ExnerNl, xz_ExnerBl) pz_VelXNl = AsselinFilter(pz_VelXAl, pz_VelXNl, pz_VelXBl ) xr_VelZNl = AsselinFilter(xr_VelZAl, xr_VelZNl, xr_VelZBl ) xz_PotTempNl = AsselinFilter(xz_PotTempAl, xz_PotTempNl, xz_PotTempBl) xz_KmNl = AsselinFilter(xz_KmAl, xz_KmNl, xz_KmBl) xza_MixRtNl = xza_AsselinFilter_xza_xza_xza( & & xza_MixRtAl, xza_MixRtNl, xza_MixRtBl & & ) !------------------------------------------------------------- ! スポンジ層 !------------------------------------------------------------- ! 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) ) !! xz_ExnerAl = xz_ExnerAl - DelTimeLFrog(t) * xz_DampSponge_xz(xz_ExnerBl) !! pz_VelXAl = pz_VelXAl - DelTimeLFrog(t) * pz_DampSponge_pz(pz_VelXBl ) !! xr_VelZAl = xr_VelZAl - DelTimeLFrog(t) * xr_DampSponge_xr(xr_VelZBl ) !! xz_PotTempAl = xz_PotTempAl - DelTimeLFrog(t) * xz_DampSponge_xz(xz_PotTempBl) !! xz_KmAl = xz_KmAl - DelTimeLFrog(t) * xz_DampSponge_xz(xz_KmBl) !-------------------------------------------------------------- ! 混合比がゼロ以下にならないための処置 !--------------------------------------------------------------- xza_MixRtWork = xza_MixRtAl xza_MixRtAl = max( - xza_MixRtBasicZ, xza_MixRtAl ) call Store2Fill2( (xza_MixRtAl - xza_MixRtWork) / DelTimeLFrog(t) ) !------------------------------------------------------------- ! 境界条件 ! ここで境界条件を適用する理由がないのでコメントアウト !------------------------------------------------------------- ! xz_KmAl = xz_BoundaryXCyc_xz( xz_KmAl ) ! xz_KmAl = xz_BoundaryZSym_xz( xz_KmAl ) ! xz_PotTempAl = xz_BoundaryXCyc_xz( xz_PotTempAl ) ! xz_PotTempAl = xz_BoundaryZSym_xz( xz_PotTempAl ) ! xza_MixRtAl = xza_BoundaryXCyc_xza( xza_MixRtAl ) ! xza_MixRtAl = xza_BoundaryZSym_xza( xza_MixRtAl ) ! pz_VelXAs = pz_BoundaryXCyc_pz( pz_VelXAl ) ! pz_VelXAs = pz_BoundaryZSym_pz( pz_VelXAl ) ! xz_ExnerAs = xz_BoundaryXCyc_xz( xz_ExnerAl ) ! xz_ExnerAs = xz_BoundaryZSym_xz( xz_ExnerAl ) ! xr_VelZAs = xr_BoundaryXCyc_xr( xr_VelZAl ) ! xr_VelZAs = xr_BoundaryZAntiSym_xr( xr_VelZAl ) !---------------------------------------------------------------- ! 安定度の計算. 値は StoreSet で保持 !---------------------------------------------------------------- call ECCM_Stab( xz_PotTempAl, xz_ExnerAl, xza_MixRtAl, & & xz_Stab, xz_StabTemp, xz_StabMolWt ) call StoreStab( xz_Stab, xz_StabTemp, xz_StabMolWt ) call StoreMixRt( 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 !---------------------------------------------------------------- ! ファイル出力 !---------------------------------------------------------------- if ( mod( anint(Time), TimeDisp ) == 0 ) then ! 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 & & ) !積算値のクリア call StoreClean call Store2Clean end if end do LFrog !---------------------------------------------------------------- ! 出力ファイルのクローズ !---------------------------------------------------------------- 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), & ! & xz_Stab(DimXMin:DimXMax, DimZMin:DimZMax), & & xz_StabTemp(DimXMin:DimXMax, DimZMin:DimZMax), & & xz_StabMolWt(DimXMin:DimXMax, DimZMin:DimZMax), & ! & 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