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

program arare
                                                                 !=begin
  !== Dependency
  use gt4_history
  use dc_trace, only: BeginSub, EndSub, DbgMessage
  use dc_message    
  use debugset, only: debugset_init
  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, TimeDisp
  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_ForceVelX_nl(:,:)
  real(8), allocatable  :: sf_ForceVelZ_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(:,:)

  real(8), allocatable  :: ss_Km_bl(:,:)
  real(8), allocatable  :: ss_Km_nl(:,:)
  real(8), allocatable  :: ss_Km_al(:,:)
  real(8), allocatable  :: ss_Kh_bl(:,:)
  real(8), allocatable  :: ss_Kh_nl(:,:)
  real(8), allocatable  :: ss_Kh_al(:,:)

  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 debugset_init(cfgfile)

  !--- 格子点情報の初期化
  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 ステップ目は Euler スキーム
  Time = Time + DelTimeLong
     
  !--- 渦粘性係数, 渦拡散係数の計算
  call EddyViscosity(fs_VelX_bl, sf_VelZ_bl, ss_PotTemp_bl, &
       & ss_Km_bl, ss_Kh_bl)    
  ss_Km_al = ss_Km_bl; ss_Km_nl = ss_Km_bl
  ss_Kh_al = ss_Kh_bl; ss_Kh_nl = ss_Kh_bl
  
  !--- 長い時間ステップの移流項, 拡散項, 粘性拡散項, 数値粘性項の計算
  call ArareMomentumTimeLong( &
       & fs_VelX_nl, sf_VelZ_nl, ss_PotTemp_nl, &
       & fs_VelX_bl, sf_VelZ_bl, ss_Km_bl, &
       & fs_ForceVelX_nl, sf_ForceVelZ_nl )
  
  !--- 長い時間ステップでのスカラー量の計算
