!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2004. All rights reserved.
!---------------------------------------------------------------------
                                                                 !=begin
!= Program Arare
!
!   * Developer: SUGIYAMA Ko-ichiro (sugiyama@gfd-dennou.org)
!   * Version: $Id$ 
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview 
!
!非静力学モデル deepconv/arare. 
!
!== Error Handling
!
!== Known Bugs
!
!== Note
!
!方程式系は準圧縮系.
!
!== Future Plans
!
                                                                 !=end

program arare
                                                                 !=begin
  !== Dependency
  use gt4_history
  use gridset,  only: gridset_init, DimXMin, DimXMax, DimZMin, DimZMax, &
       &              RegZMin, RegZMax
  use fileset,  only: fileset_init, cfgfile, basicfile, distrbfile
  use timeset,  only: timeset_init, NstepShort, NstepLong, &
       &              DelTimeShort, DelTimeLong
  use expname,  only: expname_init
  use linlib,   only: linlib_init
  use bcset,    only: bcset_init
  use basicset, only: basicset_init
  use distrbset,only: distrbset_init
  use arareprm, only: arareprm_init
                                                                 !=end
  implicit none

  real(8), allocatable  :: ss_Exner_as(:,:)
  real(8), allocatable  :: fs_VelX_as(:,:)
  real(8), allocatable  :: sf_VelZ_as(:,:)
  real(8), allocatable  :: ss_Exner_ns(:,:)
  real(8), allocatable  :: fs_VelX_ns(:,:)
  real(8), allocatable  :: sf_VelZ_ns(:,:)

  real(8), allocatable  :: ss_DivVel_ns(:,:)
  real(8), allocatable  :: fs_ForceX_nl(:,:)
  real(8), allocatable  :: sf_ForceZ_nl(:,:)
  
  real(8), allocatable  :: ss_PotTemp_bl(:,:)
  real(8), allocatable  :: ss_Exner_bl(:,:)
  real(8), allocatable  :: fs_VelX_bl(:,:)
  real(8), allocatable  :: sf_VelZ_bl(:,:)
  real(8), allocatable  :: ss_PotTemp_nl(:,:)
  real(8), allocatable  :: ss_Exner_nl(:,:)
  real(8), allocatable  :: fs_VelX_nl(:,:)
  real(8), allocatable  :: sf_VelZ_nl(:,:)
  real(8), allocatable  :: ss_PotTemp_al(:,:)
  real(8), allocatable  :: ss_Exner_al(:,:)
  real(8), allocatable  :: fs_VelX_al(:,:)
  real(8), allocatable  :: sf_VelZ_al(:,:)

  !--- test
  real(8), allocatable  :: ss_Km_al(:,:)
  real(8), allocatable  :: ss_Kh_al(:,:)
  real(8), allocatable  :: fs_DiffuseX_nl(:,:)
  real(8), allocatable  :: sf_DiffuseZ_nl(:,:)

  real(8)               :: Time
  integer               :: t, tau

  type(gt_history)      :: ss_reg, fs_reg, sf_reg
  type(gt_history)      :: ss_dim, fs_dim, sf_dim
  

!!!
!!!初期化
!!!
  !--- I/O ファイル名の初期化
  call fileset_init
  
  !--- 格子点情報の初期化
  call gridset_init(cfgfile)
  
  !--- 時間積分の設定の初期化
  call timeset_init(cfgfile)
  
  !--- 実験名の初期化
  call expname_init(cfgfile)
  
  !--- 計算の設定を読み込む
  call arareprm_init(cfgfile)

  !--- 境界条件を設定
  call bcset_init(cfgfile)

  !--- 基本場の値を設定
  call basicset_init(basicfile)

  !--- 擾乱場の値を設定
  call distrbset_init(distrbfile)
  
  !--- 行列計算ライブラリの初期化. 
  !--- 引数はスカラーの z 方向の計算領域の大きさ
  call linlib_init(cfgfile, RegZMax - RegZMin)
  
  !--- 出力ファイルの定義
  call ArareOpenRegFile
  call ArareOpenDimFile

