!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004. All rights reserved.
!---------------------------------------------------------------------
                                                                 !=begin
!= Program Arare
!
!   * Developer: SUGIYAMA Ko-ichiro, KITAMORI Taichi
!   * Version: $Id: arare.f90,v 1.1.1.1.2.20 2005/11/30 11:59:00 kitamo Exp $ 
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview 
!
!非静力学モデル deepconv/arare. 
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!  * 方程式系は準圧縮系.
!
!== Future Plans
!
                                                                 !=end
!-- (2005/04/21 小高正嗣) --
! 全体的なコメント
! - プログラム全体の流れはこれでよいと思われる.
! - 命名法(手続き名, 変数, 仮引数)は検討の余地あり
! - コメント文をどの程度書き込むか?
!   arare.f90 からコード解説の「メインプログラムの流れ」が自動生成
!   されるようになっているとよいだろう. 
!  

program arare
                                                                 !=begin

  !== Dependency
    !-- (2005/04/21 小高正嗣) --
    ! use される変数にもコメントを付けたほうがよいと思うが, うるさい? 
    ! コメントは各サブルーチン/モジュール内の記述と一致させる.
    !
  use dc_trace, only: BeginSub, &      ! デバッグ出力サブルーチン  
       &              EndSub,   &      ! デバッグ出力サブルーチン  
       &              DbgMessage       ! デバッグ出力サブルーチン  
  use dc_message                       ! only 属性不要 ?
  use fileset,  only: fileset_init,  & ! I/O ファイル名初期化サブルーチン  
       &              cfgfile,       & ! 設定ファイル
       &              ReStartFile,   & ! 入力用リスタートファイル
       &              NewReStartFile   ! 出力用リスタートファイル
  use debugset, only: debugset_init, & ! デバッグ設定初期化サブルーチン
       &              DebugOn          ! デバッグ出力スイッチ
  use gridset,  only: gridset_init,  & ! 格子点情報初期化サブルーチン
       &              DimXMin,       & ! x 方向の配列添字の下限 
       &              DimXMax,       & ! x 方向の配列添字の上限
       &              DimZMin,       & ! z 方向の配列添字の下限
       &              DimZMax,       & ! z 方向の配列添字の上限
       &              RegXMin,       & ! x 方向の物理領域を表す配列添字の下限
       &              RegXMax,       & ! x 方向の物理領域を表す配列添字の上限 
       &              RegZMin,       & ! z 方向の物理領域を表す配列添字の下限
       &              RegZMax,       & ! z 方向の物理領域を表す配列添字の上限 
       &              DelX,          & ! x 方向の格子間隔
       &              DelZ,          & ! z 方向の格子間隔
       &              s_Z              ! z 座標
  use timeset,  only: timeset_init,  & ! 時間積分設定初期化サブルーチン
       &              NstepShort,    & ! 長いタイムステップのステップ数
       &              NstepLong,     & ! 短いタイムステップのステップ数
       &              DelTimeShort,  & ! 長いタイムステップ
       &              DelTimeLong,   & ! 短いタイムステップ
       &              TimeDisp         ! ファイル出力する時間間隔
  use nameset,  only: nameset_init     ! 実験名初期化サブルーチン
  use physset,  only: physset_init,  & ! 物理定数の初期化サブルーチン
       &              PressSfc,      & ! 地表面圧力
       &              GasR,          & ! 気体定数
       &              Grav             ! 重力
  use linlib,   only: linlib_init      ! 行列計算ライブラリ初期化サブルーチン
  use bcset,    only: bcset_init       ! 境界条件設定サブルーチン
  use basicset, only: basicset_file_init, &  ! 基本場設定サブルーチン
       &              ss_DensBasicZ,      &  ! 基本場の密度
       &              ss_PotTempBasicZ,   &  ! 基本場の密度
       &              ss_ExnerBasicZ,     &  ! 基本場の密度
       &              ss_CpBasicZ,        &  ! 定圧比熱 
       &              ss_VelSoundBasicZ      ! 音速
  use arareset, only: arareset_init    ! 計算設定を読み込むサブルーチン
  use average,  only: ss_avr_fs,     & 
                        ! x 方向フラックス格子点からスカラー格子点への平均操作
       &              ss_avr_sf,     & 
                        ! z 方向フラックス格子点からスカラー格子点への平均操作
       &              sf_avr_ss,     &
                        ! スカラー格子点から z 方向フラックス格子点への平均操作
       &              fs_avr_ss
                        ! スカラー格子点から z 方向フラックス格子点への平均操作
  use summation,only: sumxz, sumx            ! 領域全体の和を計算
  use differentiate_center4,only: ss_dx_fs, ss_dz_sf, sf_dz_ss  ! 微分演算
  use gt4_history
                                                                 !=end
  !== 暗黙の型宣言禁止
  implicit none

  !== 変数配列
  real(8), allocatable  :: fs_VelX_as(:,:)    ! 水平速度      [t=(n-1)Δτ]
  real(8), allocatable  :: fs_VelX_ns(:,:)    ! 水平速度      [t=nΔτ]
  real(8), allocatable  :: sf_VelZ_as(:,:)    ! 鉛直速度      [t=(n-1)Δτ]
  real(8), allocatable  :: sf_VelZ_ns(:,:)    ! 鉛直速度      [t=nΔτ]
  real(8), allocatable  :: ss_Exner_ns(:,:)   ! エクスナー関数[t=nΔτ]
  real(8), allocatable  :: ss_Exner_as(:,:)   ! エクスナー関数[t=(n-1)Δτ]

  real(8), allocatable  :: ss_PotTemp_bl(:,:) ! 温位 [t=(n-1)Δt]
  real(8), allocatable  :: ss_PotTemp_nl(:,:) ! 温位 [t=nΔt    ]
  real(8), allocatable  :: ss_PotTemp_al(:,:) ! 温位 [t=(n+1)Δt]
  real(8), allocatable  :: ss_Exner_bl(:,:)   ! エクスナー関数 [t=(n-1)Δt]
  real(8), allocatable  :: ss_Exner_nl(:,:)   ! エクスナー関数 [t=nΔt    ]
  real(8), allocatable  :: ss_Exner_al(:,:)   ! エクスナー関数 [t=(n+1)Δt]
  real(8), allocatable  :: fs_VelX_bl(:,:)    ! 水平速度 [t=(n-1)Δt]
  real(8), allocatable  :: fs_VelX_nl(:,:)    ! 水平速度 [t=nΔt    ]
  real(8), allocatable  :: fs_VelX_al(:,:)    ! 水平速度 [t=(n+1)Δt]
  real(8), allocatable  :: sf_VelZ_bl(:,:)    ! 鉛直速度 [t=(n-1)Δt]
  real(8), allocatable  :: sf_VelZ_nl(:,:)    ! 鉛直速度 [t=nΔt    ]
  real(8), allocatable  :: sf_VelZ_al(:,:)    ! 鉛直速度 [t=(n+1)Δt]

  real(8), allocatable  :: fs_AdvDiffX_nl(:,:)! 水平方向の運動方程式の移流項と拡散項   [t=nΔt]
  real(8), allocatable  :: sf_AdvDiffZ_nl(:,:)! 鉛直方向の運動方程式の移流項と拡散項   [t=nΔt]
  
  real(8), allocatable  :: ss_Km_bl(:,:)      ! 渦粘性係数 [t=(n-1)Δt]
  real(8), allocatable  :: ss_Km_nl(:,:)      ! 渦粘性係数 [t=nΔt    ]
  real(8), allocatable  :: ss_Km_al(:,:)      ! 渦粘性係数 [t=(n+1)Δt]
  real(8), allocatable  :: ss_Kh_bl(:,:)      ! 渦拡散係数 [t=(n-1)Δt]
  real(8), allocatable  :: ss_Kh_nl(:,:)      ! 渦拡散係数 [t=nΔt    ]
  real(8), allocatable  :: ss_Kh_al(:,:)      ! 渦拡散係数 [t=(n+1)Δt]

  real(8), allocatable  :: sf_SfcTemp_al(:)      ! 地表面温度 [t=(n+1)Δt    ]
  real(8), allocatable  :: sf_SfcHeatFlux(:)     ! 地表面熱フラックス
  real(8), allocatable  :: ff_SfcMomentumFlux(:) ! 地表面運動量フラックス
  real(8), allocatable  :: ss_QRad(:,:)          ! 放射加熱
  real(8), allocatable  :: ss_QDis(:,:)          ! 散逸加熱

  real(8), allocatable  :: ss_PotTempTendAdv_nl(:,:)  ! 移流による温位変化
  real(8), allocatable  :: ss_PotTempTendSrc_nl(:,:)  ! 生成による温位変化
  real(8), allocatable  :: ss_PotTempTendTurb_nl(:,:) ! 乱流拡散による温位変化
  real(8), allocatable  :: s_PotTempTendAdvZ(:)  ! 移流による温位変化(水平平均)
  real(8), allocatable  :: s_PotTempTendSrcZ(:)  ! 生成による温位変化(水平平均)
  real(8), allocatable  :: s_PotTempTendTurbZ(:) 
                                      ! 乱流拡散による温位変化(水平平均)
  real(8), allocatable  :: s_PotTempTendRadZ(:)
                                      ! 放射加熱による温位変化(水平平均)
  real(8), allocatable  :: s_PotTempTendDisZ(:)
                                      ! 散逸加熱による温位変化(水平平均)
  real(8), allocatable  :: s_PotTempTendSfcZ(:) 
                                      ! 地表面加熱による温位変化(水平平均)
  real(8), allocatable  :: sf_VelZTendPressGrad_ns(:,:)
                                      ! 圧力勾配による鉛直速度変化
  real(8), allocatable  :: fs_VelXTendAdv_nl(:,:)  ! 移流による水平速度変化
  real(8), allocatable  :: fs_VelXTendTurb_nl(:,:) ! 乱流拡散による水平速度変化
  real(8), allocatable  :: fs_VelXTendNumDiff_nl(:,:) ! 数値拡散による水平速度変化
  real(8), allocatable  :: sf_VelZTendAdv_nl(:,:)  ! 移流による鉛直速度変化
  real(8), allocatable  :: sf_VelZTendTurb_nl(:,:) ! 乱流拡散による鉛直速度変化
  real(8), allocatable  :: sf_VelZTendBuoy_nl(:,:) ! 浮力による鉛直速度変化
  real(8), allocatable  :: sf_VelZTendNumDiff_nl(:,:) ! 数値拡散による鉛直速度変化
  real(8), allocatable  :: f_VelZTendPressGradZ(:)
                                     ! 圧力勾配による鉛直速度変化(水平平均)
  real(8), allocatable  :: f_VelZTendAdvZ(:) 
                                     ! 移流による鉛直速度変化(水平平均)
  real(8), allocatable  :: f_VelZTendTurbZ(:)
                                     ! 乱流拡散による鉛直速度変化(水平平均)
  real(8), allocatable  :: f_VelZTendBuoyZ(:)
                                     ! 浮力による鉛直速度変化(水平平均)

  real(8), allocatable  :: ss_KineticEnergy_nl(:,:)   ! 運動エネルギー
  real(8), allocatable  :: ss_PotentialEnergy_nl(:,:) ! ポテンシャルエネルギー
  real(8), allocatable  :: ss_AvailPotentialEnergy_nl(:,:) ! 有効ポテンシャルエネルギー
  real(8), allocatable  :: ss_ElasticEnergy_nl(:,:) ! ポテンシャルエネルギー
  real(8)               :: MassTotal             ! 領域全体の質量
  real(8)               :: KineticEnergyTotal    ! 領域全体の運動エネルギー
  real(8)               :: PotentialEnergyTotal  ! 領域全体のポテンシャルエネルギー
  real(8)               :: AvailPotentialEnergyTotal  ! 領域全体の有効ポテンシャルエネルギー
  real(8)               :: ElasticEnergyTotal    ! 領域全体の弾性エネルギー

  real(8), allocatable  :: ss_EnergyTendRad_nl(:,:) ! 放射加熱によるエネルギー変化
  real(8), allocatable  :: ss_EnergyTendDis_nl(:,:) ! 放射加熱によるエネルギー変化
  real(8), allocatable  :: ss_EnergyTendVelTurb_nl(:,:)
                                    ! 運動量拡散によるエネルギー変化
  real(8), allocatable  :: ss_EnergyTendPotTempTurb_nl(:,:)
                                    ! 熱拡散によるエネルギー変化
  real(8), allocatable  :: ss_EnergyTendSfcHeatFlux_nl(:,:)
                                    ! 地表面熱フラックスによるエネルギー変化
  real(8), allocatable  :: ss_EnergyTendNonconsv_nl(:,:)
                                    ! 準圧縮系のエネルギー非保存項
  real(8), allocatable  :: ss_EnergyTendPotTempBasic_nl(:,:)
                                    ! 基本場温位によるエネルギー変化
  real(8)               :: EnergyTendRadTotal 
                                    ! 放射加熱による全領域のエネルギー変化
  real(8)               :: EnergyTendDisTotal
                                    ! 散逸加熱による全領域のエネルギー変化
  real(8)               :: EnergyTendVelTurbTotal
                                    ! 運動量拡散による全領域のエネルギー変化
  real(8)               :: EnergyTendPotTempTurbTotal
                                    ! 熱拡散による全領域のエネルギー変化
  real(8)               :: EnergyTendSfcHeatFluxTotal
                                    ! 地表面熱フラックスによるエネルギー変化
  real(8)               :: EnergyTendSfcMomentumFluxTotal
                                    ! 地表面運動量フラックスによるエネルギー変化
  real(8)               :: EnergyTendNonconsvTotal
                                    ! 全領域の準圧縮系のエネルギー非保存項
  real(8)               :: EnergyTendPotTempBasicTotal
                                    ! 全領域の基本場温位によるエネルギー変化
  real(8)               :: PotTempMean        ! 温位擾乱の領域全体平均値
                                              ! 診断量の計算のためだけに使用

  real(8)               :: Time               ! 計算継続時間
  real(8)               :: DelTime            ! 時間ステップ
  real(8)               :: ReStartTime(2)     ! リスタートファイルの時刻
  integer               :: t, tau             ! ループを回すための変数
  integer               :: NStep_leapfrog     ! ループを回す回数

  real(8), allocatable  :: ss_Z(:,:)          ! 2-D 座標(診断量計算で使う)

  integer               :: StartRunTime(8)    ! 計算開始時刻
  integer               :: EndRunTime(8)      ! 計算終了時刻
  integer               :: RunTime(8)         ! 実行時間
  integer               :: i,j,k                ! ループを回すための変数

  type(GT_HISTORY)   :: hst_pred, hst_diag    ! gt4_history のヒストリ名

  !== 初期化

  !=== 開始時刻の記録
  !fortran90 の組み込み関数 date_and_time を利用.
  call date_and_time(values = StartRunTime)

  !=== I/O ファイル名の初期化
  !NAMELIST ファイル名を指定し, deepconv/arare の
  !出力ファイル名を NAMELIST から得る
  call fileset_init

  !=== デバグ設定の初期化
  !NAMELIST から情報を得て, デバッグ出力スイッチの切替えを行う.  
  call debugset_init(cfgfile)

  !=== 格子点情報の初期化
  !NAMELIST から情報を得て, 格子点を計算する
  call gridset_init(cfgfile)
  
  !=== 時間積分の設定の初期化
  !NAMELIST から必要な情報を読み取り, 時間関連の変数の設定を行う. 
  call timeset_init(cfgfile)
  
  !=== 実験名の初期化
  !netCDF に書き込む実験名を NAMELIST から読み込む
  call nameset_init(cfgfile)
  
  !=== 物理定数の初期化
  !NAMELIST から物理定数を読み込む
  call physset_init(cfgfile)
  
  !=== 計算の設定を読み込む
  !NameList ファイルから arare 実行に必要なパラメタを取得する
  call arareset_init(cfgfile)

  !=== 境界条件を設定
  !NAMELIST から境界条件のタイプを設定する. 
  call bcset_init(cfgfile)

  !=== 出力ファイルの定義
  call OpenDataFile(hst_pred)
  call OpenDiagDataFile(hst_diag)

