!= Program Arare ! ! Authors:: SUGIYAMA Ko-ichiro, ODAKA Masatsugu ! Version:: $Id: arare_081006.f90,v 1.1 2009-03-05 05:39:41 yamasita Exp $ ! Tag Name:: $Name: arare4-20100305 $ ! 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 basicset, only : basicset_init, & & xz_DensBasicZ, xza_MixRtBasicZ, xz_PotTempBasicZ, & & xz_VelSoundBasicZ, xz_ExnerBasicZ ! 積算値管理モジュール use StorePotTemp, only : StorePotTemp_init, StorePotTempClean, & & StorePotTempCond 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 DynImpFuncMarscond, 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_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 ! (2008/05/16 山下達也追加) 凝結過程モジュール ! 凝結量計算, 凝結物質濃度計算, 凝結による潜熱の計算 ! use LatentHeatPerMass, only : xz_LatHeatPerMassNl ! use MassCondense, only : xz_LatHeatPerMassNl, & ! & xz_SatRatioNs, & ! & xz_PotTempNs, & ! & xz_ExnerNs, & ! & xz_DensCloudNs, & ! & xz_MassCondNs ! use DensityCloud, only : xz_DensCloudNs, & ! & DelTimeShort, & ! & pz_VelXNl, & ! & xr_VelZNl, & ! & xz_DensClousNl, & ! & xz_MassCondNs, & ! & xz_DensCloudAs ! use LatentHeat, only : xz_MassCondNs, & ! & xz_LatHeatPerMassNl, & ! & ss_QCond ! 湿潤気塊の浮力計算用モジュール 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(:,:,:) real(8), allocatable :: xz_TendPotTempNl(:,:) ! 長い時間ステップで評価すべき温位の項を格納する配列 ! 2008/06/20 山下達也 追加 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_MassCondNs(:,:) real(8), allocatable :: xz_MassCondNl(:,:) ! 凝結熱による温位変化率 real(8), allocatable :: xz_QCond(:,:) ! 作業変数 real(8), allocatable :: worknum(:,:) !コマンドライン引数の解釈 ! 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 StoreMixRt_init( ) call StoreBuoy_init( ) call StoreStab_init( ) write(*,*) "OK" !内部変数の初期化. とりあえずゼロを入れて値を確定させておく. call ArareAlloc write(*,*) "OK" !予報変数の初期化 ! ReStartFile が設定されている場合にはファイルを読み込み, ! 設定されていない場合には, デフォルトの基本場と擾乱場を作る. write(*,*) InitFile if (trim(InitFile) /= '') then !基本場, 擾乱場の初期値を netCDF ファイルから取得する. call ReStartFile_Get( & & ReStartTime, & & xz_PotTempBl, xz_ExnerBl, pz_VelXBl, xr_VelZBl, & & xza_MixRtBl, xz_KmBl, xz_KhBl, & & xz_DensCloudBl, xz_SatRatioBl, & & xz_PotTempNl, xz_ExnerNl, pz_VelXNl, xr_VelZNl, & & xza_MixRtNl, xz_KmNl, xz_KhNl, & & xz_DensCloudNl, xz_SatRatioNl ) write(*,*) "OK" else !デフォルト設定の基本場, 擾乱場を作成する. call BasicEnv() call DisturbEnv(cfgfile, & & xz_PotTempBl, xz_ExnerBl, pz_VelXBl, xr_VelZBl, & & xza_MixRtBl, xz_KmBl, xz_KhBl, & & xz_DensCloudBl, xz_SatRatioBl ) ! 1 ループ目だけは, bl と nl の値を同じにしておく. ! 1 ステップ目はオイラー法で解く必要があるが, ! 1 ステップ目とそれ以外のステップを別々にコーディングしたくない為 xz_PotTempNl = xz_PotTempBl xz_ExnerNl = xz_ExnerBl pz_VelXNl = pz_VelXBl xr_VelZNl = xr_VelZBl xza_MixRtNl = xza_MixRtBl xz_KmNl = xz_KmBl xz_KhNl = xz_KhBl xz_DensCloudNl = xz_DensCloudBl xz_SatRatioNl = xz_SatRatioBl xz_PotTempSum = xz_PotTempBl + xz_PotTempBasicZ xz_ExnerSum = xz_ExnerBl + xz_ExnerBasicZ xz_TempSum = xz_PotTempSum * xz_ExnerSum 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 ) 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 ) 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_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) write(*,*) "OK" 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) ! write(*,*) "OK" !------------------------------------------------------------- ! 短い時間ステップの初期値作成 !------------------------------------------------------------- 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 LatentHeatPerMass( xz_LatHeatPerMassNl ) !(out) 単位質量あたりの潜熱 ! write(*,*) "OK" !=== 凝結量の計算 call MassCondense( & & xz_LatHeatPerMassNl, & !(in) 単位質量あたりの潜熱 & xz_SatRatioNs, & !(in) 過飽和度 & xz_PotTempNs, & !(in) 温位 & xz_ExnerNs, & !(in) 無次元圧力 & xz_DensCloudNs, & !(in) 雲密度 & xz_MassCondNs) !(out)凝結量 ! write(*,*) "OK" !=== 凝結物質濃度の計算 call DensityCloud( & & xz_DensCloudNs, & !(in) 現在の時間の雲密度 & DelTimeShort, & !(in) 時間間隔 & pz_VelXNl, xr_VelZNl, xz_DensCloudNl, & !(in) フラックス項の計算で用いる値 & xz_MassCondNs, & !(in) 凝結量 & xz_DensCloudAs) !(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_QCond) !(out)凝結熱による温位変化率 ! write(*,*) "OKOK" ! ! (2008/06/16 山下達也) ! !=== 過飽和度の計算 ! call SaturationRatio( & ! & xz_ExnerAs, & !(in) ! & xz_PotTempAs, & !(in) ! & xz_SatRatioAs) !(out) ! ! do k = DimZMin, DimZMax ! write(*,*) "satratio",xz_SatRatioAs(1,k) ! end do ! ! 温位の移流計算. !(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_ExnerBl ) & ! & + 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 ) ! write(*,*) "OK" !------------------------------------------------------------- ! 速度 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) !境界条件 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) ! do k = DimZMin, DimZMax ! write(*,*) "satratio",xz_SatRatioAs(1,k) ! end do !------------------------------------------------------------- ! 短い時間ステップのループを回すための処置 !------------------------------------------------------------- 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 !---------------------------------------------------------------- ! 最終的な短い時間ステップでの値を長い時間ステップでの値とみなす !---------------------------------------------------------------- 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_xza(xza_MixRtAl, xza_MixRtNl, xza_MixRtBl) ! call AsselinFilter_xz(xz_DensCloudAl, xz_DensCloudNl, xz_DensCloudBl) 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(山下 達也) ! 基本場と擾乱場の和を出力する為に追加 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_PotTempSum, & & xz_TempSum, & & xz_ExnerNl, & & pz_VelXNl, & & xr_VelZNl, & & xza_MixRtNl, & & xz_KmNl, & & xz_KhNl, & & xz_DensCloudNl, & & xz_SatRatioNl & & ) !積算値のクリア call StorePotTempClean call StoreMixRtClean call StoreBuoyClean call StoreStabClean 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_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), & & worknum(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_SatRatioBl = 0.0d0; xz_SatRatioNl = 0.0d0 xz_SatRatioAl = 0.0d0 xz_SatRatioNs = 0.0d0; xz_SatRatioAs = 0.0d0 xz_QCond = 0.0d0 DelTimeEular = 0.0d0 DelTimeLFrog = 0.0d0; NStepEular = 0.0d0 NstepLFrog = 0.0d0 worknum = 0.0d0 end subroutine ArareAlloc end program arare