!!!
!!! CFL 条件のチェック
!!!
  call cflchk


!!!
!!!変数の初期化
!!!  
  !--- 配列の割り当て, 初期化
  call ArareAlloc


!!!
!!!数値積分
!!!
  !--- 初期状態をファイル出力
  Time = 0.0d0
  call ArareOutputRegFile
  call ArareOutputDimFile
  
  !--- 長い時間ステップ
  !--- 最初の 1 ステップだけは前進差分
  Time = Time + DelTimeLong
  
  !長い時間ステップの移流項, 拡散項, 温位の計算
  call ArareNonAcousticTerm1st( fs_VelX_nl, sf_VelZ_nl, ss_PotTemp_nl, &
       & fs_ForceX_nl, sf_ForceZ_nl, ss_PotTemp_al) 
  
  !短いタイムステップの初期値
  fs_VelX_ns  = fs_VelX_nl
  sf_VelZ_ns  = sf_VelZ_nl
  ss_Exner_ns = ss_Exner_nl 
  
  !--- 短いタイムステップでの積分
  !--- NstepShort は偶数である必要あり
  do tau = 1, NstepShort / 2.0d0
     
     call ArareAcousticTerm( &
          & fs_VelX_ns, sf_VelZ_ns, ss_Exner_ns,  &
          & fs_ForceX_nl, sf_ForceZ_nl,             & 
          & fs_VelX_as, sf_VelZ_as, ss_Exner_as )

     !--- ループを回すための処置
     call ArareLoopTimeShort
  end do

  !--- 最終的な速度を長い時間ステップでのものとみなす
  call ArareVariableTimeLong

  !--- アッセリンの時間フィルタ
  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 )

  !--- 粘性拡散係数
  call TurbStab(fs_VelX_al, sf_VelZ_al, ss_PotTemp_al, ss_Km_al, ss_Kh_al)
  
  !--- ファイル出力
  call ArareOutputRegFile
  call ArareOutputDimFile
  
  !--- ループを回すための処置
  call ArareLoopTimeLong
  
  
  !--- 長い時間ステップ
  !--- 2 ステップ目以降は leap-frog スキーム
  do t = 1, NstepLong 
     Time = Time + DelTimeLong
     write(*,*) Time
     
     !長い時間ステップの移流項, 拡散項, 拡散項の計算
     call ArareNonAcousticTerm( fs_VelX_nl, sf_VelZ_nl, &
          & ss_PotTemp_nl, ss_PotTemp_bl, &
          & fs_ForceX_nl, sf_ForceZ_nl, ss_PotTemp_al )
  
     !短いタイムステップの初期値
     fs_VelX_ns  = fs_VelX_bl
     sf_VelZ_ns  = sf_VelZ_bl
     ss_Exner_ns = ss_Exner_bl 
     
     !--- 短いタイムステップでの積分
     do tau = 1, NstepShort
        
        call ArareAcousticTerm( &
             & fs_VelX_ns, sf_VelZ_ns, ss_Exner_ns, &
             & fs_ForceX_nl, sf_ForceZ_nl,               &
             & fs_VelX_as, sf_VelZ_as, ss_Exner_as )
        
        !--- ループを回すための処置
        call ArareLoopTimeShort
     end do
     
     !--- 最終的な速度を長い時間ステップでのものとみなす
     call ArareVariableTimeLong

     !--- test
     !--- 音波に関連しない項 
     call DiffuseX(fs_VelX_nl, sf_VelZ_nl, ss_PotTemp_nl, fs_DiffuseX_nl)
     call DiffuseZ(fs_VelX_nl, sf_VelZ_nl, ss_PotTemp_nl, sf_DiffuseZ_nl)

     !--- ファイル出力
     call ArareOutputRegFile
     call ArareOutputDimFile
     
     !--- ループを回すための処置
     call ArareLoopTimeLong
  end do

  !--- 出力ファイルのクローズ
  call ArareCloseRegFile
  call ArareCloseDimFile