!  call ArarePotTempTimeLong( Time, &
!       & fs_VelX_nl, sf_VelZ_nl, ss_PotTemp_nl, &
!       & ss_PotTemp_bl, ss_Kh_bl, &
!       & ss_PotTemp_al )

  !--- 長い時間ステップでのスカラー量の計算(MPDATA)
  call ArarePotTempTimeLongMPDATA( Time, &
       & fs_VelX_nl, sf_VelZ_nl, ss_PotTemp_nl, &
       & ss_PotTemp_bl, ss_Kh_bl, &
       & ss_PotTemp_al )

  !--- 短いタイムステップの初期値を作る
  call ArareVariableTimeShort( Time )
     
  !--- 短いタイムステップでの積分
  do tau = 1, NstepShort / 2
     call ArareMomentumTimeShort( &
          & fs_VelX_ns, sf_VelZ_ns, ss_Exner_ns, &
          & fs_ForceVelX_nl, sf_ForceVelZ_nl,               &
          & fs_VelX_as, sf_VelZ_as, ss_Exner_as )
     
     !--- ループを回すための処置
     call ArareLoopTimeShort
  end do
  
  !--- 最終的な速度を長い時間ステップでのものとみなす
  call ArareVariableTimeLong
  
  !--- ループを回すための処置
  call ArareLoopTimeLong
  
  !--- ファイル出力
  if ( mod( Time, TimeDisp ) == 0 ) then 
     call DbgMessage(fmt="Time = %f", d=(/Time/))
     call ArareOutputRegFile
     call ArareOutputDimFile
  end if
  
  
  !--- 長い時間ステップ
  !--- 2 ステップ目以降は leap-frog スキーム
  do t = 2, NstepLong 
     Time = Time + DelTimeLong
     
     !--- 渦粘性係数, 渦拡散係数の計算
     call ArareEddyViscosity(                                &
          & 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 ArareMomentumTimeLong( &
          & fs_VelX_nl, sf_VelZ_nl, ss_PotTemp_nl, &
          & fs_VelX_bl, sf_VelZ_bl, ss_Km_bl, &
          & fs_ForceVelX_nl, sf_ForceVelZ_nl )

     !--- 長い時間ステップでのスカラー量の計算
!     call ArarePotTempTimeLong( Time, &
!          & fs_VelX_nl, sf_VelZ_nl, ss_PotTemp_nl, &
!          & ss_PotTemp_bl, ss_Kh_bl, &
!          & ss_PotTemp_al )

    !--- 長い時間ステップでのスカラー量の計算(MPDATA)
     call ArarePotTempTimeLongMPDATA( Time, &
          & fs_VelX_nl, sf_VelZ_nl, ss_PotTemp_nl, &
          & ss_PotTemp_bl, ss_Kh_bl, &
          & ss_PotTemp_al )
     
     !--- 短いタイムステップの初期値を作る
     call ArareVariableTimeShort( Time )

     !--- 短いタイムステップでの積分
     do tau = 1, NstepShort
        call ArareMomentumTimeShort( &
             & fs_VelX_ns, sf_VelZ_ns, ss_Exner_ns, &
             & fs_ForceVelX_nl, sf_ForceVelZ_nl,               &
             & fs_VelX_as, sf_VelZ_as, ss_Exner_as )
        
        !--- ループを回すための処置
        call ArareLoopTimeShort
     end do
     
     !--- 最終的な速度を長い時間ステップでのものとみなす
     call ArareVariableTimeLong
     
     !--- アッセリンの時間フィルタ
     call ArareTimeFilter

     !--- ループを回すための処置
     call ArareLoopTimeLong

     !--- ファイル出力
     if ( mod( Time, TimeDisp ) == 0 ) then 
        call DbgMessage(fmt="Time = %f", d=(/Time/))
        call MessageNotify('Message', 'arare', "Time = %f", d=(/Time/))
        call ArareOutputRegFile
        call ArareOutputDimFile
     end if
     
  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 ArareVariableTimeShort( Time )
    use timeset, only: DelTimeLong

    real(8), intent(in) :: Time

    call BeginSub("ArareVariableTimeShort", &
&                fmt="Time = %f, %c",       &
&                d=(/Time/),  &
&                c1="Initialize variables for integration by short time step.")

    if ( Time == DelTimeLong ) then 
       fs_VelX_ns  = fs_VelX_nl
       sf_VelZ_ns  = sf_VelZ_nl
       ss_Exner_ns = ss_Exner_nl 
    else
       fs_VelX_ns  = fs_VelX_bl
       sf_VelZ_ns  = sf_VelZ_bl
       ss_Exner_ns = ss_Exner_bl 
    end if

    call EndSub("ArareVariableTimeShort")

  end subroutine ArareVariableTimeShort
  
  
!!!
!!!長い時間ステップの計算でのループを回すための操作
!!!
  subroutine ArareLoopTimeLong

    call BeginSub("ArareLoopTimeLong", &
&                 "%c",                &
&         c1="Replace variables for time integration by using long time step.")

    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

    call EndSub("ArareLoopTimeLong")

  end subroutine ArareLoopTimeLong


!!!
!!!短い時間ステップの計算でのループを回すための操作
!!!
  subroutine ArareLoopTimeShort

    call BeginSub("ArareLoopTimeShort", &
&                 "%c",                 &
&        c1="Replace variables for time integration by using short time step.")

    ss_Exner_ns  = ss_Exner_as 
    fs_VelX_ns   = fs_VelX_as
    sf_VelZ_ns   = sf_VelZ_as

    call EndSub("ArareLoopTimeShort")

  end subroutine ArareLoopTimeShort
  

!!!
!!!アッセリンの時間フィルタ
!!!
  subroutine ArareTimeFilter

    call BeginSub("ArareTimeFilter", &
&                 fmt="%c",          &
&                 c1="Asselin time filter for leapfrog scheme is adapted.")

    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 EndSub("ArareTimeFilter")

  end subroutine ArareTimeFilter

  
!!!
!!!短いタイムステップでの u, w, pi の計算
!!!
  subroutine ArareMomentumTimeShort( fs_VelX_ns, sf_VelZ_ns, ss_Exner_ns, &
       & fs_ForceVelX_nl, sf_ForceVelZ_nl,                              &
       & fs_VelX_as, sf_VelZ_as, ss_Exner_as )

    use bcset, only: fs_BC, sf_BC, ss_BC

    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_ForceVelX_nl(:,:)
    real(8), intent(in)  :: sf_ForceVelZ_nl(:,:)
    real(8), intent(out) :: fs_VelX_as(:,:)
    real(8), intent(out) :: sf_VelZ_as(:,:)
    real(8), intent(out) :: ss_Exner_as(:,:)
    
    !--- 内部変数
    real(8) ::  fs_PresGradVelX_ns(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8) ::  sf_PresGradVelZ_as(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8) ::  fs_ForceVelX_ns(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8) ::  sf_ForceVelZ_as(DimXMin:DimXMax, DimZMin:DimZMax)

    call BeginSub("ArareMomentumTimeShort", & 
&                 fmt="%c",                 &
&                 c1="Time integration of acoustic wave component.")


    !--- 速度の収束
    call DivergenceVel(fs_VelX_ns, sf_VelZ_ns, ss_DivVel_ns)
    
    !--- 速度 u の計算 (圧力傾度は前進差分で評価)
    call PressureGradVelX(ss_Exner_ns, ss_DivVel_ns, &
         &               fs_PresGradVelX_ns)    
    
    fs_ForceVelX_ns = fs_PresGradVelX_ns + fs_ForceVelX_nl 

    call Integral(fs_ForceVelX_ns, DelTimeShort, fs_VelX_ns, fs_VelX_as)
    call boundary(fs_BC, fs_VelX_as)
    
    
    !--- エクスナー関数の計算 (陰解法)
    call Exner_HeVi( ss_DivVel_ns, sf_ForceVelZ_nl, fs_VelX_as, &
         &           sf_VelZ_ns, ss_Exner_ns, ss_Exner_as)
    call boundary(ss_BC, ss_Exner_as)
    

    !--- 速度 w の計算 (圧力傾度は後退差分で評価)
    call PressureGradVelZ(ss_Exner_ns, ss_Exner_as, ss_DivVel_ns, &
         &                sf_PresGradVelZ_as)

    sf_ForceVelZ_as = sf_PresGradVelZ_as + sf_ForceVelZ_nl 

    call Integral(sf_ForceVelZ_as, DelTimeShort, sf_VelZ_ns, sf_VelZ_as)
    call boundary(sf_BC, sf_VelZ_as)

    call EndSub("ArareMomentumTimeShort")

  end subroutine ArareMomentumTimeShort


!!!
!!!長いタイムステップにおいて移流項, 浮力項, 粘性拡散項, 数値粘性を評価
!!!
  subroutine ArareMomentumTimeLong( &
       & fs_VelX_nl, sf_VelZ_nl, ss_PotTemp_nl,  &
       & fs_VelX_bl, sf_VelZ_bl, ss_Km_bl, &
       & fs_ForceVelX_nl, sf_ForceVelZ_nl)
    
    use gridset, only: DimXMin, DimXMax, DimZMin, DimZMax
    use timeset, only: DelTimeLong
    use bcset,   only: ss_BC

    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)  :: fs_VelX_bl(:,:)
    real(8), intent(in)  :: sf_VelZ_bl(:,:)
    real(8), intent(in)  :: ss_Km_bl(:,:)
    real(8), intent(out) :: fs_ForceVelX_nl(:,:)
    real(8), intent(out) :: sf_ForceVelZ_nl(:,:)
    
    real(8)  :: fs_AdvVelX_nl(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)  :: sf_AdvVelZ_nl(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)  :: sf_Buoy_nl(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)  :: fs_TurbVelX_bl(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)  :: sf_TurbVelZ_bl(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)  :: fs_NumDiffVelX_bl(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)  :: sf_NumDiffVelZ_bl(DimXMin:DimXMax, DimZMin:DimZMax)

    call BeginSub("ArareMomentumTimeLong", &
&                 "%c",                    &
&                  c1="Time integration of advection component.")

    !--- 移流項
    call AdvectVelX(fs_VelX_nl, sf_VelZ_nl, fs_AdvVelX_nl)
    call AdvectVelZ(fs_VelX_nl, sf_VelZ_nl, sf_AdvVelZ_nl)

    !--- 浮力項 
    call  Buoyancy(ss_PotTemp_nl, sf_Buoy_nl)
    
    !--- 乱流拡散項 (時刻は t - Δt で評価)
    call TurbulenceVelX(fs_VelX_bl, sf_VelZ_bl, ss_Km_bl, fs_TurbVelX_bl)
    call TurbulenceVelZ(fs_VelX_bl, sf_VelZ_bl, ss_Km_bl, sf_TurbVelZ_bl)
    
    !--- 数値粘性項 (時刻は t - Δt で評価)
    call NumDiffusionVelX(fs_VelX_bl, fs_NumDiffVelX_bl)
    call NumDiffusionVelZ(sf_VelZ_bl, sf_NumDiffVelZ_bl)
    
    !--- 項をまとめて出力
    fs_ForceVelX_nl = &
         & - fs_AdvVelX_nl + fs_TurbVelX_bl + fs_NumDiffVelX_bl 
    sf_ForceVelZ_nl = &
         & - sf_AdvVelZ_nl + sf_TurbVelZ_bl + sf_NumDiffVelZ_bl + sf_Buoy_nl

    call EndSub("ArareMomentumTimeLong")
    
  end subroutine ArareMomentumTimeLong

