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

program arare
  !
  !非静力学モデル deepconv/arare. 
  !

  !モジュール読み込み
  use chemcalc

  use fileset
  use debugset
  use timeset
  use gridset
  use basicset

  use damping
  use timefilter
  use boundary
  use cflcheck

!  use RestartFileIO
  use HistoryFileIO

  use DynFunc
  use DynImpFunc

  use NumDiffusion
  use Turbulence

  use MoistAdjust
  use WarmRainPrm 
  use MoistBuoyancy

  !暗黙の型宣言禁止
  implicit none

  !内部変数
  character(80) :: cfgfile
  real(8), allocatable :: pz_VelXBl(:,:)
  real(8), allocatable :: pz_VelXNl(:,:)
  real(8), allocatable :: pz_VelXAl(:,:)
  real(8), allocatable :: pz_VelXNs(:,:)
  real(8), allocatable :: pz_VelXAs(:,:)
  real(8), allocatable :: xr_VelZBl(:,:)
  real(8), allocatable :: xr_VelZNl(:,:)
  real(8), allocatable :: xr_VelZAl(:,:)
  real(8), allocatable :: xr_VelZNs(:,:)
  real(8), allocatable :: xr_VelZAs(:,:)
  real(8), allocatable :: xz_ExnerBl(:,:)
  real(8), allocatable :: xz_ExnerNl(:,:)
  real(8), allocatable :: xz_ExnerAl(:,:)
  real(8), allocatable :: xz_ExnerNs(:,:)
  real(8), allocatable :: xz_ExnerAs(:,:)
  real(8), allocatable :: xz_PotTempBl(:,:)
  real(8), allocatable :: xz_PotTempNl(:,:)
  real(8), allocatable :: xz_PotTempAl(:,:)
  real(8), allocatable :: xz_KmBl(:,:)
  real(8), allocatable :: xz_KmNl(:,:)
  real(8), allocatable :: xz_KmAl(:,:)
  real(8), allocatable :: xz_KhBl(:,:)
  real(8), allocatable :: xz_KhNl(:,:)
  real(8), allocatable :: xz_KhAl(:,:)
  real(8), allocatable :: xz_MixRtBl(:,:,:)
  real(8), allocatable :: xz_MixRtNl(:,:,:)
  real(8), allocatable :: xz_MixRtAl(:,:,:)

  real(8), allocatable :: pz_AccelVelXNl(:,:)
  real(8), allocatable :: xr_AccelVelZNl(:,:)

  real(8)              :: Time
  real(8)              :: ReStartTime(2)
  real(8), allocatable :: DelTimeLFrog(:)
  real(8)              :: DelTimeEular
  integer, allocatable :: NStepEular(:)
  integer              :: NStepLFrog

  integer              :: t, tau
  integer              :: spc

 
  !----------------------------------------------------------------------
  ! 変数参照型モジュールの初期化
  !----------------------------------------------------------------------
  !NAMELIST ファイルの取得
  call getarg( 1, cfgfile )
  write(*,*) "Input NAMELIST file: ", cfgfile
  
  !I/O ファイル名の初期化
  !  NAMELIST ファイル名を指定し, deepconv/arare の
  !  出力ファイル名を NAMELIST から得る
  call fileset_init(cfgfile)
  
  !デバグ設定の初期化
  !  NAMELIST から情報を得て, デバッグ出力スイッチの切替えを行う.  
  call debugset_init(cfgfile)
  
  !時刻に関する設定の初期化
  !  NAMELIST から必要な情報を読み取り, 時間関連の変数の設定を行う. 
  call timeset_init(cfgfile)

  !格子点情報の初期化
  !  NAMELIST から情報を得て, 格子点を計算する
  call gridset_init(cfgfile)
  
  !基本場の情報の初期化
  !  NAMELIST から情報を得て, 基本場を設定する.
  call basicset_init(cfgfile)
  
  !内部変数の初期化. とりあえずゼロを入れて値を確定させておく. 
  call ArareAlloc
  
  !予報変数の初期化
  !  ReStartFile が設定されている場合にはファイルを読み込み, 
  !  設定されていない場合には, デフォルトの基本場と擾乱場を作る. 