contains

!!!
!!!短いタイムステップで計算した値を長いタイムステップでの値に代入
!!!
  subroutine ArareVariableTimeLong
    ss_Exner_al  = ss_Exner_as 
    fs_VelX_al   = fs_VelX_as
    sf_VelZ_al   = sf_VelZ_as     
  end subroutine ArareVariableTimeLong
  
  
!!!
!!!長い時間ステップの計算でのループを回すための操作
!!!
  subroutine ArareLoopTimeLong
    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_PotTemp_nl = ss_PotTemp_al
    ss_Exner_nl  = ss_Exner_al 
    fs_VelX_nl   = fs_VelX_al
    sf_VelZ_nl   = sf_VelZ_al
  end subroutine ArareLoopTimeLong


!!!
!!!短い時間ステップの計算でのループを回すための操作
!!!
  subroutine ArareLoopTimeShort
    ss_Exner_ns  = ss_Exner_as 
    fs_VelX_ns = fs_VelX_as
    sf_VelZ_ns = sf_VelZ_as
  end subroutine ArareLoopTimeShort
  
  
!!!
!!!短いタイムステップでの u, w, pi の計算
!!!
  subroutine ArareAcousticTerm( fs_VelX_ns, sf_VelZ_ns, ss_Exner_ns, &
       & fs_ForceX_nl, sf_ForceZ_nl,                              &
       & fs_VelX_as, sf_VelZ_as, ss_Exner_as )

    implicit none
    real(8), intent(in)  :: fs_VelX_ns(:,:)
    real(8), intent(in)  :: sf_VelZ_ns(:,:)
    real(8), intent(in)  :: ss_Exner_ns(:,:)
    real(8), intent(in)  :: fs_ForceX_nl(:,:)
    real(8), intent(in)  :: sf_ForceZ_nl(:,:)
    real(8), intent(out) :: fs_VelX_as(:,:)
    real(8), intent(out) :: sf_VelZ_as(:,:)
    real(8), intent(out) :: ss_Exner_as(:,:)

    !--- 速度の収束
    call DivVel(fs_VelX_ns, sf_VelZ_ns, ss_DivVel_ns)
     
    !--- 速度 u の計算
    call VelX_Euler(ss_Exner_ns, ss_DivVel_ns, fs_ForceX_nl, fs_VelX_ns, &
         & fs_VelX_as)
    
    !--- エクスナー関数の計算
    call Exner_HeVi(   &
         & ss_DivVel_ns, sf_ForceZ_nl, fs_VelX_as, sf_VelZ_ns, ss_Exner_ns, &
         & ss_Exner_as)
    
    !--- 速度 w の計算
    call VelZ_Backward(    &
         & ss_Exner_ns, ss_Exner_as, ss_DivVel_ns, &
         & sf_ForceZ_nl, sf_VelZ_ns, &
         & sf_VelZ_as)

  end subroutine ArareAcousticTerm


