!----------------------------------------------------------------------
!     Copyright (c) 2005 Masatsugu Odaka. All rights reserved.
!----------------------------------------------------------------------
! Sample program for gt4f90io and ISPACK  2005/01/25 M. Odaka
! 
! Solving a non-linear 2-D shallow water system on a sphere. 
!     
!     du/dt - 2\Omega\sin\phi v = -g/a\cos\phi dh/d\lambda + D(u),
!     dv/dt + 2\Omega\sin\phi u = -g/a  dh/d\phi + D(v)
!     dH/dt + ( 1/\cos\phi d([H-hs] u)/d\lambda 
!               + 1/\cos\phi d([H-hs] v\cos\phi)/d\phi ) = D(h).
!
! where hs is hight of topography, D(a) is hyper viscosity. 
!
! The time integration is performed by using 3 step Adams-Bashforth 
! scheme. 
! 
program shallow_zd

!== モジュール引用宣言 ================================================

  use w_module 
  use gt4_history, only : HistoryCreate, HistoryPut, HistoryGet, HistoryAddVariable, HistoryClose
  use dc_trace,    only : SetDebug, BeginSub, EndSub, DbgMessage
  use dc_message,  only : MessageNotify
  implicit none

!== 宣言部 ============================================================

!----------------------------------------------------------------------
!  NAMELIST ファイル名
!----------------------------------------------------------------------

  character(len=100) :: NameListFile='shallow_zd.nml'

!----------------------------------------------------------------------
!  変数(格子点データ)
!----------------------------------------------------------------------

  real(8), allocatable :: xy_VelLon_n(:,:)      ! 速度経度成分 (t)
  real(8), allocatable :: xy_VelLat_n(:,:)      ! 速度緯度成分 (t)
  real(8), allocatable :: xy_Hsfc_n(:,:)        ! 表面変位 (t)
  real(8), allocatable :: xy_Vor_n(:,:)         ! 渦度 (t)
  real(8), allocatable :: xy_HDiv_n(:,:)        ! 発散 (t)

  real(8), allocatable :: xy_VelLon_b(:,:)      ! 速度経度成分 (t-Δt)
  real(8), allocatable :: xy_VelLat_b(:,:)      ! 速度緯度成分 (t-Δt)
  real(8), allocatable :: xy_Hsfc_b(:,:)        ! 表面変位 (t-Δt)
  real(8), allocatable :: xy_Vor_b(:,:)         ! 渦度 (t-Δt)
  real(8), allocatable :: xy_HDiv_b(:,:)        ! 発散 (t-Δt)

  real(8), allocatable :: xy_VelLon_bb(:,:)     ! 速度経度成分 (t-2Δt)
  real(8), allocatable :: xy_VelLat_bb(:,:)     ! 速度緯度成分 (t-2Δt)
  real(8), allocatable :: xy_Hsfc_bb(:,:)       ! 表面変位 (t-2Δt)
  real(8), allocatable :: xy_Vor_bb(:,:)        ! 渦度 (t-2Δt)
  real(8), allocatable :: xy_HDiv_bb(:,:)       ! 発散 (t-2Δt)

  real(8), allocatable :: xy_Coli(:,:)          ! 惑星渦度
  real(8), allocatable :: xy_Topo(:,:)          ! 表面地形

!----------------------------------------------------------------------
!  変数(スペクトルデータ)
!----------------------------------------------------------------------

  real(8), allocatable :: w_Hsfc_a(:)           ! 変位 (t+Δt)
  real(8), allocatable :: w_Vor_a(:)            ! 渦度 (t+Δt)
  real(8), allocatable :: w_HDiv_a(:)           ! 発散 (t+Δt)

  real(8), allocatable :: w_Hsfc_n(:)           ! 変位 (t)
  real(8), allocatable :: w_Vor_n(:)            ! 渦度 (t)
  real(8), allocatable :: w_HDiv_n(:)           ! 発散 (t)

  real(8), allocatable :: w_Hsfc_b(:)           ! 変位 (t-Δt)
  real(8), allocatable :: w_Vor_b(:)            ! 渦度 (t-Δt)
  real(8), allocatable :: w_HDiv_b(:)           ! 発散 (t-Δt)

  real(8), allocatable :: w_dHsfc_n(:)          ! 変位増分(t)
  real(8), allocatable :: w_dVor_n(:)           ! 渦度増分(t)
  real(8), allocatable :: w_dHDiv_n(:)          ! 発散増分(t)

  real(8), allocatable :: w_dHsfc_b(:)          ! 変位増分 (t-Δt)
  real(8), allocatable :: w_dVor_b(:)           ! 渦度増分 (t-Δt)
  real(8), allocatable :: w_dHDiv_b(:)          ! 発散増分 (t-Δt)

  real(8), allocatable :: w_dHsfc_bb(:)         ! 変位増分 (t-2Δt)
  real(8), allocatable :: w_dVor_bb(:)          ! 渦度増分 (t-2Δt)
  real(8), allocatable :: w_dHDiv_bb(:)         ! 発散増分 (t-2Δt)

!----------------------------------------------------------------------
!  固定パラメタ
!----------------------------------------------------------------------

  real(8), parameter  :: pi = 3.141592653589793D0 ! 円周率

!----------------------------------------------------------------------
!  NAMELIST 変数 (実験設定等)
!----------------------------------------------------------------------

  character(len=100) :: ExpTitle=''              ! 実験名
  character(len=100) :: ExpInst =''              ! 実行者所属
  character(len=100) :: ExpSrc  =''              ! 実行者名
  character(len=100) :: ExpModel=''              ! モデル設定名

  namelist /expset/  ExpTitle, ExpInst, ExpSrc, ExpModel

!----------------------------------------------------------------------
!  NAMELIST 変数 (時間単位)
!----------------------------------------------------------------------

  integer            :: UnitYear  = 360          ! 1 年の長さ(日)
  integer            :: UnitDay   = 24           ! 1 日の長さ(時間)
  integer            :: UnitHour  = 60           ! 1 時間の長さ(分)    
  integer            :: UnitMinit = 60           ! 1 分の長さ(秒)

  namelist /dateset/ UnitYear, UnitDay, UnitHour, UnitMinit