!  if (InitFile /= '') then    
!    !リスタートファイルの読み込み
!    call basicset_file_init(InitFile)
!
!    !擾乱場の初期値を netCDF ファイルから取得する.
!    call ReStartFile_Get(                        &
!      & ReStartTime,                             &
!      & pz_VelXBl,     xr_VelZBl,  xz_ExnerBl,   &
!      & xz_PotTempBl,  xz_KmBl,    xz_MixRtBl,   &
!      & pz_VelXNl,     xr_VelZNl,  xz_ExnerNl,   &
!      & xz_PotTempNl,  xz_KmNl,    xz_MixRtNl  )    
!  else
    !デフォルト設定の基本場, 擾乱場を作成する. 
    call BasicEnv()
    call DisturbEnv(cfgfile,                      &
      & pz_VelXBl,   xr_VelZBl,                   &
      & xz_ExnerBl,  xz_PotTempBl,                &
      & xz_MixRtBl )
!  end if
  
  
  !----------------------------------------------------------------------
  ! パッケージ型モジュールの初期化
  !   デフォルトの値から変更する必要のあるルーチンのみ初期化
  !----------------------------------------------------------------------
  call Damping_Init( cfgfile )      !波の減衰係数の初期化
  call DynImpFunc_Init()            !陰解法の初期化
  call NumDiffusion_Init()          !数値拡散項の初期化
  call Turbulence_Init()            !乱流計算の初期化
  call MoistAdjust_Init()           !湿潤飽和調整法の初期化
  call WarmRainPrm_Init()           !暖かい雨のパラメタリゼーションの初期化
  call MoistBuoy_Init()             !分子量に対する浮力計算ルーチンの初期化

  !----------------------------------------------------------------------
  ! 時刻とループ回数の初期化 
  !----------------------------------------------------------------------
!  if (InitFile /= '') then    
!    Time = ReStartTime(2)               !リスタート開始時刻
!    NstepLFrog   = NstepLong - 1
!    NstepEular   = NstepShort 
!    DelTimeLFrog = DelTimeLong * 2.0d0 
!    DelTimeEular = DelTimeShort
!  else
    Time = 0.0d0                          !計算開始時刻
    NstepLFrog      = NstepLong
    NstepEular      = NstepShort 
    NstepEular(1)   = NstepShort * 5.0d-1 ! 1 ループ目だけ
    DelTimeLFrog    = DelTimeLong * 2.0d0     
    DelTimeLFrog(1) = DelTimeLong         ! 1 ループ目だけ
    DelTimeEular    = DelTimeShort
    
    ! 1 ループ目だけは, bl と nl の値を同じにしておく. 
    ! 1 ステップ目はオイラー法で解く必要があるが, 
    ! 1 ステップ目とそれ以外のステップを別々にコーディングしたくない為
    pz_VelXNl    = pz_VelXBl
    xr_VelZNl    = xr_VelZBl
    xz_ExnerNl   = xz_ExnerBl
    xz_PotTempNl = xz_PotTempBl
    xz_KmNl      = xz_KmBl
    xz_KhNl      = xz_KhBl
    xz_MixRtNl   = xz_MixRtBl
!  end if
  
  !----------------------------------------------------------------
  ! ファイル入出力
  !----------------------------------------------------------------
  !ファイルオープン
  call HistoryFile_Open( )