!!!
!!!長いタイムステップでの u, w, pi の計算
!!!
  subroutine ArareNonAcousticTerm1st( fs_VelX_nl, sf_VelZ_nl, &
       & ss_PotTemp_nl, fs_ForceX_nl, sf_ForceZ_nl, ss_PotTemp_al)
    use gridset, only: DimXMin, DimXMax, DimZMin, DimZMax

    implicit none
    real(8), intent(in)  :: fs_VelX_nl(:,:)
    real(8), intent(in)  :: sf_VelZ_nl(:,:)
    real(8), intent(in)  :: ss_PotTemp_nl(:,:)
    real(8), intent(out) :: fs_ForceX_nl(:,:)
    real(8), intent(out) :: sf_ForceZ_nl(:,:)
    real(8), intent(out) :: ss_PotTemp_al(:,:)

    real(8)  ::  fs_AdvX_nl(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)  ::  sf_AdvZ_nl(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)  ::  sf_Buoy_nl(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)  ::  fs_DiffuseX_nl(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)  ::  sf_DiffuseZ_nl(DimXMin:DimXMax, DimZMin:DimZMax)


    !--- 移流項
    call AdvectX(fs_VelX_nl, sf_VelZ_nl, fs_AdvX_nl)
    call AdvectZ(fs_VelX_nl, sf_VelZ_nl, sf_AdvZ_nl)

    !--- 浮力項
    call  Buoyancy(ss_PotTemp_nl, sf_Buoy_nl)
    
    !--- 音波に関連しない項 
    call DiffuseX(fs_VelX_nl, sf_VelZ_nl, ss_PotTemp_nl, fs_DiffuseX_nl)
    call DiffuseZ(fs_VelX_nl, sf_VelZ_nl, ss_PotTemp_nl, sf_DiffuseZ_nl)
    
    !--- 項をまとめる
    fs_ForceX_nl = - fs_AdvX_nl + fs_DiffuseX_nl
    sf_ForceZ_nl = - sf_AdvZ_nl + sf_DiffuseZ_nl + sf_Buoy_nl

    !温位の数値積分. Euler 法を利用
    call PotTemp_Euler(fs_VelX_nl, sf_VelZ_nl, ss_PotTemp_nl, ss_PotTemp_al)
    
  end subroutine ArareNonAcousticTerm1st


!!!
!!!長いタイムステップでの u, w, pi の計算
!!!
  subroutine ArareNonAcousticTerm( fs_VelX_nl, sf_VelZ_nl, &
       & ss_PotTemp_nl, ss_PotTemp_bl, &
       & fs_ForceX_nl, sf_ForceZ_nl, ss_PotTemp_al)
    use gridset, only: DimXMin, DimXMax, DimZMin, DimZMax

    implicit none
    real(8), intent(in)  :: fs_VelX_nl(:,:)
    real(8), intent(in)  :: sf_VelZ_nl(:,:)
    real(8), intent(in)  :: ss_PotTemp_nl(:,:)
    real(8), intent(in)  :: ss_PotTemp_bl(:,:)
    real(8), intent(out) :: fs_ForceX_nl(:,:)
    real(8), intent(out) :: sf_ForceZ_nl(:,:)
    real(8), intent(out) :: ss_PotTemp_al(:,:)
    
    real(8)  ::  fs_AdvX_nl(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)  ::  sf_AdvZ_nl(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)  ::  sf_Buoy_nl(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)  ::  fs_DiffuseX_nl(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)  ::  sf_DiffuseZ_nl(DimXMin:DimXMax, DimZMin:DimZMax)


    !--- 移流項
    call AdvectX(fs_VelX_nl, sf_VelZ_nl, fs_AdvX_nl)
    call AdvectZ(fs_VelX_nl, sf_VelZ_nl, sf_AdvZ_nl)

    !--- 浮力項
    call  Buoyancy(ss_PotTemp_nl, sf_Buoy_nl)
    
    !--- 音波に関連しない項 
    call DiffuseX(fs_VelX_nl, sf_VelZ_nl, ss_PotTemp_nl, fs_DiffuseX_nl)
    call DiffuseZ(fs_VelX_nl, sf_VelZ_nl, ss_PotTemp_nl, sf_DiffuseZ_nl)
    
    !--- 項をまとめる
    fs_ForceX_nl = - fs_AdvX_nl + fs_DiffuseX_nl
    sf_ForceZ_nl = - sf_AdvZ_nl + sf_DiffuseZ_nl + sf_Buoy_nl
    
    !温位の数値積分. Leap-Frog 法を利用
    call PotTemp_LeapFrog(fs_VelX_nl, sf_VelZ_nl, &
         & ss_PotTemp_bl, ss_PotTemp_nl, &
         & ss_PotTemp_al)
    
  end subroutine ArareNonAcousticTerm