!----------------------------------------------------------------------
!  NAMELIST 変数 (時間積分パラメタ)
!----------------------------------------------------------------------

  real(8)            :: TimeIntYear 		 ! 積分時間(年)
  real(8)            :: TimeIntDay               ! 積分時間(日)
  real(8)            :: TimeIntHour              ! 積分時間(時)
  real(8)            :: TimeIntMinit      	 ! 積分時間(分)
  real(8)            :: TimeIntSec               ! 積分時間(秒)
  real(8)            :: OutputYear  		 ! 出力間隔(年)
  real(8)            :: OutputDay                ! 出力間隔(日)
  real(8)            :: OutputHour               ! 出力間隔(時)
  real(8)            :: OutputMinit              ! 出力間隔(分)
  real(8)            :: OutputSec                ! 出力間隔(秒)
  real(8)            :: DelTime                  ! 時間刻み(秒)
  character(len=100) :: IntScheme=''             ! 時間積分法

  namelist /timeset/ &
&           TimeIntYear, TimeIntDay, TimeIntHour, TimeIntMinit, TimeIntSec, &
&           OutputYear,  OutputDay,  OutputHour,  OutputMinit,  OutputSec, &
&           DelTime, IntScheme  

!----------------------------------------------------------------------
!  NAMELIST 変数 (入出力ファイル名)
!----------------------------------------------------------------------

  character(len=100) :: InputRstFile=''          ! 入力リスタートファイル名  
  character(len=100) :: OutputRstfile=''         ! 出力リスタートファイル名  
  character(len=100) :: OutputFile=''            ! 出力ファイル名  

  namelist /fileset/ InputRstFile, OutputRstfile, OutputFile

!----------------------------------------------------------------------
!  NAMELIST 変数(空間解像度設定)
!----------------------------------------------------------------------

  integer :: nm=21             ! 切断全波数
  integer :: im=64             ! 経度方向格子点数 ( > 3*nm + 1)
  integer :: jm=32             ! 緯度方向格子点数 ( > 3*nm/2 )

  namelist /gridset/ nm, im, jm 

!----------------------------------------------------------------------
!  NAMELIST 変数(物理パラメタ)
!----------------------------------------------------------------------

  real(8)            :: Rplanet  = 6.371D6       ! 惑星半径
  real(8)            :: Omega    = 7.292D-5      ! 回転角速度
  real(8)            :: Grav     = 9.8D0         ! 重力加速度
  real(8)            :: TfilCoef = 0.05          ! 時間フィルター係数
  integer            :: VisOrder = 2             ! 超粘性の次数
  real(8)            :: VisCoef  = 2.338D16      ! 超粘性係数

  namelist /paramset/ Rplanet, Omega, Grav, &
&                     TfilCoef, VisOrder, VisCoef

!----------------------------------------------------------------------
!  NAMELIST 変数 (デバッグ出力制御)
!----------------------------------------------------------------------

  logical            :: DebugOn = .false.        ! デバッグ出力制御

  namelist /debugset/ DebugOn

!----------------------------------------------------------------------
!  内部変数
!----------------------------------------------------------------------

  real(8)            :: delta                   ! 浅水/順圧切替え制御
  real(8)            :: alpha                   ! 陰解/陽解法切替え制御
  real(8)            :: GeoPotAvr               ! 平均ポテンシャル
  real(8)            :: HsfcAvr                 ! 平均ポテンシャル

  real(8)            :: Time                    ! 時間(秒)
  integer            :: NStep                   ! 時間ステップ数
  integer            :: NStepAll                ! 全時間ステップ数
  integer            :: NStepOut                ! 出力ステップ間隔

  real(8)            :: InitTime                ! 積分開始時刻
  real(8)            :: TimeIntAll              ! 全積分時間(秒)
  real(8)            :: OutputAll               ! 全出力間隔(秒)

  real(8)            :: MassAll                 ! 全質量
  real(8)            :: EnergyAll               ! 全エネルギー
  real(8)            :: AMomentAll              ! 全角運動量エネルギー

  character(len=20) :: DbgMessageFmt='*** DbgMESSAGE ***  '

!== 初期化 ============================================================

!----------------------------------------------------------------------
!  NAMELIST 読み込み 
!----------------------------------------------------------------------

  open(10,file=trim(NameListFile))
  read(10,nml=expset)  
  read(10,nml=dateset) 
  read(10,nml=timeset) 
  read(10,nml=fileset) 
  read(10,nml=gridset) 
  read(10,nml=paramset) 
  read(10,nml=debugset) 
  close(10)

!----------------------------------------------------------------------
!  デバッグ出力制御設定 
!  DebugOn が .true. なら SetDebug を呼ぶ
!----------------------------------------------------------------------

  if (DebugOn) then
    call SetDebug
  end if  

!----------------------------------------------------------------------
!  動的割り付け配列変数の初期化
!----------------------------------------------------------------------

  allocate(xy_VelLon_n(im,jm) , xy_VelLat_n(im,jm) , xy_Hsfc_n(im,jm) , &
&          xy_Vor_n(im,jm)    , xy_HDiv_n(im,jm)   ,                    & 
&          xy_VelLon_b(im,jm) , xy_VelLat_b(im,jm) , xy_Hsfc_b(im,jm) , &
&          xy_Vor_b(im,jm)    , xy_HDiv_b(im,jm)   ,                    &
&          xy_VelLon_bb(im,jm), xy_VelLat_bb(im,jm), xy_Hsfc_bb(im,jm), &
&          xy_Vor_bb(im,jm)   , xy_HDiv_bb(im,jm))

  allocate(xy_Coli(im,jm), xy_Topo(im,jm)) 
  
  allocate(w_Hsfc_a((nm+1)*(nm+1)), &
&          w_Vor_a((nm+1)*(nm+1)) , &
&          w_HDiv_a((nm+1)*(nm+1)), &
&          w_Hsfc_n((nm+1)*(nm+1)), &
&          w_Vor_n((nm+1)*(nm+1)) , &
&          w_HDiv_n((nm+1)*(nm+1)), &
&          w_Hsfc_b((nm+1)*(nm+1)), &
&          w_Vor_b((nm+1)*(nm+1)) , &
&          w_HDiv_b((nm+1)*(nm+1)))

  allocate(w_dHsfc_n((nm+1)*(nm+1)) , &
&          w_dVor_n((nm+1)*(nm+1))  , &
&          w_dHDiv_n((nm+1)*(nm+1)) , &
&          w_dHsfc_b((nm+1)*(nm+1)) , &
&          w_dVor_b((nm+1)*(nm+1))  , &
&          w_dHDiv_b((nm+1)*(nm+1)) , &
&          w_dHsfc_bb((nm+1)*(nm+1)), &
&          w_dVor_bb((nm+1)*(nm+1)) , &
&          w_dHDiv_bb((nm+1)*(nm+1))) 