!== 変数の初期化

  !=== 配列の割り当て, 初期化
  call ArareAlloc

  !=== 2-D 座標配列の初期化
  do k = DimZmin, DimZMax
     ss_Z(:,k) = s_Z(k)
  end do

  !=== 基本場の値を設定
  !基本場の値を netCDF ファイルから取得する. 
  call basicset_file_init( ReStartFile )

  !=== 擾乱場の値を設定
  !擾乱場の初期値を netCDF ファイルから取得する. 
  call getdisturbvar( ReStartFile, ReStartTime,   &
       & fs_VelX_bl,     sf_VelZ_bl,  ss_Exner_bl,  &
       & ss_PotTemp_bl,  ss_Km_bl,    ss_Kh_bl,     &
       & fs_VelX_nl,     sf_VelZ_nl,  ss_Exner_nl,  &
       & ss_PotTemp_nl,  ss_Km_nl,    ss_Kh_nl      )

  !=== 行列計算ライブラリの初期化. 
  !利用する線形計算パッケージを NAMELIST から設定し, 
  !配列の大きさを決定する
  !  * 引数はスカラーの z 方向の計算領域の大きさ
  call linlib_init(cfgfile, RegZMax - RegZMin)
  
  !=== 時刻の初期化
  Time = ReStartTime(2)          !リスタート開始時刻
  DelTime = 2.0 * DelTimeLong    !leap-frog スキームでの積分時間

  !=== ループ回数の初期化
  !時刻 0 秒からの計算, すなわち ReStartTime(1) = 0 では, 
  !時刻 0 秒から積分時間 TimeInt まで積分する. 
  !このとき leap-frog スキームは NstepLong - 1 回ループさせることとなる. 
  !
  !継続計算, すなわち ReStartTime(1) /= 0 では, 
  !RestartTime(2) 秒から積分時間 TimeInt まで積分する. 
  !このとき leap-frog スキームは NstepLong1 回ループさせることとなる. 
  !継続計算の場合にはリスタートファイルの時刻から計算を行いたいからである.

  if ( ReStartTime(1) == 0.0 ) then 
     Nstep_leapfrog = NstepLong - 1
  else
     Nstep_leapfrog = NstepLong
  end if

  !=== 時刻 0 からの計算では, 必要に応じて初期値をファイル出力
  if ( ReStartTime(1) == 0.0 ) then 
     call OutputDataFile( hst_pred, ReStartTime(1), &
          & fs_VelX_bl, sf_VelZ_bl, ss_Exner_bl, &
          & ss_PotTemp_bl, ss_Km_bl, ss_Kh_bl )

     if ( mod( ReStartTime(2), TimeDisp ) == 0 ) then 
        call OutputDataFile( hst_pred,  ReStartTime(2), &
             & fs_VelX_nl, sf_VelZ_nl, ss_Exner_nl, &
             & ss_PotTemp_nl, ss_Km_nl, ss_Kh_nl )
     end if
  end if


