!--------------------------------------------------------------------- ! Copyright (C) GFD Dennou Club, 2004, 2005. All rights reserved. !--------------------------------------------------------------------- != Program Arare ! ! * Developer: SUGIYAMA Ko-ichiro ! * Version: $ ! * Tag Name: $Name: arare3j $ ! * Change History: ! !== Overview ! !非静力学モデル deepconv/arare. ! !== Error Handling ! !== Known Bugs ! !== Note ! ! * 方程式系は準圧縮系. ! !== Future Plans ! program arare ! !非静力学モデル deepconv/arare. ! !モジュール読み込み use chemcalc use fileset use debugset use timeset use gridset use basicset use damping use timefilter use boundary use cflcheck ! use RestartFileIO use HistoryFileIO use DynFunc use DynImpFunc use NumDiffusion use Turbulence use MoistAdjust use WarmRainPrm use MoistBuoy !暗黙の型宣言禁止 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_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 :: xz_MixRtBl(:,:,:) real(8), allocatable :: xz_MixRtNl(:,:,:) real(8), allocatable :: xz_MixRtAl(:,:,:) real(8), allocatable :: pz_AccelVelXNl(:,:) real(8), allocatable :: xr_AccelVelZNl(:,:) real(8) :: Time real(8) :: ReStartTime(2) real(8), allocatable :: DelTimeLFrog(:) real(8) :: DelTimeEular integer, allocatable :: NStepEular(:) integer :: NStepLFrog integer :: t, tau integer :: spc !---------------------------------------------------------------------- ! 変数参照型モジュールの初期化 !---------------------------------------------------------------------- !NAMELIST ファイルの取得 call getarg( 1, cfgfile ) write(*,*) "Input NAMELIST file: ", cfgfile !I/O ファイル名の初期化 ! NAMELIST ファイル名を指定し, deepconv/arare の ! 出力ファイル名を NAMELIST から得る call fileset_init(cfgfile) !デバグ設定の初期化 ! NAMELIST から情報を得て, デバッグ出力スイッチの切替えを行う. call debugset_init(cfgfile) !時刻に関する設定の初期化 ! NAMELIST から必要な情報を読み取り, 時間関連の変数の設定を行う. call timeset_init(cfgfile) !格子点情報の初期化 ! NAMELIST から情報を得て, 格子点を計算する call gridset_init(cfgfile) !基本場の情報の初期化 ! NAMELIST から情報を得て, 基本場を設定する. call basicset_init(cfgfile) !内部変数の初期化. とりあえずゼロを入れて値を確定させておく. call ArareAlloc !予報変数の初期化 ! ReStartFile が設定されている場合にはファイルを読み込み, ! 設定されていない場合には, デフォルトの基本場と擾乱場を作る. ! if (InitFile /= '') then ! !リスタートファイルの読み込み ! call basicset_file_init(InitFile) ! ! !擾乱場の初期値を netCDF ファイルから取得する. ! call ReStartFile_Get( & ! & ReStartTime, & ! & pz_VelXBl, xr_VelZBl, xz_ExnerBl, & ! & xz_PotTempBl, xz_KmBl, xz_MixRtBl, & ! & pz_VelXNl, xr_VelZNl, xz_ExnerNl, & ! & xz_PotTempNl, xz_KmNl, xz_MixRtNl ) ! else !デフォルト設定の基本場, 擾乱場を作成する. call BasicEnv() call DisturbEnv(cfgfile, & & pz_VelXBl, xr_VelZBl, & & xz_ExnerBl, xz_PotTempBl, & & xz_MixRtBl ) ! end if !---------------------------------------------------------------------- ! パッケージ型モジュールの初期化 ! デフォルトの値から変更する必要のあるルーチンのみ初期化 !---------------------------------------------------------------------- call Damping_Init( cfgfile ) !波の減衰係数の初期化 call DynImpFunc_Init() !陰解法の初期化 call NumDiffusion_Init() !数値拡散項の初期化 call Turbulence_Init() !乱流計算の初期化 call MoistAdjust_Init() !湿潤飽和調整法の初期化 call WarmRainPrm_Init() !暖かい雨のパラメタリゼーションの初期化 call MoistBuoy_Init() !分子量に対する浮力計算ルーチンの初期化 !---------------------------------------------------------------------- ! 時刻とループ回数の初期化 !---------------------------------------------------------------------- ! if (InitFile /= '') then ! Time = ReStartTime(2) !リスタート開始時刻 ! NstepLFrog = NstepLong - 1 ! NstepEular = NstepShort ! DelTimeLFrog = DelTimeLong * 2.0d0 ! DelTimeEular = DelTimeShort ! else Time = 0.0d0 !計算開始時刻 NstepLFrog = NstepLong NstepEular = NstepShort NstepEular(1) = NstepShort * 5.0d-1 ! 1 ループ目だけ DelTimeLFrog = DelTimeLong * 2.0d0 DelTimeLFrog(1) = DelTimeLong ! 1 ループ目だけ DelTimeEular = DelTimeShort ! 1 ループ目だけは, bl と nl の値を同じにしておく. ! 1 ステップ目はオイラー法で解く必要があるが, ! 1 ステップ目とそれ以外のステップを別々に書きたくないので. pz_VelXNl = pz_VelXBl xr_VelZNl = xr_VelZBl xz_ExnerNl = xz_ExnerBl xz_PotTempNl = xz_PotTempBl xz_KmNl = xz_KmBl xz_KhNl = xz_KhBl xz_MixRtNl = xz_MixRtBl ! end if !---------------------------------------------------------------- ! ファイル入出力 !---------------------------------------------------------------- !ファイルオープン call HistoryFile_Open( ) ! if (InitFile /= '') then ! call HistoryFile_Output( ReStartTime(1), & ! & xz_CNcrNl, xz_CLcrNl, xz_EVrvNl, xz_PRrNl,& ! & pz_VelXBl, xr_VelZBl, xz_ExnerBl, & ! & xz_PotTempBl, xz_KmBl, xz_MixRtBl ) ! ! call HistoryFile_Output( ReStartTime(2), & ! & xz_CNcrNl, xz_CLcrNl, xz_EVrvNl, xz_PRrNl,& ! & pz_VelXNl, xr_VelZNl, xz_ExnerNl, & ! & xz_PotTempNl, xz_KmNl, xz_MixRtNl ) ! end if !時刻 0 の場合には出力 if ( Time == 0 ) then call HistoryFile_Output( & & Time, & & xz_PotTempNl, & & xz_ExnerNl, & & pz_VelXNl, & & xr_VelZNl, & & xz_MixRtNl, & & xz_KmNl ) call HistoryFile_OutputAnal( & & xz_ExnerNl, & & xz_PotTempNl, & & xz_MixRtNl ) end if !---------------------------------------------------------------------- ! 設定のチェック !---------------------------------------------------------------------- !CFL 条件のチェック call CFLCheckTimeShort( xz_VelSoundBasicZ ) !---------------------------------------------------------------------- ! 数値積分 !---------------------------------------------------------------------- LFrog: do t = 1, NstepLFrog !時刻の設定 Time = Time + DelTimeLong write(*,*) Time !---------------------------------------------------------------- ! 渦粘性係数, 渦拡散係数を求める. !---------------------------------------------------------------- xz_KmAl = & & xz_KmBl & & + DelTimeLFrog(t) & & * ( & & + xz_AdvScalar(xz_KmNl, pz_VelXNl, xr_VelZNl) & & + xz_BuoyKm(xz_PotTempBl) & & + xz_MoistBuoyKm(xz_ExnerBl, xz_PotTempBl, xz_MixRtBl) & & + xz_ShearKm(xz_KmBl, pz_VelXBl, xr_VelZBl) & & + xz_NumDiffScalar(xz_KmBl) & & + xz_DispKm(xz_KmBl) & & ) !負の処理 where (xz_KmAl < 0.0d0) xz_KmAl = 0.0d0 end where !値の上限 (中島健介, 1994, 学位論文参照) where (xz_KmAl > 800.0d0) xz_KmAl = 800.0d0 end where !境界条件 xz_KmAl = xz_Boundary_Cyc_Sym(xz_KmAl) !アッセリンの時間フィルタ. Nl の値をフィルタリング xz_KmNl = AsselinFilter(xz_KmAl, xz_KmNl, xz_KmBl) !渦拡散定数を決める xz_KhAl = 3.0d0 * xz_KmAl !---------------------------------------------------------------- ! 温位の移流計算. !---------------------------------------------------------------- xz_PotTempAl = & & xz_PotTempBl & & + DelTimeLFrog(t) & & * ( & & + xz_AdvScalar( xz_PotTempNl, pz_VelXNl, xr_VelZNl) & & + xz_TurbScalar(xz_PotTempBl, xz_KhBl) & & + xz_NumDiffScalar(xz_PotTempBl) & & + xz_AdvScalar( xz_PotTempBasicZ, pz_VelXNl, xr_VelZNl)& & + xz_TurbScalar(xz_PotTempBasicZ, xz_KhBl) & & + xz_DispQ( xz_KmBl ) & ! & + xz_RadiQ( xz_ExnerBl ) & ! & + xz_HFluxSfc(pz_VelXBl, xz_PotTempBl, xz_ExnerBl) & & ) !境界条件 xz_PotTempAl = xz_Boundary_Cyc_Sym( xz_PotTempAl ) !アッセリンの時間フィルタ. Nl の値をフィルタリング xz_PotTempNl = AsselinFilter(xz_PotTempAl, xz_PotTempNl, xz_PotTempBl) !---------------------------------------------------------------- ! 凝縮成分の混合比の移流計算. !---------------------------------------------------------------- do spc = 1, SpcNum !時間積分 xz_MixRtAl(:,:,spc) = & & xz_MixRtBl(:,:,spc) & & + DelTimeLFrog(t) & & * ( & & + xz_AdvScalar(xz_MixRtNl(:,:,spc), pz_VelXNl, xr_VelZNl) & & + xz_TurbScalar(xz_MixRtBl(:,:,spc), xz_KhBl) & & + xz_NumDiffScalar(xz_MixRtBl(:,:,spc)) & & + xz_AdvScalar(xz_MixRtBasicZ(:,:,spc), xr_VelZNl, xr_VelZNl)& & + xz_TurbScalar(xz_MixRtBasicZ(:,:,spc), xz_KhBl) & ! & + xz_FallRain(xz_MixRtBl(:,:,spc), spc) & & ) !境界条件 xz_MixRtAl(:,:,spc) = xz_Boundary_Cyc_Sym(xz_MixRtAl(:,:,spc)) !アッセリンの時間フィルタ. Nl の値をフィルタリング xz_MixRtNl(:,:,spc) & & = AsselinFilter(xz_MixRtAl(:,:,spc), xz_MixRtNl(:,:,spc), xz_MixRtBl(:,:,spc)) !移流によって負になった部分を埋める call FillNegative( xz_DensBasicZ, xz_MixRtAl(:,:,spc) ) xz_MixRtAl(:,:,spc) = xz_Boundary_Cyc_Sym(xz_MixRtAl(:,:,spc)) end do !------------------------------------------------------------- ! 暖かい雨のパラメタリゼーション. ! 雲<-->雨 の変換を行う !------------------------------------------------------------- ! xz_MixRtAl(:,:,:) = & ! & xz_MixRtAl(:,:,:) & ! & + DelTimeLFrog(t) & ! & * ( & ! & + xz_Cloud2Rain( xz_MixRtAl(:,:,:) ) & !! & + xz_Cloud2RainNH4SH( xz_MixRtAl(:,:,:) ) & ! & ) ! !境界条件 ! do spc = 1, SpcNum ! xz_MixRtAl(:,:,spc) = xz_Boundary_Cyc_Sym(xz_MixRtAl(:,:,spc)) ! end do !------------------------------------------------------------- ! 湿潤飽和調節法 ! 水蒸気<-->雲の変換を行う. !------------------------------------------------------------- call MoistAdjustSvapPress( xz_ExnerNl, xz_PotTempAl, xz_MixRtAl(:,:,:) ) ! call MoistAdjustSvapPress2( xz_ExnerNl, xz_PotTempAl, xz_MixRtAl(:,:,:) ) ! call MoistAdjustNH4SH2( xz_ExnerNl, xz_PotTempAl, xz_MixRtAl(:,:,:) ) !------------------------------------------------------------- ! 暖かい雨のパラメタリゼーション. ! 蒸気<-->雨 の変換を行う !------------------------------------------------------------- ! xz_PotTempAl = & ! & xz_PotTempAl & ! & + DelTimeLFrog(t) & ! & * ( & ! & + xz_LatentHeatRain2Gas( & ! & xz_ExnerNl, xz_PotTempAl, xz_MixRtAl & ! & ) & !! & + xz_LatentHeatRain2Gas( & !! & xz_ExnerNl, xz_PotTempAl, xz_MixRtAl & !! & ) & ! & ) / ( xz_ExnerBasicZ * CpDry ) ! !境界条件 ! xz_PotTempAl = xz_Boundary_Cyc_Sym( xz_PotTempAl ) ! xz_MixRtAl = & ! & xz_MixRtAl & ! & + DelTimeLFrog(t) & ! & * ( & ! & + xz_Rain2Gas(xz_ExnerNl, xz_PotTempAl, xz_MixRtAl) & !! & + xz_Rain2GasNH4SH(xz_ExnerNl, xz_PotTempAl, xz_MixRtAl) & ! & ) ! !境界条件 ! do spc = 1, SpcNum ! xz_MixRtAl(:,:,spc) = xz_Boundary_Cyc_Sym(xz_MixRtAl(:,:,spc)) ! end do !------------------------------------------------------------- ! 長い時間ステップでの, 速度の移流, 拡散, 数値粘性, 浮力 !------------------------------------------------------------- pz_AccelVelXNl = & & + pz_AdvVelX(pz_VelXNl, 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_VelZNl) & & + xr_Buoy(xz_PotTempNl) & & + xr_MoistBuoy(xz_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_MFluxSfc(pz_VelXNs) & & ) pz_VelXAs = pz_Boundary_Cyc_Sym(pz_VelXAs) !------------------------------------------------------------- ! エクスナー関数の計算 !------------------------------------------------------------- xz_ExnerAs = xz_Exner( & & xr_AccelVelZNl, & & pz_VelXNs, & & pz_VelXAs, & & xr_VelZNs, & & xz_ExnerNs) xz_ExnerAs = xz_Boundary_Cyc_Sym(xz_ExnerAs) !------------------------------------------------------------- ! 速度 w の計算 !------------------------------------------------------------- xr_VelZAs = & & xr_VelZNs & & + DelTimeEular & & * ( & & - xr_GradPi(xz_ExnerAs,xz_ExnerNs,pz_VelXNs,xr_VelZNs) & & + xr_AccelVelZNl & & ) xr_VelZAs = xr_Boundary_Cyc_AntiSym(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 !------------------------------------------------------------- ! スポンジ層 !------------------------------------------------------------- xz_ExnerAs = xz_ExnerAs - DelTimeLFrog(t) * xz_DampSponge_xz(xz_ExnerBl) pz_VelXAs = pz_VelXAs - DelTimeLFrog(t) * pz_DampSponge_pz(pz_VelXBl ) xr_VelZAs = xr_VelZAs - DelTimeLFrog(t) * xr_DampSponge_xr(xr_VelZBl ) ! xz_PotTempAl = xz_PotTempAl - DelTimeLFrog(t) * xz_DampSponge_xz(xz_PotTempBl) !---------------------------------------------------------------- ! 長い時間ステップのループを回すための処置 !---------------------------------------------------------------- xz_PotTempBl = xz_PotTempNl xz_MixRtBl = xz_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 xz_MixRtNl = xz_MixRtAl xz_ExnerNl = xz_ExnerAl pz_VelXNl = pz_VelXAl xr_VelZNl = xr_VelZAl xz_KmNl = xz_KmAl xz_KhNl = xz_KhAl !---------------------------------------------------------------- ! ファイル出力 !---------------------------------------------------------------- if ( mod( Time, TimeDisp ) == 0 ) then ! CFL 条件のチェック call CFLCheckTimeLongVelX( pz_VelXNl ) call CFLCheckTimeLongVelZ( xr_VelZNl ) !ヒストリファイル出力 call HistoryFile_Output( & & Time, & & xz_PotTempNl, & & xz_ExnerNl, & & pz_VelXNl, & & xr_VelZNl, & & xz_MixRtNl, & & xz_KmNl ) call HistoryFile_OutputAnal( & & xz_ExnerNl, & & xz_PotTempNl, & & xz_MixRtNl ) end if end do LFrog !---------------------------------------------------------------- ! 出力ファイルのクローズ !---------------------------------------------------------------- call HistoryFile_Close !---------------------------------------------------------------- ! リスタートファイルの作成 !---------------------------------------------------------------- ! call ReStartFile_Open( ) ! call ReStartFile_OutPut( & ! & Time, & ! & pz_VelXBl, xr_VelZBl, xz_ExnerBl, & ! & xz_PotTempBl, xz_KmBl, xz_MixRtBl ) ! call ReStartFile_OutPut( & ! & Time, & ! & pz_VelXNl, xr_VelZNl, xz_ExnerNl, & ! & xz_PotTempNl, xz_KmNl, xz_MixRtNl ) ! 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_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), & ! & xz_MixRtBl(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum), & & xz_MixRtNl(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum), & & xz_MixRtAl(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum), & ! & pz_AccelVelXNl(DimXMin:DimXMax, DimZMin:DimZMax), & & xr_AccelVelZNl(DimXMin:DimXMax, DimZMin:DimZMax), & ! & 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 xz_MixRtBl = 0.0d0; xz_MixRtNl = 0.0d0; xz_MixRtAl = 0.0d0 pz_AccelVelXNl = 0.0d0 xr_AccelVelZNl = 0.0d0 DelTimeEular = 0.0d0 DelTimeLFrog = 0.0d0; NStepEular = 0.0d0 NstepLFrog = 0.0d0 end subroutine ArareAlloc end program arare