!----------------------------------------------------------------------
!  spml ライブラリの初期化
!  座標変数 xy_Lon, xy_Lat x_Lon, y_Lat はここで値が確定する.
!----------------------------------------------------------------------

  call DbgMessage(fmt='call %c', c1='w_initial') 
  call w_initial(nm,im,jm)                       ! ISPACK初期化

!----------------------------------------------------------------------
!  リスタートファイル読み込み
!  Hsfc = Hsfc - HsfcAvr とする. 
!----------------------------------------------------------------------

  xy_Coli = 2 * Omega * sin(xy_Lat) 

  call InputRestart_gtool4(InputRstFile, InitTime, xy_Topo, &
&              xy_VelLon_n, xy_VelLat_n, xy_Hsfc_n, xy_Vor_n, xy_HDiv_n, &
&              xy_VelLon_b, xy_VelLat_b, xy_Hsfc_b, xy_Vor_b, xy_HDiv_b, &
&              xy_VelLon_bb, xy_VelLat_bb, xy_Hsfc_bb, xy_Vor_bb, xy_HDiv_bb  )
 
  HsfcAvr   = AvrLonLat_xy(xy_Hsfc_n)
  GeoPotAvr = Grav*HsfcAvr 

  xy_Hsfc_n  = xy_Hsfc_n  - HsfcAvr
  xy_Hsfc_b  = xy_Hsfc_b  - HsfcAvr
  xy_Hsfc_bb = xy_Hsfc_bb - HsfcAvr

  w_Hsfc_n = w_xy(xy_Hsfc_n) ; w_Hsfc_b = w_xy(xy_Hsfc_b) ; w_Hsfc_a = w_Hsfc_n
  w_Vor_n  = w_xy(xy_Vor_n)  ; w_Vor_b  = w_xy(xy_Vor_b)  ; w_Vor_a  = w_Vor_n
  w_HDiv_n = w_xy(xy_HDiv_n) ; w_HDiv_b = w_xy(xy_HDiv_b) ; w_HDiv_a = w_HDiv_n

!----------------------------------------------------------------------
!  出力ファイル初期化
!----------------------------------------------------------------------

  call Output_gtool4_init                        

!----------------------------------------------------------------------
!  浅水/順圧切替え
!  変数 ExpModel で切替え. デフォルトは浅水モデル. 
!----------------------------------------------------------------------

  select case (ExpModel)
  case ('baro') 
    delta = 0.0D0
    call MessageNotify('Message',     &
&                      'shallow_zd',  &
&                      'delta=%f; %c', &
&                       d=(/delta/),    &
&                       c1='barotoropic vorticity equation') 
  case default
    delta = 1.0D0
    call MessageNotify('Message',     &
&                      'shallow_zd',  &
&                      'delta=%f; %c', &
&                       d=(/delta/),    &
&                       c1='shallow water equation') 
  end select

!----------------------------------------------------------------------
!  時間積分スキームの切替えパラメータの設定
!  変数 IntScheme で切替え. デフォルトは semi-implicit
!----------------------------------------------------------------------

  select case (IntScheme)
  case ('explicit') 
    alpha = 0.0D0
    call MessageNotify('Message',     &
&                      'shallow_zd',  &
&                      'alpha=%f; %c', &
&                       d=(/alpha/),    &
&                       c1='explicit time integration' )
  case default
    alpha = 1.0D0
    call MessageNotify('Message',     &
&                      'shallow_zd',  &
&                      'alpha=%f; %c', &
&                       d=(/alpha/),    &
&                       c1='semi-implicit time integration') 
  end select

!----------------------------------------------------------------------
!  時間積分設定
!  NAMELIST で入力した全積分時間と出力時間間隔を年月日分秒から秒へ変換し, 
!  DelTime を用いて全計算ステップ数 NstepAll と出力ステップ数 NstepOut
! を計算する.
!----------------------------------------------------------------------

  Time = InitTime

  if (DelTime == 0.0) then
    call MessageNotify('Error',      &
&                      'shallow_zd', &
&                      'DelTime is 0.')
  end if 
 
  TimeIntAll = (TimeIntYear * UnitYear * UnitDay * UnitHour * UnitMinit &
&             + TimeIntDay             * UnitDay * UnitHour * UnitMinit &
&             + TimeIntHour                      * UnitHour * UnitMinit &
&             + TimeIntMinit                                * UnitMinit &
&             + TimeIntSec) 

  if (mod(TimeIntAll, DelTime) /= 0) then
    call MessageNotify('Error',      &
&                      'shallow_zd', &
&                      'mod(TimeIntAll, DelTime) is not 0.')
  else
    NstepAll = INT(TimeIntAll / DelTime)
  end if 

  OutputAll = (OutputYear * UnitYear * UnitDay * UnitHour * UnitMinit &
&            + OutputDay             * UnitDay * UnitHour * UnitMinit &
&            + OutputHour                      * UnitHour * UnitMinit &
&            + OutputMinit                                * UnitMinit &
&            + OutputSec) 

  if (mod(OutputAll, DelTime) /= 0) then
    call MessageNotify('Error',      &
&                      'shallow_zd', &
&                      'mod(OutputAll, DelTime) is not 0.')
  else
    NstepOut = INT(OutputAll / DelTime)
  end if 

  call DbgMessage(fmt='%c TimeIntAll=%f', d=(/TimeIntALL/), c1=DbgMessageFmt)
  call DbgMessage(fmt='%c OutPutAll=%f',  d=(/OutPutALL/), c1=DbgMessageFmt)
  call DbgMessage(fmt='%c NstepAll=%d', i=(/NstepALL/), c1=DbgMessageFmt)
  call DbgMessage(fmt='%c NstepOut=%d', i=(/NstepOut/), c1=DbgMessageFmt)


!== 時間積分 ==========================================================

  call DbgMessage(fmt='%c %c', &
&                 c1=DbgMessageFmt, &
&                 c2='Time integration starts.') 

!----------------------------------------------------------------------
!  時刻 t-2Δt における増分の計算
!----------------------------------------------------------------------
 
  w_dVor_bb = w_dVor_xy(xy_Coli + xy_Vor_bb,               &
&                       xy_VelLon_bb, xy_VelLat_bb)       
     
  w_dHDiv_bb = w_dHDiv_xy(xy_Coli + xy_Vor_bb,             &
                          xy_VelLon_bb, xy_VelLat_bb, xy_Hsfc_bb)
&
  w_dHsfc_bb = w_dHsfc_xy(xy_VelLon_bb, xy_VelLat_bb,      &
&                         xy_Hsfc_bb - xy_Topo, xy_HDiv_bb)