!!!
!!!長いタイムステップでのスカラー量の計算
!!!
  subroutine ArarePotTempTimeLong( Time, &
       & fs_VelX_nl, sf_VelZ_nl, ss_Scalar_nl, &
       & ss_Scalar_bl, ss_Kh_bl,  &
       & ss_Scalar_al)
    
    use gridset, only: DimXMin, DimXMax, DimZMin, DimZMax
    use timeset, only: DelTimeLong
    use bcset,   only: ss_BC
    use differentiate, only: sf_dz_ss
    use average,       only: ss_avr_sf
    use basicset,      only: ss_PotTempBasicZ

    implicit none
    real(8), intent(in)  :: Time
    real(8), intent(in)  :: fs_VelX_nl(:,:)
    real(8), intent(in)  :: sf_VelZ_nl(:,:)
    real(8), intent(in)  :: ss_Scalar_nl(:,:)
    real(8), intent(in)  :: ss_Scalar_bl(:,:)
    real(8), intent(in)  :: ss_Kh_bl(:,:)
    real(8), intent(out) :: ss_Scalar_al(:,:)

    !--- 内部変数
    real(8)  :: ss_AdvScalar_nl(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)  :: ss_SrcScalar_nl(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)  :: ss_TurbScalar_bl(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)  :: ss_NumDiffScalar_bl(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)  :: ss_ForceScalar_nl(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)  :: DelTime

    call BeginSub("ArarePotTempTimeLong", &
&                 fmt="%c",              &
&                 c1="Time integration of scalar component.")

    !--- 移流項
    call AdvectScalar(fs_VelX_nl, sf_VelZ_nl, ss_Scalar_nl, &
         &             ss_AdvScalar_nl)

    !--- 生成項
    ss_SrcScalar_nl = ss_avr_sf( sf_VelZ_nl * sf_dz_ss( ss_PotTempBasicZ ) ) 

    !--- 乱流拡散項 (時刻は t - Δt で評価)
    call TurbulenceScalar(ss_Scalar_bl, ss_Kh_bl, ss_TurbScalar_bl)
    
    !--- 数値粘性項 (時刻は t - Δt で評価)
    call NumDiffusionScalar(ss_Scalar_bl, ss_NumDiffScalar_bl)
    
    !--- 項をまとめる
    ss_ForceScalar_nl = &
         & - ss_AdvScalar_nl  &
         & - ss_SrcScalar_nl  &
         & + ss_TurbScalar_bl & 
         & + ss_NumDiffScalar_bl 
    
    !--- 温位の時間積分
    if ( Time == DelTimeLong ) then 
       DelTime = DelTimeLong
       call Integral(ss_ForceScalar_nl, DelTime, ss_Scalar_nl, ss_Scalar_al)
    else
       DelTime = 2.0 * DelTimeLong
       call Integral(ss_ForceScalar_nl, DelTime, ss_Scalar_bl, ss_Scalar_al)
    end if
    
    !--- 境界条件
    call boundary(ss_BC, ss_Scalar_al)
    
    call EndSub("ArarePotTempTimeLong")
    
  end subroutine ArarePotTempTimeLong