!!!
!!!配列の割り当て, 初期化
!!!
  subroutine ArareAlloc
    use distrbset, only: fs_VelXDistrb, sf_VelZDistrb, &
         &               ss_ExnerDistrb, ss_PotTempDistrb

    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_ForceX_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
         & sf_ForceZ_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_DivVel_ns(DimXMin:DimXMax, DimZMin:DimZMax), &
         !--- test
         & fs_DiffuseX_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
         & sf_DiffuseZ_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_Km_al(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_ForceX_nl = 0.0d0
    sf_ForceZ_nl = 0.0d0

    ss_DivVel_ns = 0.0d0

    Time = 0.0d0

    !速度場の初期値
    fs_VelX_al  = fs_VelXDistrb
    fs_VelX_nl  = fs_VelXDistrb
    sf_VelZ_al = sf_VelZDistrb
    sf_VelZ_nl = sf_VelZDistrb

    !無次元圧力の初期値
    ss_Exner_al = ss_ExnerDistrb
    ss_Exner_nl = ss_ExnerDistrb

    !温位の初期値
    ss_PotTemp_al = ss_PotTempDistrb
    ss_PotTemp_nl = ss_PotTempDistrb

  end subroutine ArareAlloc


!!!
!!!出力ファイルの定義. 計算領域のみ
!!! 
  subroutine ArareOpenRegFile
    use gridset, only: RegXMin, RegXMax, RegZMin, RegZMax, &
         & NX, NZ, f_X, f_Z, s_X, s_Z
    use expname, only: exptitle, expsrc, expinst
    use fileset, only: fs_RegFile, sf_RegFile, ss_RegFile
    
    ! ヒストリー作成
    call HistoryCreate(                                    &
         & file = fs_RegFile,                               &
         & title = exptitle,                               &
         & source = expsrc,                                &
         & institution = expinst,                          &
         & dims=(/'x','z','t'/),                           &
         & dimsizes=(/NX+1, NZ, 0/),                       &
         & longnames=(/'X-coordinate',                     &
         &             'Z-coordinate',                     &
         &             'Time        '/),                   &
         & units=(/'m','m','s'/), origin=0.0,              &
         & interval=0.0, &
         & history = fs_reg)
    
    ! 変数出力
    call HistoryPut('x', f_X(RegXMin: RegXMax), fs_reg)
    call HistoryPut('z', s_Z(RegXMin+1: RegXMax), fs_reg)
    
    ! ヒストリー作成
    call HistoryCreate(                                    &
         & file = sf_RegFile,                               &
         & title = exptitle,                               &
         & source = expsrc,                                &
         & institution = expinst,                          &
         & dims=(/'x','z','t'/),                           &
         & dimsizes=(/NX, NZ+1, 0/),                       &
         & longnames=(/'X-coordinate',                     & 
         &             'Z-coordinate',                     &
         &             'Time        '/),                   &
         & units=(/'m','m','s'/),                          &
         & origin=0.0,                                     &
         & interval=0.0,                                   &
         & history = sf_reg)
    
    ! 変数出力
    call HistoryPut('x', s_X(RegXMin+1: RegXMax), sf_reg)
    call HistoryPut('z', f_Z(RegXMin: RegXMax), sf_reg)

    ! ヒストリー作成
    call HistoryCreate(                                  &
         & file = ss_RegFile,                            &         
         & title = exptitle,                             &
         & source = expsrc,                              &
         & institution = expinst,                        &
         & dims=(/'x','z','t'/),                         &
         & dimsizes=(/NX, NZ, 0/),                       &
         & longnames=(/'X-coordinate',                   &
         &             'Z-coordinate',                   &
         &             'Time        '/),                 &
         & units=(/'m','m','s'/), origin=0.0,            &
         & interval=0.0,                                 &
         & history = ss_reg)
    
    ! 変数出力
    call HistoryPut('x', s_X(RegXMin+1: RegXMax), ss_reg)
    call HistoryPut('z', s_Z(RegXMin+1: RegXMax), ss_reg)
    
    ! 無次元圧力
    call HistoryAddVariable(                             &
         & varname='Exner', dims=(/'x','z','t'/),        &
         & longname='nondimensional pressure', units='1',&
         & xtype='double',                               & 
         & history = ss_reg)    
    ! 仮温位
    call HistoryAddVariable(                            &
         & varname='PotTemp', dims=(/'x','z','t'/),      &
         & longname='virtual potential temperature',    &
         & units='K', xtype='double',                 & 
         & history = ss_reg)    
    ! 速度
    call HistoryAddVariable(                            &
         & varname='VelX', dims=(/'x','z','t'/),             &
         & longname='zonal velocity',         &
         & units='m/s', xtype='double',                & 
         & history = fs_reg)        
    ! 速度
    call HistoryAddVariable(                            &
         & varname='VelZ', dims=(/'x','z','t'/),             &
         & longname='vertical velocity',         &
         & units='m/s', xtype='double',                & 
         & history = sf_reg)    
    
  end subroutine ArareOpenRegFile


!!!
!!!ファイル出力. 長い時間ステップで出力. 計算領域のみ
!!!
  subroutine ArareOutputRegFile
    use gridset, only: RegXMin, RegXMax, RegZMin, RegZMax

    call HistoryPut('t', Time, ss_reg)
    call HistoryPut('t', Time, fs_reg)
    call HistoryPut('t', Time, sf_reg)

    call HistoryPut('Exner', &
         & ss_Exner_al(RegXMin+1:RegXMax, RegZMin+1:RegZMax), ss_reg)
    call HistoryPut('PotTemp', &
         & ss_PotTemp_al(RegXMin+1:RegXMax, RegZMin+1:RegZMax), ss_reg)
    call HistoryPut('VelX', &
         & fs_VelX_al(RegXMin:RegXMax, RegZMin+1:RegZMax), fs_reg)
    call HistoryPut('VelZ', &
         & sf_VelZ_al(RegXMin+1:RegXMax, RegZMin:RegZMax), sf_reg)
  end subroutine ArareOutputRegFile


!!!
!!!出力ファイルのクローズ. 計算領域のみ
!!!  
  subroutine ArareCloseRegFile
    call HistoryClose( ss_reg )
    call HistoryClose( fs_reg )
    call HistoryClose( sf_reg )
  end subroutine ArareCloseRegFile


!!!
!!!出力ファイルの定義. 全領域
!!!
  subroutine ArareOpenDimFile
    use gridset, only: f_X, f_Z, s_X, s_Z
    use expname, only: exptitle, expsrc, expinst
    use fileset, only: fs_DimFile, sf_DimFile, ss_DimFile
    
    integer::   N
    
    N = size(f_X, 1)

    ! ヒストリー作成
    call HistoryCreate(                                    &
         & file = fs_DimFile,                              &
         & title = exptitle,                               &
         & source = expsrc,                                &
         & institution = expinst,                          &
         & dims=(/'x','z','t'/),                           &
         & dimsizes=(/N, N, 0/),                           &
         & longnames=(/'X-coordinate',                     &
         &             'Z-coordinate',                     &
         &             'Time        '/),                   &
         & units=(/'m','m','s'/), origin=0.0,              &
         & interval=0.0, &
         & history = fs_dim)
    
    ! 変数出力
    call HistoryPut('x', f_X, fs_dim)
    call HistoryPut('z', s_Z, fs_dim)
    
    ! ヒストリー作成
    call HistoryCreate(                                    &
         & file = sf_DimFile,                              &
         & title = exptitle,                               &
         & source = expsrc,                                &
         & institution = expinst,                          &
         & dims=(/'x','z','t'/),                           &
         & dimsizes=(/N, N, 0/),                           &
         & longnames=(/'X-coordinate',                     & 
         &             'Z-coordinate',                     &
         &             'Time        '/),                   &
         & units=(/'m','m','s'/),                          &
         & origin=0.0,                                     &
         & interval=0.0,                                   &
         & history = sf_dim)
    
    ! 変数出力
    call HistoryPut('x', s_X, sf_dim)
    call HistoryPut('z', f_Z, sf_dim)

    ! ヒストリー作成
    call HistoryCreate(                                  &
         & file = ss_DimFile,                            &
         & title = exptitle,                             &
         & source = expsrc,                              &
         & institution = expinst,                        &
         & dims=(/'x','z','t'/),                         &
         & dimsizes=(/N, N, 0/),                         &
         & longnames=(/'X-coordinate',                   &
         &             'Z-coordinate',                   &
         &             'Time        '/),                 &
         & units=(/'m','m','s'/), origin=0.0,            &
         & interval=0.0,                                 &
         & history = ss_dim)
    
    ! 変数出力
    call HistoryPut('x', s_X, ss_dim)
    call HistoryPut('z', s_Z, ss_dim)
    
    ! 無次元圧力
    call HistoryAddVariable(                             &
         & varname='Exner', dims=(/'x','z','t'/),        &
         & longname='nondimensional pressure', units='1',&
         & xtype='double',                               & 
         & history = ss_dim)    
    ! 仮温位
    call HistoryAddVariable(                            &
         & varname='PotTemp', dims=(/'x','z','t'/),      &
         & longname='virtual potential temperature',    &
         & units='K', xtype='double',                   & 
         & history = ss_dim)    
    ! 速度
    call HistoryAddVariable(                            &
         & varname='VelX', dims=(/'x','z','t'/),        &
         & longname='zonal velocity',                   &
         & units='m/s', xtype='double',                 & 
         & history = fs_dim)        
    ! 速度
    call HistoryAddVariable(                            &
         & varname='VelZ', dims=(/'x','z','t'/),        &
         & longname='vertical velocity',                &
         & units='m/s', xtype='double',                 & 
         & history = sf_dim)    

    !--- test
    ! 粘性拡散係数
    call HistoryAddVariable(                            &
         & varname='Km', dims=(/'x','z','t'/),        &
         & longname='viscosity',                &
         & units='1', xtype='double',                 & 
         & history = ss_dim)    
    ! 拡散項
    call HistoryAddVariable(                            &
         & varname='Km', dims=(/'x','z','t'/),        &
         & longname='diffusion',                &
         & units='1', xtype='double',                 & 
         & history = fs_dim)    
    ! 拡散項
    call HistoryAddVariable(                            &
         & varname='Km', dims=(/'x','z','t'/),        &
         & longname='diffusion',                &
         & units='1', xtype='double',                 & 
         & history = sf_dim)    
  end subroutine ArareOpenDimFile


!!!
!!!ファイル出力. 長い時間ステップの値を出力. 全領域
!!!
  subroutine ArareOutputDimFile
    call HistoryPut('t', Time, ss_dim)
    call HistoryPut('t', Time, fs_dim)
    call HistoryPut('t', Time, sf_dim)

    call HistoryPut('Exner', ss_Exner_al, ss_dim)
    call HistoryPut('PotTemp', ss_PotTemp_al, ss_dim)
    call HistoryPut('VelX',  fs_VelX_al, fs_dim)
    call HistoryPut('VelZ',  sf_VelZ_al, sf_dim)

    !--- test
    call HistoryPut('Km', ss_Km_al, ss_dim)
    call HistoryPut('Km', fs_DiffuseX_nl, fs_dim)
    call HistoryPut('Km', sf_DiffuseZ_nl, sf_dim)
  end subroutine ArareOutputDimFile


!!!
!!!出力ファイルのクローズ. 全領域
!!!  
  subroutine ArareCloseDimFile
    call HistoryClose( ss_dim )
    call HistoryClose( fs_dim )
    call HistoryClose( sf_dim )
  end subroutine ArareCloseDimFile


end program arare