!----------------------------------------------------------------------
!  時刻 t-Δt における増分の計算
!----------------------------------------------------------------------
 
  w_dVor_b = w_dVor_xy(xy_Coli + xy_Vor_b,              &
&                      xy_VelLon_b, xy_VelLat_b)

  w_dHDiv_b = w_dHDiv_xy(xy_Coli + xy_Vor_b,             &
                         xy_VelLon_b, xy_VelLat_b, xy_Hsfc_b)
&
  w_dHsfc_b = w_dHsfc_xy(xy_VelLon_b, xy_VelLat_b, &
&                        xy_Hsfc_b - xy_Topo, xy_HDiv_b)

!----------------------------------------------------------------------
!  最初の 2 ステップ分の時間積分開始
!----------------------------------------------------------------------

  do Nstep = 1, 2

    Time = Time + DelTime
    call DbgMessage(fmt='%c Time=%f', d=(/Time/), c1=DbgMessageFmt) 

!----------------------------------------------------------------------
!  時刻 t における増分の計算
!----------------------------------------------------------------------
 
    w_dVor_n = w_dVor_xy(xy_Coli + xy_Vor_n,               & 
&                        xy_VelLon_n, xy_VelLat_n)

    w_dHDiv_n = w_dHDiv_xy(xy_Coli + xy_Vor_n,             &
&                           xy_VelLon_n, xy_VelLat_n, xy_Hsfc_n)

    w_dHsfc_n = w_dHsfc_xy(xy_VelLon_n, xy_VelLat_n,       &
&                          xy_Hsfc_n - xy_Topo, xy_HDiv_n)

!----------------------------------------------------------------------
!  InitTime = 0.0 の場合 2 段 Adams-Bashforth スキームで計算
!  (1 ステップ目はオイラー法), その他の場合は 3 段 Adams-Bashforth 
!  スキームで計算
!----------------------------------------------------------------------

    if (InitTime == 0.0) then

      call DbgMessage(fmt='%c InitTime=%f : %c', &
&                     d=(/InitTime/),            &
&                     c1=DbgMessageFmt,          &
&                     c2='2 step Adams-Bashforth scheme is adapted.') 

      w_Vor_a  = w_Vor_n  &
&                + w_AdamsBashforth2_w(w_dVor_n, w_dVor_b) &
&                + DelTime*w_NumVis_w(w_Vor_n)

      w_HDiv_a = w_HDiv_n & 
&                + w_AdamsBashforth2_w(w_dHDiv_n, w_dHDiv_b) &
&                + DelTime*w_NumVis_w(w_HDiv_n)

      w_Hsfc_a = w_Hsfc_n &
&                + w_AdamsBashforth2_w(w_dHsfc_n, w_dHsfc_b) & 
&                + DelTime*w_NumVisScaler_w(w_Hsfc_n)

    else

      w_Vor_a  = w_Vor_n  &
&                + w_AdamsBashforth3_w(w_dVor_n, w_dVor_b, w_dVor_bb ) &
&                + DelTime*w_NumVis_w(w_Vor_n)

      w_HDiv_a = w_HDiv_n & 
&                + w_AdamsBashforth3_w(w_dHDiv_n, w_dHDiv_b, w_dHDiv_bb) &
&                + DelTime*w_NumVis_w(w_HDiv_n)

      w_Hsfc_a = w_Hsfc_n &
&                + w_AdamsBashforth3_w(w_dHsfc_n, w_dHsfc_b, w_dHsfc_bb ) &
&                + DelTime*w_NumVisScaler_w(w_Hsfc_n)

    end if
 
!----------------------------------------------------------------------
!  変数と増分を 1 ステップ更新する.
!----------------------------------------------------------------------

    call RoopVariable_w(w_Vor_a , w_Vor_n , w_Vor_b )
    call RoopVariable_w(w_HDiv_a, w_HDiv_n, w_HDiv_b)
    call RoopVariable_w(w_Hsfc_a, w_Hsfc_n, w_Hsfc_b)

    call RoopVariable_w(w_dVor_n , w_dVor_b , w_dVor_bb )
    call RoopVariable_w(w_dHDiv_n, w_dHDiv_b, w_dHDiv_bb)
    call RoopVariable_w(w_dHsfc_n, w_dHsfc_b, w_dHsfc_bb)

    call RoopVariable_xy(xy_Vor_n   , xy_Vor_b   , xy_Vor_bb   )
    call RoopVariable_xy(xy_HDiv_n  , xy_HDiv_b  , xy_HDiv_bb  )
    call RoopVariable_xy(xy_Hsfc_n  , xy_Hsfc_b  , xy_Hsfc_bb  )
    call RoopVariable_xy(xy_VelLon_n, xy_VelLon_b, xy_VelLon_bb)
    call RoopVariable_xy(xy_VelLat_n, xy_VelLat_b, xy_VelLat_bb)

!----------------------------------------------------------------------
!  格子点データの計算
!----------------------------------------------------------------------

    call uv(w_Vor_n, w_HDiv_n, xy_VelLon_n, xy_VelLat_n) 

    xy_Vor_n  = xy_w(w_Vor_n)
    xy_HDiv_n = xy_w(w_HDiv_n)
    xy_Hsfc_n = xy_w(w_Hsfc_n)

!----------------------------------------------------------------------
!  最初の 2 ステップ分の時間積分終了
!----------------------------------------------------------------------

  end do  

!----------------------------------------------------------------------
!  3 ステップ以降の時間積分開始
!----------------------------------------------------------------------

  do Nstep = 3, NstepALL

    Time = Time + DelTime
    call DbgMessage(fmt='%c Time=%f', d=(/Time/), c1=DbgMessageFmt) 

!----------------------------------------------------------------------
!  時刻 t における増分の計算
!----------------------------------------------------------------------

    w_dVor_n = w_dVor_xy(xy_Coli + xy_Vor_n,              &
&                        xy_VelLon_n, xy_VelLat_n)

    w_dHDiv_n = w_dHDiv_xy(xy_Coli + xy_Vor_n,             &
&                          xy_VelLon_n, xy_VelLat_n, xy_Hsfc_n)

    w_dHsfc_n = w_dHsfc_xy(xy_VelLon_n, xy_VelLat_n, &
&                          xy_Hsfc_n - xy_Topo, xy_HDiv_n)

!----------------------------------------------------------------------
!  1 ステップ計算
!----------------------------------------------------------------------

    w_Vor_a  = w_Vor_n  &