!  if (InitFile /= '') then    
!    call HistoryFile_Output( ReStartTime(1),           &
!      &  xz_CNcrNl, xz_CLcrNl, xz_EVrvNl, xz_PRrNl,&
!      &  pz_VelXBl,    xr_VelZBl,   xz_ExnerBl,  &
!      &  xz_PotTempBl, xz_KmBl,     xz_MixRtBl  )
!
!    call HistoryFile_Output( ReStartTime(2),           &
!      &  xz_CNcrNl, xz_CLcrNl, xz_EVrvNl, xz_PRrNl,&
!      &  pz_VelXNl,    xr_VelZNl,   xz_ExnerNl,  &
!      &  xz_PotTempNl, xz_KmNl,     xz_MixRtNl  )
!  end if

  !時刻 0 の場合には出力
  if ( Time == 0 ) then
    call HistoryFile_Output(  &
      &  Time,                &
      &  xz_PotTempNl,        &
      &  xz_ExnerNl,          &
      &  pz_VelXNl,           &
      &  xr_VelZNl,           &
      &  xz_MixRtNl,          &
      &  xz_KmNl              )
    call HistoryFile_OutputAnal(  &
      &  xz_PotTempNl,            &
      &  xz_ExnerNl,              &
      &  xz_MixRtNl               )
  end if
  
  !----------------------------------------------------------------------
  ! 設定のチェック
  !----------------------------------------------------------------------
  !CFL 条件のチェック
  call CFLCheckTimeShort( xz_VelSoundBasicZ )

  !----------------------------------------------------------------------
  ! 数値積分
  !----------------------------------------------------------------------
  LFrog: do t = 1, NstepLFrog

    !時刻の設定
    Time = Time + DelTimeLong
    write(*,*) Time

    !----------------------------------------------------------------
    ! 渦粘性係数, 渦拡散係数を求める.
    !----------------------------------------------------------------
    xz_KmAl =                                                      &
      & xz_KmBl                                                    &
      & + DelTimeLFrog(t)                                          &
      &   * (                                                      &
      &     + xz_AdvScalar(xz_KmNl, pz_VelXNl, xr_VelZNl)           &
      &     + xz_BuoyKm(xz_PotTempBl)                               &
!      &     + xz_MoistBuoyKm(xz_ExnerBl, xz_PotTempBl, xz_MixRtBl) &
      &     + xz_ShearKm(xz_KmBl, pz_VelXBl, xr_VelZBl)             &   
      &     + xz_NumDiffScalar(xz_KmBl)                             &
      &     + xz_DispKm(xz_KmBl)                                   &
      &    )
    
    !負の処理
    where (xz_KmAl < 0.0d0)
      xz_KmAl = 0.0d0
    end where
   
    !値の上限 (中島健介, 1994, 学位論文参照)
    where (xz_KmAl > 800.0d0)
     xz_KmAl = 800.0d0
    end where

    !境界条件
    xz_KmAl = xz_Boundary_Cyc_Sym(xz_KmAl)
   
    !アッセリンの時間フィルタ. Nl の値をフィルタリング
    xz_KmNl = AsselinFilter(xz_KmAl, xz_KmNl, xz_KmBl)
  
    !渦拡散定数を決める
    xz_KhAl = 3.0d0 * xz_KmAl

    
    !----------------------------------------------------------------
    ! 温位の移流計算.
    !----------------------------------------------------------------
    xz_PotTempAl =                                                &
      &   xz_PotTempBl                                            &
      & + DelTimeLFrog(t)                                         &
      &   * (                                                     &
      &     + xz_AdvScalar( xz_PotTempNl, pz_VelXNl, xr_VelZNl)    &
      &     + xz_TurbScalar(xz_PotTempBl, xz_KhBl)                 &
      &     + xz_NumDiffScalar(xz_PotTempBl)                       &
      &     + xz_AdvScalar( xz_PotTempBasicZ, pz_VelXNl, xr_VelZNl)&
      &     + xz_TurbScalar(xz_PotTempBasicZ, xz_KhBl)             &
      &     + xz_DispQ( xz_KmBl )                                 &
!      &     + xz_RadiQ( xz_ExnerBl )                              &
!      &     + xz_HFluxSfc(pz_VelXBl, xz_PotTempBl, xz_ExnerBl)    &
      &      )
    
    !境界条件
    xz_PotTempAl = xz_Boundary_Cyc_Sym( xz_PotTempAl )
    
    !アッセリンの時間フィルタ. Nl の値をフィルタリング
    xz_PotTempNl = AsselinFilter(xz_PotTempAl, xz_PotTempNl, xz_PotTempBl)

    
    !----------------------------------------------------------------
    ! 凝縮成分の混合比の移流計算.
    !----------------------------------------------------------------
    do spc = 1, SpcNum
      !時間積分
      xz_MixRtAl(:,:,spc) =                                               & 
        &   xz_MixRtBl(:,:,spc)                                           &
        & + DelTimeLFrog(t)                                               &
        &   * (                                                           &
        &    + xz_AdvScalar(xz_MixRtNl(:,:,spc), pz_VelXNl, xr_VelZNl)    &
        &    + xz_TurbScalar(xz_MixRtBl(:,:,spc), xz_KhBl)                &
        &    + xz_NumDiffScalar(xz_MixRtBl(:,:,spc))                      &
        &    + xz_AdvScalar(xz_MixRtBasicZ(:,:,spc), xr_VelZNl, xr_VelZNl)& 
        &    + xz_TurbScalar(xz_MixRtBasicZ(:,:,spc), xz_KhBl)            &
        &    + xz_FallRain(xz_MixRtBl(:,:,spc), spc)                      &
        &    )
      !境界条件
      xz_MixRtAl(:,:,spc) = xz_Boundary_Cyc_Sym(xz_MixRtAl(:,:,spc))
      
      !アッセリンの時間フィルタ. Nl の値をフィルタリング
      xz_MixRtNl(:,:,spc) &
        & = AsselinFilter(xz_MixRtAl(:,:,spc), xz_MixRtNl(:,:,spc), xz_MixRtBl(:,:,spc))

      !移流によって負になった部分を埋める
      call FillNegative( xz_DensBasicZ, xz_MixRtAl(:,:,spc) )
      xz_MixRtAl(:,:,spc) = xz_Boundary_Cyc_Sym(xz_MixRtAl(:,:,spc))
      
    end do
    
    
    !-------------------------------------------------------------
    ! 暖かい雨のパラメタリゼーション.
    !   雲<-->雨 の変換を行う
    !-------------------------------------------------------------
    xz_MixRtAl(:,:,:) =                                 & 
     &   xz_MixRtAl(:,:,:)                             &
      & + DelTimeLFrog(t)                               &
      &   * (                                           &
      &       + xz_Cloud2Rain( xz_MixRtAl(:,:,:) )      &
      &      + xz_Cloud2RainNH4SH( xz_MixRtAl(:,:,:) ) &
      &      )
    !境界条件
    do spc = 1, SpcNum
     xz_MixRtAl(:,:,spc) = xz_Boundary_Cyc_Sym(xz_MixRtAl(:,:,spc))
    end do
    

    !-------------------------------------------------------------
    ! 湿潤飽和調節法
    !   水蒸気<-->雲の変換を行う.
    !-------------------------------------------------------------