!!!
!!!長いタイムステップでのスカラー量の計算(MPDATA)
!!!
  subroutine ArarePotTempTimeLongMPDATA( Time, &
       & fs_VelX_nl, sf_VelZ_nl, ss_Scalar_nl, &
       & ss_Scalar_bl, ss_Kh_bl,  &
       & ss_Scalar_al)
    
    use gridset, only: DimXMin, DimXMax, DimZMin, DimZMax
    use timeset, only: DelTimeLong
    use bcset,   only: ss_BC
    use differentiate, only: sf_dz_ss
    use average,       only: ss_avr_sf
    use basicset,      only: ss_PotTempBasicZ

    implicit none
    real(8), intent(in)  :: Time
    real(8), intent(in)  :: fs_VelX_nl(:,:)
    real(8), intent(in)  :: sf_VelZ_nl(:,:)
    real(8), intent(in)  :: ss_Scalar_nl(:,:)
    real(8), intent(in)  :: ss_Scalar_bl(:,:)
    real(8), intent(in)  :: ss_Kh_bl(:,:)
    real(8), intent(out) :: ss_Scalar_al(:,:)

    !--- 内部変数
    real(8)  :: ss_AdvScalar_nl(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)  :: ss_SrcScalar_nl(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)  :: ss_TurbScalar_bl(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)  :: ss_ForceScalar_nl(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)  :: DelTime

    call BeginSub("ArarePotTempTimeLongMPDATA", &
&                 fmt="%c",              &
&                 c1="Time integration of scalar component.")

    !--- 移流項

    if ( Time == DelTimeLong ) then 
       DelTime = DelTimeLong
       call AdvectScalar_MPDATA(fs_VelX_nl, sf_VelZ_nl, &
&                              ss_Scalar_nl, DelTime, ss_AdvScalar_nl)
    else
       DelTime = 2.0*DelTimeLong
       call AdvectScalar_MPDATA(fs_VelX_nl, sf_VelZ_nl, &
&                              ss_Scalar_bl, DelTime, ss_AdvScalar_nl)
    end if

    !--- 生成項
    ss_SrcScalar_nl = ss_avr_sf( sf_VelZ_nl * sf_dz_ss( ss_PotTempBasicZ ) ) 

    !--- 乱流拡散項 (時刻は t - Δt で評価)
    call TurbulenceScalar(ss_Scalar_bl, ss_Kh_bl, ss_TurbScalar_bl)
    
    !--- 項をまとめる
    ss_ForceScalar_nl = &
         & - ss_AdvScalar_nl  &
         & - ss_SrcScalar_nl  &
         & + ss_TurbScalar_bl  
    
    !--- 温位の時間積分
    if ( Time == DelTimeLong ) then 
       DelTime = DelTimeLong
       call Integral(ss_ForceScalar_nl, DelTime, ss_Scalar_nl, ss_Scalar_al)
    else
       DelTime = 2.0 * DelTimeLong
       call Integral(ss_ForceScalar_nl, DelTime, ss_Scalar_bl, ss_Scalar_al)
    end if
    
    !--- 境界条件
    call boundary(ss_BC, ss_Scalar_al)
    
    call EndSub("ArarePotTempTimeLongMPDATA")
    
  end subroutine ArarePotTempTimeLongMPDATA