&              + w_AdamsBashforth3_w(w_dVor_n, w_dVor_b, w_dVor_bb ) &
&                + DelTime*w_NumVis_w(w_Vor_n)

    w_HDiv_a = w_HDiv_n & 
&              + w_AdamsBashforth3_w(w_dHDiv_n, w_dHDiv_b, w_dHDiv_bb) &
&                + DelTime*w_NumVis_w(w_HDiv_n)

    w_Hsfc_a = w_Hsfc_n &
&              + w_AdamsBashforth3_w(w_dHsfc_n, w_dHsfc_b, w_dHsfc_bb) &
&                + DelTime*w_NumVisScaler_w(w_Hsfc_n)

!----------------------------------------------------------------------
!  変数と増分を 1 ステップ更新する.
!----------------------------------------------------------------------

    call RoopVariable_w(w_Vor_a , w_Vor_n , w_Vor_b )
    call RoopVariable_w(w_HDiv_a, w_HDiv_n, w_HDiv_b)
    call RoopVariable_w(w_Hsfc_a, w_Hsfc_n, w_Hsfc_b)

    call RoopVariable_w(w_dVor_n , w_dVor_b , w_dVor_bb )
    call RoopVariable_w(w_dHDiv_n, w_dHDiv_b, w_dHDiv_bb)
    call RoopVariable_w(w_dHsfc_n, w_dHsfc_b, w_dHsfc_bb)

    call RoopVariable_xy(xy_Vor_n   , xy_Vor_b   , xy_Vor_bb   )
    call RoopVariable_xy(xy_HDiv_n  , xy_HDiv_b  , xy_HDiv_bb  )
    call RoopVariable_xy(xy_Hsfc_n  , xy_Hsfc_b  , xy_Hsfc_bb  )
    call RoopVariable_xy(xy_VelLon_n, xy_VelLon_b, xy_VelLon_bb)
    call RoopVariable_xy(xy_VelLat_n, xy_VelLat_b, xy_VelLat_bb)

!----------------------------------------------------------------------
!  格子点データの計算
!----------------------------------------------------------------------

    call uv(w_Vor_n, w_HDiv_n, xy_VelLon_n, xy_VelLat_n) 

    xy_Vor_n  = xy_w(w_Vor_n)
    xy_HDiv_n = xy_w(w_HDiv_n)
    xy_Hsfc_n = xy_w(w_Hsfc_n)

!----------------------------------------------------------------------
!  出力
!  Nstep が NstepOut の整数倍になったら出力を行う.
!----------------------------------------------------------------------

    if(mod(Nstep, NStepOut) == 0) then             

      MassAll   = &
&       IntLonLat_xy(xy_Hsfc_n + HsfcAvr - xy_Topo)

      EnergyAll =                                           &
&       IntLonLat_xy((xy_Hsfc_n + HsfcAvr - xy_Topo) * (    &
&                     xy_VelLon_n**2 + xy_VelLat_n**2 +     &
&                     xy_Hsfc_n + HsfcAvr - xy_Topo)*0.5D0) 

      AMomentAll = &
&       IntLonLat_xy((xy_Hsfc_n + HsfcAvr - xy_Topo)*cos(xy_Lat)*Rplanet &
&                     *(xy_VelLon_n + Omega*cos(xy_Lat)*Rplanet))

      call MessageNotify('Message',      &
&                        'shallow_zd',   &
&                        'Time=%f', d=(/Time/))
      call MessageNotify('Message',      &
&                        'shallow_zd',   &
&                        'Total Mass=%f', d=(/MassAll/))
      call MessageNotify('Message',      &
&                        'shallow_zd',   &
&                        'Total Energy=%f', d=(/EnergyAll/))
      call MessageNotify('Message',      &
&                        'shallow_zd',   &
&                        'Total Angular Momentum=%f', d=(/AMomentAll/))
      call Output_gtool4
    endif

  enddo

  call DbgMessage(fmt='%c %c',      &
&                 c1=DbgMessageFmt, &
&                 c2='Time integration end.') 


!== 後処理 ==========================================================

!----------------------------------------------------------------------
!  ヒストリー出力後処理
!----------------------------------------------------------------------

  call Output_gtool4_close

!----------------------------------------------------------------------
!  リスタートファイル出力
!----------------------------------------------------------------------

  call OutputRestart_gtool4_init              ! リスタートファイル初期化
  call OutputRestart_gtool4                   ! リスタートファイル出力
  call Output_gtool4_close                    ! リスタートファイル出力終了

!----------------------------------------------------------------------
!  動的割り付け配列変数の解放
!----------------------------------------------------------------------

  deallocate(xy_VelLon_n , xy_VelLat_n , xy_Hsfc_n,  &
&            xy_Vor_n    , xy_HDiv_n   ,             & 
&            xy_VelLon_b , xy_VelLat_b , xy_Hsfc_b,  &
&            xy_Vor_b    , xy_HDiv_b   ,             &
&            xy_VelLon_bb, xy_VelLat_bb, xy_Hsfc_bb, &
&            xy_Vor_bb   , xy_HDiv_bb  ,             &
&            xy_Coli     , xy_Topo      )

  deallocate(w_Hsfc_a, w_Vor_a, w_HDiv_a,  &
&            w_Hsfc_n, w_Vor_n, w_HDiv_n,  &
&            w_Hsfc_b, w_Vor_b, w_HDiv_b )

  deallocate(w_dHsfc_n , w_dVor_n , w_dHDiv_n, &
&            w_dHsfc_b , w_dVor_b , w_dHDiv_b, &
&            w_dHsfc_bb, w_dVor_bb, w_dHDiv_bb)

  stop

contains

!== 内部副プログラム =================================================

!----------------------------------------------------------------------
!  w_dVor_xy
!  渦度方程式の増分 w_dVor を計算する
!----------------------------------------------------------------------

  function w_dVor_xy(xy_VorAll, xy_VelLon, xy_VelLat)

    real(8)            :: w_dVor_xy((nm+1)*(nm+1))
    real(8),intent(in) :: xy_VorAll(im,jm)
    real(8),intent(in) :: xy_VelLon(im,jm)
    real(8),intent(in) :: xy_VelLat(im,jm)

    call BeginSub('w_dVor_xy')

    w_dVor_xy = DelTime * (                       &
&     - w_Div_xy_xy(xy_VorAll*xy_VelLon,          &
&                   xy_VorAll*xy_VelLat)/Rplanet)  

    call EndSub('w_dVor_xy')

  end function w_dVor_xy