!    call MoistAdjustSvapPress( xz_ExnerNl, xz_PotTempAl, xz_MixRtAl(:,:,:) )
    call MoistAdjustSvapPress2( xz_ExnerNl, xz_PotTempAl, xz_MixRtAl(:,:,:) )
!    call MoistAdjustNH4SH( xz_ExnerNl, xz_PotTempAl, xz_MixRtAl(:,:,:) )
    call MoistAdjustNH4SH2( xz_ExnerNl, xz_PotTempAl, xz_MixRtAl(:,:,:) )


    !-------------------------------------------------------------
    ! 暖かい雨のパラメタリゼーション.
    !   蒸気<-->雨 の変換を行う
    !-------------------------------------------------------------
   xz_PotTempAl =                                        &
      &   xz_PotTempAl                                    &
      & + DelTimeLFrog(t)                                 &
      &   * (                                             &
      &      + xz_LatentHeatRain2Gas(                   &
      &           xz_ExnerNl, xz_PotTempAl, xz_MixRtAl    &
      &          )                                        &
!      &      + xz_LatentHeatRain2Gas(                   &
!      &           xz_ExnerNl, xz_PotTempAl, xz_MixRtAl    &
!      &          )                                        &
      &      ) / ( xz_ExnerBasicZ * CpDry )
    !境界条件
    xz_PotTempAl = xz_Boundary_Cyc_Sym( xz_PotTempAl )
    
   xz_MixRtAl =                                                  &  
      &   xz_MixRtAl                                              &
      & + DelTimeLFrog(t)                                         &
      &   * (                                                     &
      &     + xz_Rain2Gas(xz_ExnerNl, xz_PotTempAl, xz_MixRtAl) &
!      &    + xz_Rain2GasNH4SH(xz_ExnerNl, xz_PotTempAl, xz_MixRtAl) &
      &    )      
    !境界条件
    do spc = 1, SpcNum
      xz_MixRtAl(:,:,spc) = xz_Boundary_Cyc_Sym(xz_MixRtAl(:,:,spc))
    end do


    !-------------------------------------------------------------
    ! 長い時間ステップでの, 速度の移流, 拡散, 数値粘性, 浮力
    !-------------------------------------------------------------
    pz_AccelVelXNl =                                        &
      & + pz_AdvVelX(pz_VelXNl, pz_VelXNl, xr_VelZNl)   &
      & + pz_TurbVelX(xz_KmBl,   pz_VelXBl, xr_VelZBl)  &
      & + pz_NumDiffVelX(pz_VelXBl)
    
    xr_AccelVelZNl =                                        &
      & + xr_AdvVelZ(xr_VelZNl, pz_VelXNl, xr_VelZNl)   &
      & + xr_Buoy(xz_PotTempNl)                         &    
      & + xr_MoistBuoy(xz_MixRtNl)                     &
      & + xr_TurbVelZ(xz_KmBl,   pz_VelXBl, xr_VelZBl)  &
      & + xr_NumDiffVelZ(xr_VelZBl)                        


    !-------------------------------------------------------------
    ! 短い時間ステップの初期値作成
    !-------------------------------------------------------------
    pz_VelXNs  = pz_VelXBl
    xr_VelZNs  = xr_VelZBl
    xz_ExnerNs = xz_ExnerBl
    
    !-------------------------------------------------------------
    ! 短い時間ステップの積分. オイラー法を利用
    !-------------------------------------------------------------
    Eular: do tau = 1, NstepEular(t)

      !-------------------------------------------------------------
      ! 速度 u の計算
      !-------------------------------------------------------------
      pz_VelXAs = &
        & pz_VelXNs                                          &
        &  + DelTimeEular                                    &
        &    * (                                             &
        &     - pz_GradPi(xz_ExnerNs, pz_VelXNs, xr_VelZNs)  &
        &     + pz_AccelVelXNl                                   &
!        &     + pz_MFluxSfc(pz_VelXNs)                       &
        &     )
      pz_VelXAs = pz_Boundary_Cyc_Sym(pz_VelXAs)

      !-------------------------------------------------------------
      ! エクスナー関数の計算
      !-------------------------------------------------------------
      xz_ExnerAs = xz_Exner( &
        & xr_AccelVelZNl,     &
        & pz_VelXNs,      &
        & pz_VelXAs,      &
        & xr_VelZNs,      &
        & xz_ExnerNs)
      xz_ExnerAs = xz_Boundary_Cyc_Sym(xz_ExnerAs)

      !-------------------------------------------------------------
      ! 速度 w の計算
      !-------------------------------------------------------------
      xr_VelZAs =                                                    &
        & xr_VelZNs                                                  &
        &  + DelTimeEular                                            &
        &    * (                                                     &
        &     - xr_GradPi(xz_ExnerAs,xz_ExnerNs,pz_VelXNs,xr_VelZNs) &
        &     + xr_AccelVelZNl                                           &
        &     )
      xr_VelZAs = xr_Boundary_Cyc_AntiSym(xr_VelZAs)
      
      !-------------------------------------------------------------
      ! 短い時間ステップのループを回すための処置
      !-------------------------------------------------------------
      xz_ExnerNs  = xz_ExnerAs
      pz_VelXNs   = pz_VelXAs
      xr_VelZNs   = xr_VelZAs

    end do Eular
    
    !----------------------------------------------------------------
    ! 最終的な短い時間ステップでの値を長い時間ステップでの値とみなす
    !----------------------------------------------------------------
    xz_ExnerAl  = xz_ExnerAs
    pz_VelXAl   = pz_VelXAs
    xr_VelZAl   = xr_VelZAs

    !-------------------------------------------------------------
    ! スポンジ層
    !-------------------------------------------------------------
    xz_ExnerAs   = xz_ExnerAs   - DelTimeLFrog(t) * xz_DampSponge_xz(xz_ExnerBl)
    pz_VelXAs    = pz_VelXAs    - DelTimeLFrog(t) * pz_DampSponge_pz(pz_VelXBl )
    xr_VelZAs    = xr_VelZAs    - DelTimeLFrog(t) * xr_DampSponge_xr(xr_VelZBl )
    xz_PotTempAl = xz_PotTempAl - DelTimeLFrog(t) * xz_DampSponge_xz(xz_PotTempBl)
        

    !----------------------------------------------------------------
    ! 長い時間ステップのループを回すための処置
    !----------------------------------------------------------------
    xz_PotTempBl = xz_PotTempNl
    xz_MixRtBl   = xz_MixRtNl
    xz_ExnerBl   = xz_ExnerNl
    pz_VelXBl    = pz_VelXNl
    xr_VelZBl    = xr_VelZNl
    xz_KmBl      = xz_KmNl
    xz_KhBl      = xz_KhNl

    xz_PotTempNl = xz_PotTempAl
    xz_MixRtNl   = xz_MixRtAl
    xz_ExnerNl   = xz_ExnerAl
    pz_VelXNl    = pz_VelXAl
    xr_VelZNl    = xr_VelZAl
    xz_KmNl      = xz_KmAl
    xz_KhNl      = xz_KhAl
    
    !----------------------------------------------------------------
    ! ファイル出力
    !----------------------------------------------------------------
    if ( mod( Time, TimeDisp ) == 0 ) then

      ! CFL 条件のチェック
      call CFLCheckTimeLongVelX( pz_VelXNl )
      call CFLCheckTimeLongVelZ( xr_VelZNl )

      !ヒストリファイル出力
      call HistoryFile_Output(  &
        &  Time,                &
        &  xz_PotTempNl,        &
        &  xz_ExnerNl,          &
        &  pz_VelXNl,           &
        &  xr_VelZNl,           &
        &  xz_MixRtNl,          &
        &  xz_KmNl              )
      call HistoryFile_OutputAnal(  &
        &  xz_PotTempNl,            &
        &  xz_ExnerNl,              &
        &  xz_MixRtNl               )
    end if
    
  end do LFrog

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

  !----------------------------------------------------------------
  ! リスタートファイルの作成
  !----------------------------------------------------------------