!!!
!!!乱流エネルギーの時間発展から乱流拡散係数を計算
!!!
  subroutine ArareEddyViscosity( &
       & 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 )
    
    use gridset,  only: DimXMin, DimXMax, DimZMin, DimZMax
    use timeset,  only: DelTimeLong
    use bcset,    only: ss_BC

    implicit none
    real(8), intent(in)  :: fs_VelX_nl(:,:)
    real(8), intent(in)  :: sf_VelZ_nl(:,:)
    real(8), intent(in)  :: ss_Km_nl(:,:)
    real(8), intent(in)  :: fs_VelX_bl(:,:)
    real(8), intent(in)  :: sf_VelZ_bl(:,:)
    real(8), intent(in)  :: ss_PotTemp_bl(:,:)
    real(8), intent(in)  :: ss_Km_bl(:,:)
    real(8), intent(out) :: ss_Km_al(:,:)
    real(8), intent(out) :: ss_Kh_al(:,:)
    
    real(8)  :: ss_AdvKm_nl(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)  :: ss_BuoyKm_bl(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)  :: ss_DivVel_bl(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)  :: ss_ShearKm_bl(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)  :: ss_DiffKm_bl(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)  :: ss_DispKm_bl(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)  :: ss_ForceKm_nl(DimXMin:DimXMax, DimZMin:DimZMax)
    real(8)  :: DelTime
    

    !--- 移流項
    if ( Time == DelTimeLong ) then 
       DelTime = DelTimeLong
       call AdvectScalar_MPDATA(fs_VelX_nl, sf_VelZ_nl, &
&                              ss_Km_nl, DelTime, ss_AdvKm_nl)
    else
       DelTime = 2*DelTimeLong
       call AdvectScalar_MPDATA(fs_VelX_nl, sf_VelZ_nl, &
&                              ss_Km_bl, DelTime, ss_AdvKm_nl)
    end if

    !--- 浮力による乱流エネルギー生成
    call SubGridBuoy(ss_PotTemp_bl, ss_BuoyKm_bl) 

    !--- シアーによる乱流エネルギー生成
    call DivergenceVel( fs_VelX_bl, sf_VelZ_bl, ss_DivVel_bl )
    call SubGridShear( fs_VelX_bl, sf_VelZ_bl, ss_DivVel_bl, ss_Km_bl, &
         &             ss_ShearKm_bl )
    
    !--- 乱流エネルギーの拡散
    call SubGridDiff( ss_Km_bl, ss_DiffKm_bl )
       
    !--- 乱流エネルギーの消散
    call SubGridDisp( ss_Km_bl, ss_DispKm_bl )    

    !--- 項をまとめる
    ss_ForceKm_nl = - ss_AdvKm_nl  - ss_BuoyKm_bl  +  ss_ShearKm_bl &
         &          + ss_DiffKm_bl - ss_DispKm_bl 
    
    !--- 時間積分
    if ( Time == DelTimeLong ) then 
       DelTime = DelTimeLong
       call Integral(ss_ForceKm_nl, DelTime, ss_Km_nl, ss_Km_al)
    else
       DelTime = 2.0 * DelTimeLong
       call Integral(ss_ForceKm_nl, DelTime, ss_Km_bl, ss_Km_al)
    end if
    DelTime = 2.0 * DelTimeLong


    !--- 境界条件
    call boundary(ss_BC, ss_Km_al)

    !--- 負の領域を零に固定
    where ( ss_Km_al < 0.0d0 )
       ss_Km_al = 0.0d0
    end where
    
    !--- 渦拡散定数を決める
    ss_Kh_al = 3.0d0 * ss_Km_al
    
  end subroutine ArareEddyViscosity


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

    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_ForceVelX_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
         & sf_ForceVelZ_nl(DimXMin:DimXMax, DimZMin:DimZMax), &
         & ss_DivVel_ns(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_ForceVelX_nl = 0.0d0
    sf_ForceVelZ_nl = 0.0d0

    ss_DivVel_ns = 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

    !速度場の初期値
    fs_VelX_bl = fs_VelXDistrb
    fs_VelX_nl = fs_VelXDistrb
    sf_VelZ_bl = sf_VelZDistrb
    sf_VelZ_nl = sf_VelZDistrb

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

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

    call EndSub("ArareAlloc")

  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 BeginSub("ArareOpneRegFile", &
&                 fmt="%c",           &
&                 c1="Initialize output file name and variables in physical region.")
    
    ! ヒストリー作成
    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)    
    
    ! 渦粘性係数
    call HistoryAddVariable(                            &
         & varname='Km', dims=(/'x','z','t'/),             &
         & longname='Km',         &
         & units='1', xtype='double',                & 
         & history = ss_reg)    
    
    ! 渦拡散係数
    call HistoryAddVariable(                            &
         & varname='Kh', dims=(/'x','z','t'/),             &
         & longname='Kh',         &
         & units='1', xtype='double',                & 
         & history = ss_reg)    

    call EndSub("ArareOpneRegFile")
    
  end subroutine ArareOpenRegFile


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

    call BeginSub("ArareOutputRegFile", &
&                 fmt="%c",             &
&                 c1="Output variables in physical region.")

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

    call HistoryPut('Exner', &
         & ss_Exner_nl(RegXMin+1:RegXMax, RegZMin+1:RegZMax), ss_reg)
    call HistoryPut('PotTemp', &
         & ss_PotTemp_nl(RegXMin+1:RegXMax, RegZMin+1:RegZMax), ss_reg)
    call HistoryPut('Km', &
         & ss_Km_nl(RegXMin+1:RegXMax, RegZMin+1:RegZMax), ss_reg)
    call HistoryPut('Kh', &
         & ss_Kh_nl(RegXMin+1:RegXMax, RegZMin+1:RegZMax), ss_reg)
    call HistoryPut('VelX', &
         & fs_VelX_nl(RegXMin:RegXMax, RegZMin+1:RegZMax), fs_reg)
    call HistoryPut('VelZ', &
         & sf_VelZ_nl(RegXMin+1:RegXMax, RegZMin:RegZMax), sf_reg)

    call EndSub("ArareOutputRegFile")

  end subroutine ArareOutputRegFile