!----------------------------------------------------------------------
!  w_dHDiv_xy
!  発散方程式の増分 w_dHDiv を計算する
!  (線形項を除く)
!----------------------------------------------------------------------

  function w_dHDiv_xy(xy_VorAll, xy_VelLon, xy_VelLat, xy_Hsfc)

    real(8)            :: w_dHDiv_xy((nm+1)*(nm+1))
    real(8),intent(in) :: xy_VorAll(im,jm)
    real(8),intent(in) :: xy_VelLon(im,jm)
    real(8),intent(in) :: xy_VelLat(im,jm)
    real(8),intent(in) :: xy_Hsfc(im,jm)

    call BeginSub('w_dHDiv_xy')

    w_dHDiv_xy = DelTime*delta*(                                      &
&     w_Div_xy_xy(xy_VorAll*xy_VelLat, - xy_VorAll*xy_VelLon)/Rplanet &
&     - w_Lapla_w(Grav*w_xy(xy_Hsfc)                                  &
&                 + 0.5D0*w_xy(xy_VelLon**2 + xy_VelLat**2))/Rplanet**2 )

    call EndSub('w_dHDiv_xy')
     
  end function w_dHDiv_xy

!----------------------------------------------------------------------
!  w_dHsfc_xy
!  質量保存の式の増分 w_dHsfc を計算する
!  (線形項を除く)
!----------------------------------------------------------------------

  function w_dHsfc_xy(xy_VelLon, xy_VelLat, xy_Hsfc, xy_HDiv)
  
    real(8)            :: w_dHsfc_xy((nm+1)*(nm+1))
    real(8),intent(in) :: xy_VelLon(im,jm)
    real(8),intent(in) :: xy_VelLat(im,jm)
    real(8),intent(in) :: xy_Hsfc(im,jm)
    real(8),intent(in) :: xy_HDiv(im,jm)

    w_dHsfc_xy = DelTime*(                                         &
&     - w_Div_xy_xy(xy_Hsfc*xy_VelLon, xy_Hsfc*xy_VelLat)/Rplanet  &
&     - HsfcAvr*w_xy(xy_HDiv))

  end function w_dHsfc_xy

!----------------------------------------------------------------------
!  w_NumVis_w
!  数値粘性項の計算
!----------------------------------------------------------------------
  function w_NumVis_w(w_data)

    real(8)             :: w_NumVis_w((nm+1)*(nm+1))
    real(8), intent(in) :: w_data((nm+1)*(nm+1))

    call BeginSub('w_NumVis_w')
  
    w_NumVis_w = w_NumVisScaler_w(w_data)                       &
&     - VisCoef*(- (2.0D0/Rplanet**2)**(VisOrder/2.0D0))*w_data

    call EndSub('w_NumVis_w')

  end function w_NumVis_w

!----------------------------------------------------------------------
!  w_NumVisScaler_w
!  スカラーに対する数値粘性項の計算
!----------------------------------------------------------------------

  function w_NumVisScaler_w(w_data)

    real(8)             :: w_NumVisScaler_w((nm+1)*(nm+1))
    real(8), intent(in) :: w_data((nm+1)*(nm+1))

    call BeginSub('w_NumVisScaler_w')
  
    w_NumVisScaler_w =   &
&     - VisCoef*( (- rn(:,1)/Rplanet**2)**(VisOrder/2.0D0) )*w_data

    call EndSub('w_NumVisScaler_w')

  end function w_NumVisScaler_w

!----------------------------------------------------------------------
!  uv
!  流線関数と速度ポテンシャル(スペクトルデータ)から, 水平風(格子データ)
!  を計算する.
!----------------------------------------------------------------------

  subroutine uv(w_Vor, w_HDiv, xy_VelLon, xy_VelLat)

    real(8),intent(in)  :: w_Vor((nm+1)*(nm+1))
    real(8),intent(in)  :: w_HDiv((nm+1)*(nm+1))
    real(8),intent(out) :: xy_VelLon(im,jm)
    real(8),intent(out) :: xy_VelLat(im,jm)
    real(8)             :: w_Stream((nm+1)*(nm+1))  ! 流線関数
    real(8)             :: w_VelPot((nm+1)*(nm+1))  ! 速度ポテンシャル

    call BeginSub('uv')

    w_Stream = w_LaplaInv_w(w_Vor)*Rplanet**2
    w_VelPot = w_LaplaInv_w(w_HDiv)*Rplanet**2 

    xy_VelLon = (delta*xy_GradLon_w(w_VelPot) - xy_GradLat_w(w_Stream))/Rplanet
    xy_VelLat = (xy_GradLon_w(w_Stream) + delta*xy_GradLat_w(w_VelPot))/Rplanet

    call EndSub('uv')

  end subroutine uv

!----------------------------------------------------------------------
!  2 段 Adams-Bashforth 法による増分計算
!----------------------------------------------------------------------

  function w_AdamsBashforth2_w(w_data_n, w_data_b)

    real(8)            :: w_AdamsBashforth2_w((nm+1)*(nm+1))
    real(8),intent(in) :: w_data_n((nm+1)*(nm+1))
    real(8),intent(in) :: w_data_b((nm+1)*(nm+1))

    w_AdamsBashforth2_w = (1.5 * w_data_n - 0.5 * w_data_b)

  end function w_AdamsBashforth2_w

!----------------------------------------------------------------------
!  3 段 Adams-Bashforth 法による増分計算
!----------------------------------------------------------------------

  function w_AdamsBashforth3_w(w_data_n, w_data_b, w_data_bb)

    real(8)            :: w_AdamsBashforth3_w((nm+1)*(nm+1))
    real(8),intent(in) :: w_data_n((nm+1)*(nm+1))
    real(8),intent(in) :: w_data_b((nm+1)*(nm+1))
    real(8),intent(in) :: w_data_bb((nm+1)*(nm+1))

    w_AdamsBashforth3_w = &
                 (23.0*w_data_n - 16.0*w_data_b + 5.0*w_data_bb) / 12.0

  end function w_AdamsBashforth3_w

!----------------------------------------------------------------------
!  変数を 1 ステップ更新する(スペクトルデータ)
!----------------------------------------------------------------------
  subroutine RoopVariable_w(w_data_a, w_data_n, w_data_b)

    real(8),intent(in)    :: w_data_a((nm+1)*(nm+1))
    real(8),intent(inout) :: w_data_n((nm+1)*(nm+1))
    real(8),intent(inout) :: w_data_b((nm+1)*(nm+1))

    call BeginSub('RoopVariable_w')

    w_data_b = w_data_n
    w_data_n = w_data_a 

    call EndSub('RoopVariable_w')

  end subroutine RoopVariable_w