!  call ReStartFile_Open( )
!  call ReStartFile_OutPut(                       &
!    & Time,                                      &
!    & pz_VelXBl,    xr_VelZBl,    xz_ExnerBl,    &
!    & xz_PotTempBl, xz_KmBl,      xz_MixRtBl     )
!  call ReStartFile_OutPut(                       &
!    & Time,                                      &
!    & pz_VelXNl,    xr_VelZNl,    xz_ExnerNl,    &
!    & xz_PotTempNl, xz_KmNl,      xz_MixRtNl     )
!  call ReStartFile_Close
  
contains

  subroutine ArareAlloc
    !
    !初期化として, 配列を定義し, 値としてゼロを代入する.
    !

    !暗黙の型宣言禁止
    implicit none

    !配列割り当て
    allocate(                                                   &
      & pz_VelXBl(DimXMin:DimXMax, DimZMin:DimZMax),           &
      & pz_VelXNl(DimXMin:DimXMax, DimZMin:DimZMax),           &
      & pz_VelXAl(DimXMin:DimXMax, DimZMin:DimZMax),           &
      & pz_VelXNs(DimXMin:DimXMax, DimZMin:DimZMax),           &
      & pz_VelXAs(DimXMin:DimXMax, DimZMin:DimZMax),           &
!
      & xr_VelZBl(DimXMin:DimXMax, DimZMin:DimZMax),           &
      & xr_VelZNl(DimXMin:DimXMax, DimZMin:DimZMax),           &
      & xr_VelZAl(DimXMin:DimXMax, DimZMin:DimZMax),           &
      & xr_VelZNs(DimXMin:DimXMax, DimZMin:DimZMax),           &
      & xr_VelZAs(DimXMin:DimXMax, DimZMin:DimZMax),           &
!
      & xz_ExnerBl(DimXMin:DimXMax, DimZMin:DimZMax),          &
      & xz_ExnerNl(DimXMin:DimXMax, DimZMin:DimZMax),          &
      & xz_ExnerAl(DimXMin:DimXMax, DimZMin:DimZMax),          &
      & xz_ExnerNs(DimXMin:DimXMax, DimZMin:DimZMax),          &
      & xz_ExnerAs(DimXMin:DimXMax, DimZMin:DimZMax),          &
!
      & xz_PotTempBl(DimXMin:DimXMax, DimZMin:DimZMax),        &
      & xz_PotTempNl(DimXMin:DimXMax, DimZMin:DimZMax),        &
      & xz_PotTempAl(DimXMin:DimXMax, DimZMin:DimZMax),        &
!
      & xz_KmBl(DimXMin:DimXMax, DimZMin:DimZMax),             &
      & xz_KmNl(DimXMin:DimXMax, DimZMin:DimZMax),             &
      & xz_KmAl(DimXMin:DimXMax, DimZMin:DimZMax),             &
!
      & xz_KhBl(DimXMin:DimXMax, DimZMin:DimZMax),             &
      & xz_KhNl(DimXMin:DimXMax, DimZMin:DimZMax),             &
      & xz_KhAl(DimXMin:DimXMax, DimZMin:DimZMax),             &
!
      & xz_MixRtBl(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum), &
      & xz_MixRtNl(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum), &
      & xz_MixRtAl(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum), &
!
      & pz_AccelVelXNl(DimXMin:DimXMax, DimZMin:DimZMax),        &
      & xr_AccelVelZNl(DimXMin:DimXMax, DimZMin:DimZMax),        & 
!
      & DelTimeLFrog(NstepLong), NStepEular(NStepLong)               ) 

    pz_VelXNs = 0.0d0;     pz_VelXAs = 0.0d0    
    xr_VelZNs = 0.0d0;     xr_VelZAs = 0.0d0
    xz_ExnerNs = 0.0d0;    xz_ExnerAs = 0.0d0

    pz_VelXBl = 0.0d0;     pz_VelXNl = 0.0d0;     pz_VelXAl = 0.0d0
    xr_VelZBl = 0.0d0;     xr_VelZNl = 0.0d0;     xr_VelZAl = 0.0d0
    xz_ExnerBl = 0.0d0;    xz_ExnerNl = 0.0d0;    xz_ExnerAl = 0.0d0
    xz_KmBl = 0.0d0;       xz_KmNl = 0.0d0;       xz_KmAl = 0.0d0
    xz_KhBl = 0.0d0;       xz_KhNl = 0.0d0;       xz_KhAl = 0.0d0
    xz_PotTempBl = 0.0d0;  xz_PotTempNl = 0.0d0;  xz_PotTempAl = 0.0d0
    xz_MixRtBl = 0.0d0;    xz_MixRtNl = 0.0d0;    xz_MixRtAl = 0.0d0

    pz_AccelVelXNl = 0.0d0
    xr_AccelVelZNl = 0.0d0

    DelTimeEular = 0.0d0
    DelTimeLFrog = 0.0d0; 
    NStepEular = 0.0d0
    NstepLFrog = 0.0d0

  end subroutine ArareAlloc
      

end program arare