!!!
!!!出力ファイルのクローズ. 計算領域のみ
!!!  
  subroutine ArareCloseRegFile

    call BeginSub("ArareCloseRegFile", &
&                 fmt="%c", &
&                 c1="Close output file name for physical region.")

    call HistoryClose( ss_reg )
    call HistoryClose( fs_reg )
    call HistoryClose( sf_reg )

    call EndSub("ArareCloseRegFile")

  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

    call BeginSub("ArareOpneDimFile", &
&                  fmt="%c",         &
&     c1="Initialize output file name and variables in all computational region.")

    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)    
    ! 速度
    call HistoryAddVariable(                            &
         & varname='Km', dims=(/'x','z','t'/),        &
         & longname='Km',                &
         & units='1', xtype='double',                 & 
         & history = ss_dim)    
    ! 速度
    call HistoryAddVariable(                            &
         & varname='Kh', dims=(/'x','z','t'/),        &
         & longname='Kh',                &
         & units='1', xtype='double',                 & 
         & history = ss_dim)    

    call EndSub("ArareOpneDimFile")

  end subroutine ArareOpenDimFile


!!!
!!!ファイル出力. 長い時間ステップの値を出力. 全領域
!!!
  subroutine ArareOutputDimFile

    call BeginSub("ArareOutputDimFile", &
&                 fmt="%c",            &
&                 c1="Output variables in all computational region.")

    call HistoryPut('t', Time, ss_dim)
    call HistoryPut('t', Time, fs_dim)
    call HistoryPut('t', Time, sf_dim)

    call HistoryPut('Exner', ss_Exner_nl, ss_dim)
    call HistoryPut('PotTemp', ss_PotTemp_nl, ss_dim)
    call HistoryPut('Km', ss_Km_nl, ss_dim)
    call HistoryPut('Kh', ss_Kh_nl, ss_dim)
    call HistoryPut('VelX',  fs_VelX_nl, fs_dim)
    call HistoryPut('VelZ',  sf_VelZ_nl, sf_dim)

    call EndSub("ArareOutputDimFile")

  end subroutine ArareOutputDimFile


!!!
!!!出力ファイルのクローズ. 全領域
!!!  
  subroutine ArareCloseDimFile

    call BeginSub("ArareCloseDimFile", &
&                 fmt="%c",            &
&                  c1="Close output file name for all computational region.")

    call HistoryClose( ss_dim )
    call HistoryClose( fs_dim )
    call HistoryClose( sf_dim )

    call EndSub("ArareCloseDimFile")

  end subroutine ArareCloseDimFile


end program arare