!----------------------------------------------------------------------
!  変数を 1 ステップ更新する(格子点データ)
!----------------------------------------------------------------------
  subroutine RoopVariable_xy(xy_data_n, xy_data_b, xy_data_bb)

    real(8),intent(in)    :: xy_data_n(im,jm)
    real(8),intent(inout) :: xy_data_b(im,jm)
    real(8),intent(inout) :: xy_data_bb(im,jm)

    call BeginSub('RoopVariable_xy')

    xy_data_bb = xy_data_b
    xy_data_b  = xy_data_n 

    call EndSub('RoopVariable_xy')

  end subroutine RoopVariable_xy

!----------------------------------------------------------------------
!  InputRestart_gtool4
!  リスタートファイルの読み込み
!----------------------------------------------------------------------

  subroutine InputRestart_gtool4(filename, Time, xy_Topo,  &
&             xy_VelLon_n, xy_VelLat_n, xy_Hsfc_n, xy_Vor_n, xy_HDiv_n, &
&             xy_VelLon_b, xy_VelLat_b, xy_Hsfc_b, xy_Vor_b, xy_HDiv_b, & 
&             xy_VelLon_bb, xy_VelLat_bb, xy_Hsfc_bb, xy_Vor_bb, xy_HDiv_bb)

    character(len=*), intent(in) :: filename       ! 入力ファイル名
    real(8),intent(out)          :: Time           ! 時刻
    real(8),intent(out)          :: xy_Topo(im,jm) ! 地形

    real(8),intent(out)          :: xy_VelLon_n(im,jm) ! 速度経度成分
    real(8),intent(out)          :: xy_VelLat_n(im,jm) ! 速度緯度成分
    real(8),intent(out)          :: xy_Hsfc_n(im,jm) ! 変位
    real(8),intent(out)          :: xy_Vor_n(im,jm)  ! 渦度
    real(8),intent(out)          :: xy_HDiv_n(im,jm) ! 発散

    real(8),intent(out)          :: xy_VelLon_b(im,jm) ! 速度経度成分
    real(8),intent(out)          :: xy_VelLat_b(im,jm) ! 速度緯度成分
    real(8),intent(out)          :: xy_Hsfc_b(im,jm) ! 変位
    real(8),intent(out)          :: xy_Vor_b(im,jm)  ! 渦度
    real(8),intent(out)          :: xy_HDiv_b(im,jm) ! 発散

    real(8),intent(out)          :: xy_VelLon_bb(im,jm) ! 速度経度成分
    real(8),intent(out)          :: xy_VelLat_bb(im,jm) ! 速度緯度成分
    real(8),intent(out)          :: xy_Hsfc_bb(im,jm) ! 変位
    real(8),intent(out)          :: xy_Vor_bb(im,jm)  ! 渦度
    real(8),intent(out)          :: xy_HDiv_bb(im,jm) ! 発散

    call HistoryGet(filename, 't', Time)
    call HistoryGet(filename, 'topo', xy_Topo)

    call HistoryGet(filename, 'u_n'  , xy_VelLon_n, time=Time)
    call HistoryGet(filename, 'v_n'  , xy_VelLat_n, time=Time)
    call HistoryGet(filename, 'h_n'  , xy_Hsfc_n  , time=Time)
    call HistoryGet(filename, 'vor_n', xy_Vor_n   , time=Time)
    call HistoryGet(filename, 'div_n', xy_HDiv_n  , time=Time)

    call HistoryGet(filename, 'u_b'  , xy_VelLon_b, time=Time)
    call HistoryGet(filename, 'v_b'  , xy_VelLat_b, time=Time)
    call HistoryGet(filename, 'h_b'  , xy_Hsfc_b  , time=Time)
    call HistoryGet(filename, 'vor_b', xy_Vor_b   , time=Time)
    call HistoryGet(filename, 'div_b', xy_HDiv_b  , time=Time)

    call HistoryGet(filename, 'u_bb'  , xy_VelLon_bb, time=Time)
    call HistoryGet(filename, 'v_bb'  , xy_VelLat_bb, time=Time)
    call HistoryGet(filename, 'h_bb'  , xy_Hsfc_bb  , time=Time)
    call HistoryGet(filename, 'vor_bb', xy_Vor_bb   , time=Time)
    call HistoryGet(filename, 'div_bb', xy_HDiv_bb  , time=Time)

  end subroutine InputRestart_gtool4

!----------------------------------------------------------------------
!  Output_gtool4_init
!  出力ヒストリーファイルの初期化
!----------------------------------------------------------------------

  subroutine Output_gtool4_init
    call HistoryCreate( &                               ! ヒストリー作成
&          file=trim(OutputFile),                                & 
&          title=trim(ExpTitle),                                  &
&          source=trim(ExpSrc),                               &
&          institution=trim(ExpInst),                             &
&          dims=(/'lon','lat','t  '/), dimsizes=(/im,jm,0/),      &
&          longnames=(/'longitude','latitude ','time     '/),     &
&          units=(/'deg.','deg.','sec.'/),                        &
&          origin=real(InitTime), interval=real(NStepOut*DelTime) )

    call HistoryPut('lon',x_Lon*180/pi)              ! 変数出力
    call HistoryPut('lat',y_Lat*180/pi)              ! 変数出力

    call HistoryAddVariable( &                          ! 変数定義
&           varname='topo', dims=(/'lon','lat'/), & 
&           longname='Topography', units='m', xtype='double')

    call HistoryPut('topo',xy_Topo)

    call HistoryAddVariable( &                          ! 変数定義
&          varname='h', dims=(/'lon','lat','t  '/), & 
&          longname='surface displacement ', units='m', xtype='double')
    call HistoryAddVariable( &                          ! 変数定義
&          varname='u', dims=(/'lon','lat','t  '/), & 
&          longname='velocity(longitude) ', units='m/s', xtype='double')
    call HistoryAddVariable( &                          ! 変数定義
&          varname='v', dims=(/'lon','lat','t  '/), & 
&          longname='velocity(latitude) ', units='m/s', xtype='double')
    call HistoryAddVariable( &                          ! 変数定義
&          varname='vor', dims=(/'lon','lat','t  '/), & 
&          longname='vorticity', units='1/s', xtype='double')
    call HistoryAddVariable( &                          ! 変数定義
&          varname='div', dims=(/'lon','lat','t  '/), & 
&          longname='divergence', units='1/s', xtype='double')

  end subroutine Output_gtool4_init