!== CFL 条件のチェック

  call cflcheck

!== 長い時間ステップ. leap-frog スキームを利用
  do t = 1, Nstep_leapfrog
     Time = Time + DelTimeLong
     write(*,*) Time

     !=== 渦粘性係数, 渦拡散係数の計算
     call EddyViscosity(                                     &
          & ss_Km_bl,                                        & 
          & DelTime,                                         & !時間間隔
          & fs_VelX_nl, sf_VelZ_nl, ss_Km_nl,                & !移流
          & fs_VelX_bl, sf_VelZ_bl, ss_PotTemp_bl, ss_Km_bl, & !拡散等
          & ss_Km_al, ss_Kh_al )

     !=== 長い時間ステップでの温位の計算
     call Scalar(                                   & 
          & ss_PotTemp_bl,                          &
          & DelTime,                                & !時間間隔
          & fs_VelX_nl, sf_VelZ_nl, ss_PotTemp_nl,  & !移流
          & ss_PotTemp_bl, ss_Kh_bl,                & !乱流拡散, 等
          & ss_PotTemp_al,                          &  
          & ss_PotTempTendAdv_nl,                  & !移流
          & ss_PotTempTendSrc_nl,                  & !生成
          & ss_PotTempTendTurb_nl)                    !乱流拡散

     !=== 地表面温度の計算    
     call SurfaceTemp(sf_SfcTemp_al)

     !=== 地面フラックスの計算    
     call SurfaceHeatFlux(ss_CpBasicZ,     &
          &               ss_DensBasicZ,   &
          &               ss_PotTempBasicZ, &
          &               ss_ExnerBasicZ,   &
          &               fs_VelX_bl,       &
          &               ss_PotTemp_bl,    &
          &               ss_Exner_bl,      &
          &               sf_SfcTemp_al,    &
          &               sf_SfcHeatFlux )

     !=== 放射冷却率の計算
     call Radiation(ss_Exner_bl, ss_ExnerBasicZ, ss_QRad)

     !=== 散逸加熱の計算
     call DissipativeHeating(ss_Km_bl, ss_QDis)

     !=== フラックス, 冷却, 散逸加熱による温位変化を計算
     ss_PotTemp_al = ss_PotTemp_al + DelTime * ss_QRad

     ss_PotTemp_al = ss_PotTemp_al + DelTime * ss_QDis

     ss_PotTemp_al(:,RegZMin + 1) =   &
         &   ss_PotTemp_al(:,RegZMin + 1)        &
         & + DelTime * sf_SfcHeatFlux            & 
         &   / ( DelZ                            &
         &       * ss_DensBasicZ(:,RegZMin + 1)  & 
         &       * ss_CpBasicZ(:,RegZMin + 1)    &
         &      )

     !=== 温位変化の寄与の水平平均値を計算
     do i = RegXMin, RegXMax
        s_PotTempTendAdvZ(:) = s_PotTempTendAdvZ(:)  &
             & -  ss_PotTempTendAdv_nl(i,:) / (RegXMax - RegXMin + 1)

        s_PotTempTendSrcZ(:) = s_PotTempTendSrcZ(:)  &
             & -  ss_PotTempTendSrc_nl(i,:) / (RegXMax - RegXMin + 1)

        s_PotTempTendTurbZ(:) = s_PotTempTendTurbZ(:)  &
             & +  ss_PotTempTendTurb_nl(i,:) / (RegXMax - RegXMin + 1)

        s_PotTempTendRadZ(:) = s_PotTempTendRadZ(:)   &
             & + ss_QRad(i, :) / (RegXMax - RegXMin + 1)

        s_PotTempTendDisZ(:) = s_PotTempTendDisZ(:)   &
             & + ss_QDis(i, :) / (RegXMax - RegXMin + 1)

        s_PotTempTendSfcZ(RegZMin + 1) =   &
             & s_PotTempTendSfcZ(RegZMin + 1)   &
             & + sf_SfcHeatFlux(i)   &
             & / ( (RegXMax - RegXMin + 1)  &
             &     * DelZ   &
             &     * ss_DensBasicZ(i, RegZMin + 1)  &
             &     * ss_CpBasicZ(i, RegZMin + 1) )

     end do

     !=== 長い時間ステップでの速度 u の移流, 拡散, 粘性
     call AdvectDiffuseX(                     &
          & fs_VelX_nl, sf_VelZ_nl,           & !移流計算用入力値
          & fs_VelX_bl, sf_VelZ_bl, ss_Km_bl, & !乱流拡散, 数値粘性計算用入力値
          & fs_AdvDiffX_nl,  &                  !移流・拡散による速度変化率
          & fs_VelXTendAdv_nl,  &                  !移流による速度変化率
          & fs_VelXTendTurb_nl, &                  !乱流拡散による速度変化率
          & fs_VelXTendNumDiff_nl)                 !数値拡散による速度変化率
 
     !=== 長い時間ステップでの速度 w の移流, 拡散, 粘性
     call AdvectDiffuseZ(                          &
         & fs_VelX_nl, sf_VelZ_nl, ss_PotTemp_nl, & !移流, 浮力計算用入力値
         & fs_VelX_bl, sf_VelZ_bl, ss_Km_bl,      & !乱流拡散, 数値粘性計算用入力値
         & sf_AdvDiffZ_nl,   &
         & sf_VelZTendAdv_nl, sf_VelZTendBuoy_nl, & 
         & sf_VelZTendTurb_nl, sf_VelZTendNumDiff_nl)
     
    
     !=== 短いタイムステップの初期値を作る
     fs_VelX_ns  = fs_VelX_bl
     sf_VelZ_ns  = sf_VelZ_bl
     ss_Exner_ns = ss_Exner_bl 
    
     !=== 短いタイムステップでの積分
     do tau = 1, NstepShort 
        !==== 速度 u の計算
        call VelocityX(                                             &
             & fs_VelX_ns,                                          &
             & DelTimeShort,                                        &
             & fs_AdvDiffX_nl, fs_VelX_ns, sf_VelZ_ns, ss_Exner_ns, &
             & fs_VelX_as )
        
        !==== エクスナー関数の計算
        call Exner(                                    &
             & sf_AdvDiffZ_nl, fs_VelX_ns, fs_VelX_as, &
             &   sf_VelZ_ns, ss_Exner_ns,              &
             & ss_Exner_as )

        !==== 速度 w の計算
        call VelocityZ(                                &
             & sf_VelZ_ns,                             &
             & DelTimeShort,                           &
             & sf_AdvDiffZ_nl, fs_VelX_ns, sf_VelZ_ns, &
             &    ss_Exner_ns, ss_Exner_as,            &
             & sf_VelZ_as,  &
             & sf_VelZTendPressGrad_ns)

        !=== 地面フラックスの計算    
        call SurfaceMomentumFlux(ss_DensBasicZ,  &
                                 fs_VelX_ns,      &  ! ns ?
                                 ff_SfcMomentumFlux )

        !=== フラックスによる速度変化を計算
        fs_VelX_as(:, RegZMin + 1) = fs_VelX_as(:, RegZMin + 1)   &
             & + DelTimeShort * ff_SfcMomentumFlux                &
             &     / ( ss_DensBasicZ(:,RegZMin + 1) * DelZ )

        !=== 速度変化の寄与の水平平均値を計算
        do i = RegXMin, RegXMax
           f_VelZTendAdvZ(:) = f_VelZTendAdvZ(:)  &
                & -  sf_VelZTendAdv_nl(i,:) / (RegXMax - RegXMin + 1)

           f_VelZTendBuoyZ(:) = f_VelZTendBuoyZ(:)  &
                & -  sf_VelZTendBuoy_nl(i,:) / (RegXMax - RegXMin + 1)

           f_VelZTendTurbZ(:) = f_VelZTendTurbZ(:)  &
                & -  sf_VelZTendTurb_nl(i,:) / (RegXMax - RegXMin + 1)

           f_VelZTendPressGradZ(:) = f_VelZTendPressGradZ(:)  &
                & -  sf_VelZTendPressGrad_ns(i,:) / (RegXMax - RegXMin + 1)
        end do

        !==== 短い時間ステップのループを回すための処置
        ss_Exner_ns  = ss_Exner_as 
        fs_VelX_ns   = fs_VelX_as
        sf_VelZ_ns   = sf_VelZ_as
        
     end do
     
     !=== 最終的な短い時間ステップでの値を長い時間ステップでの値とみなす
     ss_Exner_al  = ss_Exner_as 
     fs_VelX_al   = fs_VelX_as
     sf_VelZ_al   = sf_VelZ_as     
     
     !=== アッセリンの時間フィルタ
     call AsselinTimeFilter(ss_PotTemp_al, ss_PotTemp_nl, ss_PotTemp_bl )
     call AsselinTimeFilter(ss_Exner_al,   ss_Exner_nl,   ss_Exner_bl )
     call AsselinTimeFilter(fs_VelX_al,    fs_VelX_nl,    fs_VelX_bl )
     call AsselinTimeFilter(sf_VelZ_al,    sf_VelZ_nl,    sf_VelZ_bl )

     !=== 長い時間ステップのループを回すための処置
     ss_PotTemp_bl = ss_PotTemp_nl
     ss_Exner_bl   = ss_Exner_nl 
     fs_VelX_bl    = fs_VelX_nl
     sf_VelZ_bl    = sf_VelZ_nl
     ss_Km_bl      = ss_Km_nl
     ss_Kh_bl      = ss_Kh_nl
     ss_PotTemp_nl = ss_PotTemp_al
     ss_Exner_nl   = ss_Exner_al 
     fs_VelX_nl    = fs_VelX_al
     sf_VelZ_nl    = sf_VelZ_al
     ss_Km_nl      = ss_Km_al
     ss_Kh_nl      = ss_Kh_al

     !=== ファイル出力
     if ( mod( Time, TimeDisp ) == 0 ) then 
        call DbgMessage(fmt="Time = %f", d=(/Time/))
        call OutputDataFile(hst_pred, Time, &
             &  fs_VelX_nl,    sf_VelZ_nl, ss_Exner_nl, &
             &  ss_PotTemp_nl, ss_Km_nl,   ss_Kh_nl )

        write(*,*) Time

        !=== 診断量の計算

        MassTotal =   &     ! 質量
             &   (RegXMax - RegXMin) * PressSfc / Grav  &
             &   * (    &
             &         ss_ExnerBasicZ(RegXMin,RegZMin)**(ss_CpBasicZ(RegXMin,RegZMin) / GasR)   &
             &       - ss_ExnerBasicZ(RegXMin,RegZMax)**(ss_CpBasicZ(RegXMin,RegZMax) / GasR)   &
             &   )  &
             & + PressSfc / GasR &
             &   * sumxz( &
             &              (ss_ExnerBasicZ + ss_Exner_nl)**((ss_CpBasicZ - GasR) / GasR)  &
             &            / (ss_PotTempBasicZ + ss_PotTemp_nl)  &
             &            - (ss_ExnerBasicZ)**((ss_CpBasicZ - GasR) / GasR)  &
             &            / ss_PotTempBasicZ  &
             &          )  
        
        PotTempMean = sumxz(ss_PotTemp_nl) / (RegXMax - RegXMin + 1) / (RegZMax - RegZMin + 1)
        
        ss_KineticEnergy_nl =   &   ! 運動エネルギー
             & ss_DensBasicZ  &
             & * (  &
             &       ss_avr_fs(fs_VelX_nl)**2.0d0  &
             &     + ss_avr_sf(sf_VelZ_nl)**2.0d0  &
             &   ) * DelX * DelZ
        KineticEnergyTotal = sumxz(ss_KineticEnergy_nl)
        
        ss_PotentialEnergy_nl =  &  ! 浮力によるポテンシャルエネルギー
             & - ss_DensBasicZ  &
             & * (    &
             &       Grav * ss_Z  &
             &       * PotTempMean &
             &       / ss_PotTempBasicZ  &
             &   ) * DelX * DelZ
        PotentialEnergyTotal = sumxz(ss_PotentialEnergy_nl)
   
        ss_AvailPotentialEnergy_nl =  &  ! 浮力によるポテンシャルエネルギー
             & - ss_DensBasicZ  &
             & * (    &
             &       Grav * ss_Z  &
             &       * (ss_PotTemp_nl - PotTempMean) &
             &       / ss_PotTempBasicZ  &
             &   ) * DelX * DelZ
        AvailPotentialEnergyTotal = sumxz(ss_AvailPotentialEnergy_nl)
        
        ss_ElasticEnergy_nl =  &   ! 弾性エネルギー
             & ss_DensBasicZ  &
             & * (   &
             &       ss_CpBasicZ * ss_PotTempBasicZ  &
             &     * ss_Exner_nl / ss_VelSoundBasicZ  &
             &   )**2.0d0 * 0.5d0 * DelX * DelZ
        ElasticEnergyTotal = sumxz(ss_ElasticEnergy_nl)        
        
        !=== エネルギー変化率の計算
   
        ss_EnergyTendRad_nl =   &      ! 放射冷却
             & - ss_DensBasicZ   &
             &   * Grav   &
             &   * ss_Z  & 
             &   * ss_Qrad &
             &   / ss_PotTempBasicZ  &
             &   * DelX * DelZ
        EnergyTendRadTotal = sumxz(ss_EnergyTendRad_nl)
        
        ss_EnergyTendDis_nl =   &      ! 散逸加熱
             & - ss_DensBasicZ   &
             &   * Grav   &
             &   * ss_Z  & 
             &   * ss_Qdis &
             &   / ss_PotTempBasicZ  &
             &   * DelX * DelZ
        EnergyTendDisTotal = sumxz(ss_EnergyTendDis_nl)
        
        ss_EnergyTendPotTempTurb_nl =   &   ! 乱流による熱拡散
             & - ss_DensBasicZ   &
             &   * Grav   &
             &   * ss_Z  & 
             &   * ss_PotTempTendTurb_nl &
             &   / ss_PotTempBasicZ &
             &   * DelX * DelZ
        EnergyTendPotTempTurbTotal = sumxz(ss_EnergyTendDis_nl)
        
        ss_EnergyTendVelTurb_nl =  &    ! 乱流による運動量拡散
             &   ss_DensBasicZ & 
             & * (  &
             &       ss_avr_fs(fs_VelX_nl * fs_VelXTendTurb_nl)  &
             &     + ss_avr_sf(sf_VelZ_nl * sf_VelZTendTurb_nl)  &
             &   )  &
             & * DelX * DelZ
        EnergyTendVelTurbTotal = sumxz(ss_EnergyTendVelTurb_nl)
   
        ss_EnergyTendNonconsv_nl =  &   ! 準圧縮方程式系のエネルギー非保存項
             & (    &
             &     ss_avr_fs(fs_VelX_nl)**2.0d0   &
             &   + ss_avr_sf(sf_VelZ_nl)**2.0d0   &
             &   - ss_PotTemp_nl  &
             &     * Grav  &
             &     * ss_Z  &
             &     / ss_PotTempBasicZ &
             & )    &
             & * (  &
             &       (ss_ExnerBasicZ + ss_Exner_as)**((ss_CpBasicZ - GasR) / GasR)  &
             &       / (ss_PotTempBasicZ + ss_PotTemp_nl)  &
             &     - (ss_ExnerBasicZ + ss_Exner_ns)**((ss_CpBasicZ - GasR) / GasR)  &
             &       / (ss_PotTempBasicZ + ss_PotTemp_nl)  &
             &   ) / DelTimeShort * DelX * DelZ
        EnergyTendNonconsvTotal =  &
             & sumxz(ss_EnergyTendNonconsv_nl)
   
        ss_EnergyTendSfcHeatFlux_nl(:,RegZMin + 1) = &   ! 地表面熱フラックス
             & - Grav   &
             &   * ss_Z(:,RegZMin + 1)  &
             &   * sf_SfcHeatFlux(:)    &
             &   * DelX     & 
             &   / (   ss_CpBasicZ(:,RegZMin + 1)  &
             &       * ss_PotTempBasicZ(:,RegZMin + 1) & 
             &     ) 
        EnergyTendSfcHeatFluxTotal = & 
             & sumx(ss_EnergyTendSfcHeatFlux_nl(:,RegZMin + 1))
        
        EnergyTendSfcMomentumFluxTotal =  &     ! 地表面運動量フラックス
             & sumx(ff_SfcMomentumFlux * DelX) 
        
        ss_EnergyTendPotTempBasic_nl =  &       ! 基本場温位からのエネルギー
             & ss_DensBasicZ * Grav * ss_Z  / ss_PotTempBasicZ  &
             & * ss_avr_sf(sf_VelZ_nl)  &
             & * ( 1 + ss_PotTemp_nl / ss_PotTempBasicZ )  &
             & * ss_avr_sf(sf_dz_ss(ss_PotTempBasicZ))
        EnergyTendPotTempBasicTotal = sumxz(ss_EnergyTendPotTempBasic_nl)
    
        s_PotTempTendAdvZ = s_PotTempTendAdvZ * DelTimeLong / TimeDisp
        s_PotTempTendSrcZ = s_PotTempTendSrcZ * DelTimeLong / TimeDisp
        s_PotTempTendTurbZ = s_PotTempTendTurbZ * DelTimeLong / TimeDisp
        s_PotTempTendRadZ = s_PotTempTendRadZ * DelTimeLong / TimeDisp
        s_PotTempTendDisZ = s_PotTempTendDisZ * DelTimeLong / TimeDisp
        s_PotTempTendSfcZ = s_PotTempTendSfcZ * DelTimeLong / TimeDisp

        f_VelZTendAdvZ = f_VelZTendAdvZ * DelTimeShort / TimeDisp
        f_VelZTendBuoyZ = f_VelZTendBuoyZ * DelTimeShort / TimeDisp
        f_VelZTendTurbZ = f_VelZTendTurbZ * DelTimeShort / TimeDisp
        f_VelZTendPressGradZ = f_VelZTendPressGradZ * DelTimeShort / TimeDisp

        call OutputDiagDataFile(hst_diag, &
             & Time,   &
             & s_PotTempTendAdvZ,    &
             & s_PotTempTendSrcZ,    &
             & s_PotTempTendTurbZ,   &
             & s_PotTempTendRadZ,    &
             & s_PotTempTendDisZ,    &
             & s_PotTempTendSfcZ,    &
             & f_VelZTendAdvZ,       &
             & f_VelZTendBuoyZ,      &
             & f_VelZTendTurbZ,      &
             & f_VelZTendPressGradZ, &
             & MassTotal,            &
             & KineticEnergyTotal,   &
             & PotentialEnergyTotal, &
             & AvailPotentialEnergyTotal, &
             & ElasticEnergyTotal,   &
             & EnergyTendRadTotal,   &
             & EnergyTendDisTotal,   &
             & EnergyTendPotTempTurbTotal,  &
             & EnergyTendVelTurbTotal,      &
             & EnergyTendSfcHeatFluxTotal,  &
             & EnergyTendSfcMomentumFluxTotal,  &
             & EnergyTendPotTempBasicTotal,  &
             & EnergyTendNonconsvTotal)

        s_PotTempTendAdvZ = 0.0d0
        s_PotTempTendSrcZ = 0.0d0
        s_PotTempTendTurbZ = 0.0d0
        s_PotTempTendRadZ = 0.0d0
        s_PotTempTendDisZ = 0.0d0
        s_PotTempTendSfcZ = 0.0d0
     end if
     
  end do
  
  !=== 出力ファイルのクローズ
  call CloseDataFile(hst_pred)
  call CloseDiagDataFile(hst_diag)

  !=== リスタートファイルの作成
  ReStartTime(1) = Time - DelTimeLong
  ReStartTime(2) = Time
  call MakeReStartFile( NewReStartFile, ReStartTime, &
       & fs_VelX_bl,    sf_VelZ_bl,  ss_Exner_bl,    &
       & ss_PotTemp_bl, ss_Km_bl,    ss_Kh_bl,       &
       & fs_VelX_nl,    sf_VelZ_nl,  ss_Exner_nl,    &
       & ss_PotTemp_nl, ss_Km_nl,    ss_Kh_nl        )


  !=== 終了時刻の記録とプログラム実行時間の書き出し
  !fortran90 の組み込み関数 date_and_time を利用.
  call date_and_time(values = EndRunTime)

  RunTime = EndRunTime - StartRunTime
  do i = 6,7
     if ( RunTime(i) < 0 ) then
        RunTime(i - 1) = RunTime(i - 1) - 1
        RunTime(i)     = RunTime(i) + 60
     end if
  end do
  write(*,*) RunTime(5),'hour',RunTime(6),'min',RunTime(7),'sec'

