!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004. All rights reserved.
!---------------------------------------------------------------------
                                                                 !=begin
!= Program Arare
!
!   * Developer: SUGIYAMA Ko-ichiro 
!   * Version: $ 
!   * 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 方向の配列の上限
       &              RegZMin,       & ! z 方向の物理領域の下限
       &              RegZMax          ! z 方向の物理領域の上限 
  use timeset,  only: timeset_init,  & ! 時間積分設定初期化サブルーチン
       &              NstepShort,    & ! 長いタイムステップのステップ数
       &              NstepLong,     & ! 短いタイムステップのステップ数
       &              DelTimeShort,  & ! 長いタイムステップ
       &              DelTimeLong,   & ! 短いタイムステップ
       &              TimeDisp         ! ファイル出力する時間間隔
  use nameset,  only: nameset_init     ! 実験名初期化サブルーチン
  use physset,  only: physset_init     ! 物理定数の初期化サブルーチン
  use linlib,   only: linlib_init      ! 行列計算ライブラリ初期化サブルーチン
  use bcset,    only: bcset_init       ! 境界条件設定サブルーチン
  use basicset, only: basicset_file_init ! 基本場設定サブルーチン
  use arareset, only: arareset_init    ! 計算設定を読み込むサブルーチン
                                                                 !=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)               :: Time               ! 計算継続時間
  real(8)               :: DelTime            ! 時間ステップ
  real(8)               :: ReStartTime(2)     ! リスタートファイルの時刻
  integer               :: t, tau             ! ループを回すための変数
  integer               :: NStep_leapfrog     ! ループを回す回数


  
!== 初期化
  
  !=== 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


!== 変数の初期化

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

  !=== 基本場の値を設定
  !基本場の値を 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( 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( 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)
     
     !=== 長い時間ステップでの速度 u の移流, 拡散, 粘性
     call AdvectDiffuseX(                     &
          & fs_VelX_nl, sf_VelZ_nl,           & !移流
          & fs_VelX_bl, sf_VelZ_bl, ss_Km_bl, & !乱流拡散, 数値粘性
          & fs_AdvDiffX_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 )
     
     !=== 短いタイムステップの初期値を作る
     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 )
        
        !==== 短い時間ステップのループを回すための処置
        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( Time, &
             &  fs_VelX_nl,    sf_VelZ_nl, ss_Exner_nl, &
             &  ss_PotTemp_nl, ss_Km_nl,   ss_Kh_nl )
     end if
     
  end do
  
  !=== 出力ファイルのクローズ
  call CloseDataFile

  !=== リスタートファイルの作成
  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        )

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), &
         & 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

    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

    call EndSub("ArareAlloc")

  end subroutine ArareAlloc

end program arare