!----------------------------------------------------------------------
!  Output_gtool4
!  ヒストリーファイルへの出力
!  Hsfc = Hsfc + HsfcAvr とする. 
!----------------------------------------------------------------------

  subroutine Output_gtool4
    call HistoryPut('t',Time)
    call HistoryPut('u'  ,xy_VelLon_n)
    call HistoryPut('v'  ,xy_VelLat_n)
    call HistoryPut('h'  ,xy_Hsfc_n + HsfcAvr)
    call HistoryPut('vor',xy_Vor_n)
    call HistoryPut('div',xy_HDiv_n)
  end subroutine Output_gtool4

!----------------------------------------------------------------------
!  Output_gtool4_close
!  ファイル入出力後処理
!----------------------------------------------------------------------

  subroutine Output_gtool4_close
    call HistoryClose
  end subroutine Output_gtool4_close

!----------------------------------------------------------------------
!  OutputRestart_gtool4_init
!  出力ファイルの初期化を行う
!  出力変数は xy_VelLon, xy_VelLat, xy_Hsfc, xy_Vor, xy_HDiv.
!  リープフロッグスキームを用いているため, 2 ステップ分の変数定義
!  を行う. 
!----------------------------------------------------------------------

  subroutine OutputRestart_gtool4_init

    call HistoryCreate( &                               ! ヒストリー作成
&          file=trim(OutputRstFile),                              & 
&          title=trim(ExpTitle),                                  &
&          source=trim(ExpSrc),                                   &
&          institution=trim(ExpInst),                             &
&          dims=(/'lon','lat','t  '/), dimsizes=(/im,jm,0/),      &
&          longnames=(/'longitude','latitude ','time     '/),     &
&          units=(/'deg.','deg.','sec.'/),                        &
&          origin=real(Time), interval=real(0.0) )

    call HistoryPut('lon',x_Lon*180/pi)                 ! 変数出力
    call HistoryPut('lat',y_Lat*180/pi)                 ! 変数出力

    call HistoryAddVariable( &                          ! 変数定義
&           varname='topo', dims=(/'lon','lat'/), & 
&           longname='Topography', units='m', xtype='double')

    call HistoryAddVariable( &                          ! 変数定義
&          varname='h_n', dims=(/'lon','lat','t  '/), & 
&          longname='surface displacement ', units='m', xtype='double')
    call HistoryAddVariable( &                          ! 変数定義
&          varname='u_n', dims=(/'lon','lat','t  '/), & 
&          longname='velocity(longitude) ', units='m/s', xtype='double')
    call HistoryAddVariable( &                          ! 変数定義
&          varname='v_n', dims=(/'lon','lat','t  '/), & 
&          longname='velocity(latitude) ', units='m/s', xtype='double')
    call HistoryAddVariable( &                          ! 変数定義
&          varname='vor_n', dims=(/'lon','lat','t  '/), & 
&          longname='vorticity', units='1/s', xtype='double')
    call HistoryAddVariable( &                          ! 変数定義
&          varname='div_n', dims=(/'lon','lat','t  '/), & 
&          longname='divergence', units='1/s', xtype='double')

    call HistoryAddVariable( &                          ! 変数定義
&          varname='h_b', dims=(/'lon','lat','t  '/), & 
&          longname='surface displacement ', units='m', xtype='double')
    call HistoryAddVariable( &                          ! 変数定義
&          varname='u_b', dims=(/'lon','lat','t  '/), & 
&          longname='velocity(longitude) ', units='m/s', xtype='double')
    call HistoryAddVariable( &                          ! 変数定義
&          varname='v_b', dims=(/'lon','lat','t  '/), & 
&          longname='velocity(latitude) ', units='m/s', xtype='double')
    call HistoryAddVariable( &                          ! 変数定義
&          varname='vor_b', dims=(/'lon','lat','t  '/), & 
&          longname='vorticity', units='1/s', xtype='double')
    call HistoryAddVariable( &                          ! 変数定義
&          varname='div_b', dims=(/'lon','lat','t  '/), & 
&          longname='divergence', units='1/s', xtype='double')

    call HistoryAddVariable( &                          ! 変数定義
&          varname='h_bb', dims=(/'lon','lat','t  '/), & 
&          longname='surface displacement ', units='m', xtype='double')
    call HistoryAddVariable( &                          ! 変数定義
&          varname='u_bb', dims=(/'lon','lat','t  '/), & 
&          longname='velocity(longitude) ', units='m/s', xtype='double')
    call HistoryAddVariable( &                          ! 変数定義
&          varname='v_bb', dims=(/'lon','lat','t  '/), & 
&          longname='velocity(latitude) ', units='m/s', xtype='double')
    call HistoryAddVariable( &                          ! 変数定義
&          varname='vor_bb', dims=(/'lon','lat','t  '/), & 
&          longname='vorticity', units='1/s', xtype='double')
    call HistoryAddVariable( &                          ! 変数定義
&          varname='div_bb', dims=(/'lon','lat','t  '/), & 
&          longname='divergence', units='1/s', xtype='double')

  end subroutine OutputRestart_gtool4_init

!----------------------------------------------------------------------
!  output_gtool4
!  初期値の出力を行う. 3 段の Adams-Bashforth スキームを用いているため, 
!  3 ステップ分のデータを出力する.
!  Hsfc = Hsfc + HsfcAvr とする. 
!----------------------------------------------------------------------

  subroutine OutputRestart_gtool4
    call HistoryPut('t',Time)
    call HistoryPut('topo',xy_Topo)

    call HistoryPut('u_n'  ,xy_VelLon_n)
    call HistoryPut('v_n'  ,xy_VelLat_n)
    call HistoryPut('h_n'  ,xy_Hsfc_n + HsfcAvr)
    call HistoryPut('vor_n',xy_Vor_n)
    call HistoryPut('div_n',xy_HDiv_n)

    call HistoryPut('u_b'  ,xy_VelLon_b)
    call HistoryPut('v_b'  ,xy_VelLat_b)
    call HistoryPut('h_b'  ,xy_Hsfc_b + HsfcAvr)
    call HistoryPut('vor_b',xy_Vor_b)
    call HistoryPut('div_b',xy_HDiv_b)

    call HistoryPut('u_bb'  ,xy_VelLon_bb)
    call HistoryPut('v_bb'  ,xy_VelLat_bb)
    call HistoryPut('h_bb'  ,xy_Hsfc_bb + HsfcAvr)
    call HistoryPut('vor_bb',xy_Vor_bb)
    call HistoryPut('div_bb',xy_HDiv_bb)
  end subroutine OutputRestart_gtool4

end program shallow_zd