contains

!!!
!!!配列の割り当て, 初期化
!!!
  subroutine ArareAlloc

    call BeginSub("ArareAlloc", &
&                 fmt="%c",     &
&                 c1="Allocate and initialize allocatable variables.")

    allocate(ss_Exner_as(DimXMin:DimXMax, DimZMin:DimZMax),  &
      & fs_VelX_as(DimXMin:DimXMax, DimZMin:DimZMax)  , &
      & sf_VelZ_as(DimXMin:DimXMax, DimZMin:DimZMax), &
      & ss_Exner_ns(DimXMin:DimXMax, DimZMin:DimZMax), &
      & fs_VelX_ns(DimXMin:DimXMax, DimZMin:DimZMax), &
      & sf_VelZ_ns(DimXMin:DimXMax, DimZMin:DimZMax), &
      & fs_VelX_al(DimXMin:DimXMax, DimZMin:DimZMax), &
      & sf_VelZ_al(DimXMin:DimXMax, DimZMin:DimZMax), &
      & ss_Exner_al(DimXMin:DimXMax, DimZMin:DimZMax), &
      & ss_PotTemp_al(DimXMin:DimXMax, DimZMin:DimZMax), &
      & fs_VelX_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
      & sf_VelZ_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
      & ss_Exner_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
      & ss_PotTemp_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
      & fs_VelX_bl(DimXMin:DimXMax, DimZMin:DimZMax), &
      & sf_VelZ_bl(DimXMin:DimXMax, DimZMin:DimZMax), &
      & ss_Exner_bl(DimXMin:DimXMax, DimZMin:DimZMax), &
      & ss_PotTemp_bl(DimXMin:DimXMax, DimZMin:DimZMax), &
      & fs_AdvDiffX_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
      & sf_AdvDiffZ_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
      & sf_SfcTemp_al(DimXMin:DimXMax),      &
      & sf_SfcHeatFlux(DimXMin:DimXMax),     &
      & ff_SfcMomentumFlux(DimXMin:DimXMax),   &
      & ss_QRad(DimXMin:DimXMax, DimZMin:DimZMax), &
      & ss_QDis(DimXMin:DimXMax, DimZMin:DimZMax), &
      & ss_Km_bl(DimXMin:DimXMax, DimZMin:DimZMax), &
      & ss_Km_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
      & ss_Km_al(DimXMin:DimXMax, DimZMin:DimZMax), &
      & ss_Kh_bl(DimXMin:DimXMax, DimZMin:DimZMax), &
      & ss_Kh_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
      & ss_Kh_al(DimXMin:DimXMax, DimZMin:DimZMax) )

    fs_VelX_al = 0.0d0
    fs_VelX_nl = 0.0d0
    fs_VelX_bl = 0.0d0

    sf_VelZ_al = 0.0d0
    sf_VelZ_nl = 0.0d0
    sf_VelZ_bl = 0.0d0

    ss_Exner_al = 0.0d0
    ss_Exner_nl = 0.0d0
    ss_Exner_bl = 0.0d0
    
    ss_PotTemp_al = 0.0d0
    ss_PotTemp_nl = 0.0d0
    ss_PotTemp_bl = 0.0d0

    fs_VelX_as = 0.0d0
    fs_VelX_ns = 0.0d0

    sf_VelZ_as = 0.0d0
    sf_VelZ_ns = 0.0d0

    ss_Exner_as = 0.0d0
    ss_Exner_ns = 0.0d0

    fs_AdvDiffX_nl = 0.0d0
    sf_AdvDiffZ_nl = 0.0d0

    sf_SfcTemp_al      = 0.0d0
    sf_SfcHeatFlux     = 0.0d0 
    ff_SfcMomentumFlux = 0.0d0 
    ss_QRad = 0.0d0
    ss_QDis = 0.0d0

    ss_Km_bl = 0.0d0
    ss_Km_nl = 0.0d0
    ss_Km_al = 0.0d0
    ss_Kh_bl = 0.0d0
    ss_Kh_nl = 0.0d0
    ss_Kh_al = 0.0d0

    Time = 0.0d0

    !--- 診断量出力ファイル用に使われる変数配列の割付と初期化

    allocate(ss_PotTempTendAdv_nl(DimXMin:DimXMax, DimZMin:DimZMax),   & 
      & ss_PotTempTendSrc_nl(DimXMin:DimXMax, DimZMin:DimZMax),   &
      & ss_PotTempTendTurb_nl(DimXMin:DimXMax, DimZMin:DimZMax),   &
      & s_PotTempTendAdvZ(DimZMin:DimZMax),   &
      & s_PotTempTendSrcZ(DimZMin:DimZMax),   &
      & s_PotTempTendTurbZ(DimZMin:DimZMax),   &
      & s_PotTempTendRadZ(DimZMin:DimZMax),  &
      & s_PotTempTendDisZ(DimZMin:DimZMax),  &
      & s_PotTempTendSfcZ(DimZMin:DimZMax),  &
      & fs_VelXTendAdv_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
      & fs_VelXTendTurb_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
      & fs_VelXTendNumDiff_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
      & sf_VelZTendPressGrad_ns(DimXMin:DimXMax, DimZMin:DimZMax), &
      & sf_VelZTendAdv_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
      & sf_VelZTendTurb_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
      & sf_VelZTendBuoy_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
      & sf_VelZTendNumDiff_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
      & f_VelZTendPressGradZ(DimZMin:DimZMax), &
      & f_VelZTendAdvZ(DimZMin:DimZMax),  &
      & f_VelZTendTurbZ(DimZMin:DimZMax), &
      & f_VelZTendBuoyZ(DimZMin:DimZMax), &
      & ss_KineticEnergy_nl(DimXmin:DimXMax, DimZMin:DimZMax),  &
      & ss_PotentialEnergy_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
      & ss_AvailPotentialEnergy_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
      & ss_ElasticEnergy_nl(DimXMin:DimXMax, DimZMin:DimZMax),  &
      & ss_EnergyTendRad_nl(DimXMin:DimXMax, DimZMin:DimZMax),  & 
      & ss_EnergyTendDis_nl(DimXMin:DimXMax, DimZMin:DimZMax),  &
      & ss_EnergyTendVelTurb_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
      & ss_EnergyTendPotTempTurb_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
      & ss_EnergyTendSfcHeatFlux_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
      & ss_EnergyTendNonconsv_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
      & ss_EnergyTendPotTempBasic_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
      & ss_Z(DimXMin:DimXMax, DimZMin:DimZMax) )

    ss_PotTempTendAdv_nl = 0.0d0
    ss_PotTempTendSrc_nl = 0.0d0
    ss_PotTempTendTurb_nl = 0.0d0
    s_PotTempTendAdvZ = 0.0d0
    s_PotTempTendSrcZ = 0.0d0
    s_PotTempTendTurbZ = 0.0d0
    s_PotTempTendRadZ = 0.0d0
    s_PotTempTendDisZ = 0.0d0
    s_PotTempTendSfcZ = 0.0d0

    fs_VelXTendAdv_nl       = 0.0d0
    fs_VelXTendTurb_nl      = 0.0d0
    fs_VelXTendNumDiff_nl   = 0.0d0
    sf_VelZTendPressGrad_ns = 0.0d0
    sf_VelZTendAdv_nl       = 0.0d0
    sf_VelZTendTurb_nl      = 0.0d0
    sf_VelZTendBuoy_nl      = 0.0d0
    sf_VelZTendNumDiff_nl   = 0.0d0
    f_VelZTendPressGradZ    = 0.0d0
    f_VelZTendAdvZ          = 0.0d0
    f_VelZTendTurbZ         = 0.0d0
    f_VelZTendBuoyZ         = 0.0d0 

    ss_KineticEnergy_nl = 0.0d0
    ss_PotentialEnergy_nl = 0.0d0
    ss_AvailPotentialEnergy_nl = 0.0d0
    ss_ElasticEnergy_nl = 0.0d0
    ss_EnergyTendRad_nl = 0.0d0
    ss_EnergyTendDis_nl = 0.0d0
    ss_EnergyTendVelTurb_nl = 0.0d0
    ss_EnergyTendPotTempTurb_nl = 0.0d0
    ss_EnergyTendSfcHeatFlux_nl = 0.0d0
    ss_EnergyTendNonconsv_nl = 0.0d0
    ss_EnergyTendPotTempBasic_nl = 0.0d0

    ss_Z = 0.0d0

    call EndSub("ArareAlloc")

  end subroutine ArareAlloc

end program arare

