!----------------------------------------------------------------------
!     Copyright (c) 2005 Masatsugu Odaka. All rights reserved.
!----------------------------------------------------------------------
! Sample program for gt4f90io and ISPACK  2005/05/15 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 leap-frog & trapezoidal implicit 
! scheme. 
! 
program shallow_zd_bench01

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

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

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

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

  character(len=100) :: NameListFile='shallow_zd_bench01.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_Coli(:,:)          ! 惑星渦度
  real(8), allocatable :: xy_Topo(:,:)          ! 表面地形

  real(8), allocatable :: xy_VelLonCos_n(:,:)   ! xy_VelLon_n x Cos\phi
  real(8), allocatable :: xy_VelLatCos_n(:,:)   ! xy_VelLat_n x Cos\phi

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

  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(:)            ! 変位増分
  real(8), allocatable :: w_dVor(:)             ! 渦度増分
  real(8), allocatable :: w_dHDiv(:)            ! 発散増分

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

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

!----------------------------------------------------------------------
!  gtool4 用変数(複数ファイルに分割出力するため)
!----------------------------------------------------------------------

  type(GT_HISTORY) :: hst_2d, hst_1d, hst_rst

!----------------------------------------------------------------------
!  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) :: OutputFile1D=''            ! 出力ファイル名  
  character(len=100) :: OutputFile2D=''            ! 出力ファイル名  

  namelist /fileset/ InputRstFile, OutputRstfile, OutputFile1D, OutputFile2D

!----------------------------------------------------------------------
!  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.37122D6     ! 惑星半径
  real(8)            :: Omega    = 7.292D-5      ! 回転角速度
  real(8)            :: Grav     = 9.80616D0     ! 重力加速度
  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)            :: beta                    ! 陰解法用パラメタ
  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)            :: AngMomentAll              ! 全角運動量

  real(8)            :: CpuTime                 ! 計算 CPU 時間(frt 用)

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

!----------------------------------------------------------------------
!  内部変数 (Case1 実験用)
!----------------------------------------------------------------------

  real(8), parameter :: alpha0 = 0.0D0          ! 流れの軸の傾き
!  real(8), parameter :: alpha0 = 0.05D0         ! 流れの軸の傾き
!  real(8), parameter :: alpha0 = pi*0.5 - 0.05D0! 流れの軸の傾き
!  real(8), parameter :: alpha0 = pi*0.5       ! 流れの軸の傾き

  real(8)            :: L1norm                   ! L1 ノルム
  real(8)            :: L2norm                   ! L2 ノルム
  real(8)            :: Linfnorm                 ! L∞ノルム
  real(8)            :: Lon_ref = 1.5*pi         ! 解析解の中心経度
  real(8)            :: Lat_ref = 0.0            ! 解析解の中心緯度
  real(8)            :: VelLon_ref               ! 解析解の中心経度での速度
  real(8)            :: VelLat_ref               ! 解析解の中心経度での速度
  real(8)            :: DelLon_ref               ! 解析解の中心経度の増分
  real(8)            :: DelLat_ref               ! 解析解の中心経度の増分


  real(8),allocatable :: xy_HsfcRef(:,:)         ! 表面変位(参照解)
  real(8),allocatable :: xy_HsfcError(:,:)       ! 表面変位の誤差
  real(8),allocatable :: gamma(:,:)              ! case 1 初期値計算用座標

  integer            :: it
  integer            :: Ndiv = 100          ! 解析解中心座標計算時のループ回数

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

!----------------------------------------------------------------------
!  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))

  allocate(xy_VelLonCos_n(im,jm), xy_VelLatCos_n(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((nm+1)*(nm+1)), &
&          w_dVor((nm+1)*(nm+1)) , &
&          w_dHDiv((nm+1)*(nm+1)))

  allocate(xy_HsfcRef(im,jm), xy_HsfcError(im,jm), gamma(im,jm))

!----------------------------------------------------------------------
!  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)

  HsfcAvr  = AvrLonLat_xy(xy_Hsfc_n)
  GeoPotAvr = Grav*HsfcAvr 

  xy_Hsfc_n = xy_Hsfc_n - HsfcAvr
  xy_Hsfc_b = xy_Hsfc_b - HsfcAvr

  xy_VelLonCos_n = xy_VelLon_n*cos(xy_Lat)
  xy_VelLatCos_n = xy_VelLat_n*cos(xy_Lat)

  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                        ! ヒストリー初期化

!----------------------------------------------------------------------
! 参照解の保存: 初期値を保存しておく
!----------------------------------------------------------------------

  xy_HsfcRef = xy_Hsfc_n

!----------------------------------------------------------------------
!  浅水/順圧切替え
!  変数 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
    beta  = 0.0D0

    call MessageNotify('Message',     &
&                      'shallow_zd',  &
&                      'alpha=%f; %c', &
&                       d=(/alpha/),    &
&                       c1='explicit time integration')
  case default
    alpha = 1.0D0
    beta  = GeoPotAvr*DelTime**2/Rplanet**2

    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.') 

!----------------------------------------------------------------------
!  最初の 1 ステップ分の時間積分開始
!  InitTime = 0.0 の場合オイラースキームで計算
!  その他の場合はリープフロッグスキームで計算
!----------------------------------------------------------------------

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

  if (InitTime == 0.0) then

    call DbgMessage(fmt='%c InitTime=%f : %c', &
&                   d=(/InitTime/),            &
&                   c1=DbgMessageFmt,          &
&                   c2='Euler scheme is adapted.') 

!----------------------------------------------------------------------
!  増分 w_dHDiv, w_dVor, w_dHsfc を計算
!  * Case 1 では w_dHsfc のみ計算
!----------------------------------------------------------------------
!!
!!    call get_dVor(DelTime, xy_VelLonCos_n, xy_VelLatCos_n, xy_Vor_n, w_dVor)
!!
!!    call get_dHDiv(DelTime, xy_VelLonCos_n, xy_VelLatCos_n,      &
!!&                    xy_Hsfc_n, xy_Vor_n, w_Hsfc_n, w_dHDiv) 
!!

    call get_dHsfc(DelTime, xy_VelLonCos_n, xy_VelLatCos_n,      &
&                    xy_Hsfc_n, xy_Vor_n, w_HDiv_n, w_dHsfc) 

!----------------------------------------------------------------------
!  1 ステップ計算
!  * Case 1 では w_dHsfc のみ計算
!----------------------------------------------------------------------
!
!!    w_Vor_a = w_Vor_n + w_dVor + 2.0D0*DelTime*w_NumVis_w(w_Vor_n)
!!
!!    w_HDiv_a = ((1.0D0 + beta*rn(:,1))*w_HDiv_n +                  &
!!&                 w_dHDiv + DelTime*w_NumVis_w(w_HDiv_n))   &
!!&                 /(1.0D0 - beta*rn(:,1))

    w_Hsfc_a = ((1.0D0 + beta*rn(:,1))*w_Hsfc_n +                  &
&                 w_dHsfc + DelTime*w_NumVisScaler_w(w_Hsfc_n))& 
&                 /(1.0D0 - beta*rn(:,1))

!    w_Vor_b  = w_Vor_n  ; w_Vor_n  = w_Vor_a
!    w_HDiv_b = w_HDiv_n ; w_HDiv_n = w_hDiv_a
    w_HSfc_b = w_HSfc_n ; w_Hsfc_n = w_Hsfc_a

  else

!----------------------------------------------------------------------
!  増分 w_dHDiv, w_dVor, w_dHsfc を計算
!  変数 IntScheme で切替え. デフォルトは semi-implicit
!  * Case 1 では w_dHsfc のみ計算, Intscheme は 'explicit' を指定
!----------------------------------------------------------------------
!
!!    call get_dVor(2.0*DelTime, xy_VelLonCos_n, xy_VelLatCos_n, xy_Vor_n, w_dVor)

    select case (IntScheme)
    case ('explicit') 
!!      call get_dHDiv(2.0*DelTime, xy_VelLonCos_n, xy_VelLatCos_n,      &
!!&                    xy_Hsfc_n, xy_Vor_n, w_Hsfc_n, w_dHDiv) 

      call get_dHsfc(2.0*DelTime, xy_VelLonCos_n, xy_VelLatCos_n,      &
&                    xy_Hsfc_n, xy_Vor_n, w_HDiv_n, w_dHsfc) 
    case default
!!      call get_dHDiv(2.0*DelTime, xy_VelLonCos_n, xy_VelLatCos_n,      &
!!&                    xy_Hsfc_n, xy_Vor_n, w_Hsfc_b, w_dHDiv)  

      call get_dHsfc(2.0D0*DelTime, xy_VelLonCos_n, xy_VelLatCos_n,      &
&                    xy_Hsfc_n, xy_Vor_n, w_HDiv_b, w_dHsfc) 
    end select

!----------------------------------------------------------------------
!  1 ステップ計算
!  * Case 1 では w_Hsfc のみ計算
!----------------------------------------------------------------------
!
!!    w_Vor_a = w_Vor_b + w_dVor + 2.0D0*DelTime*w_NumVis_w(w_Vor_b)
!!
!!    w_HDiv_a = ((1.0D0 + beta*rn(:,1))*w_HDiv_b +                  &
!!&                 w_dHDiv + 2.0D0*DelTime*w_NumVis_w(w_HDiv_b))   &
!!&                 /(1.0D0 - beta*rn(:,1))

    w_Hsfc_a = ((1.0D0 + beta*rn(:,1))*w_Hsfc_b +                  &
&                 w_dHsfc + 2.0D0*DelTime*w_NumVisScaler_w(w_Hsfc_b))& 
&                 /(1.0D0 - beta*rn(:,1))

!----------------------------------------------------------------------
!  時間フィルターをかけて, 変数を 1 ステップ更新する.
!  * Case 1 では w_Vor と w_Hdiv は定常
!----------------------------------------------------------------------
!
!!    call timefilter(w_Vor_a , w_Vor_n , w_Vor_b )   
!!    call timefilter(w_HDiv_a, w_HDiv_n, w_HDiv_b)      
    call timefilter(w_Hsfc_a, w_Hsfc_n, w_Hsfc_b)   

  end if

!---------------------------------------------------------------------
!  Case 1 の解析解のための計算
!  * 山の中心座標の時間発展をオイラー法で計算する.
!---------------------------------------------------------------------

  VelLon_ref =  2.0D0*pi/(12.0D0*86400.0D0) * &
&             (cos(alpha0) + tan(Lat_ref)*cos(Lon_ref)*sin(alpha0)) 

  VelLat_ref = - 2.0D0*pi/(12.0D0*86400.0D0) * &
&                   sin(Lon_ref)*sin(alpha0)

  Lon_ref = Lon_ref + DelTime*VelLon_ref 
  Lat_ref = Lat_ref + DelTime*VelLat_ref 

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

!!  call uv(w_Vor_n, w_HDiv_n, xy_VelLon_n, xy_VelLat_n) 
!!  call uvcos(w_Vor_n, w_HDiv_n, xy_VelLonCos_n, xy_VelLatCos_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)

!!  call uv(w_Vor_b, w_HDiv_b, xy_VelLon_b, xy_VelLat_b) 

!!  xy_Vor_b  = xy_w(w_Vor_b )
!!  xy_HDiv_b = xy_w(w_HDiv_b)
   xy_Hsfc_b = xy_w(w_Hsfc_b)
              
!----------------------------------------------------------------------
!  2 ステップ以降の時間積分開始
!----------------------------------------------------------------------

  do Nstep = 2, NstepALL

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

!----------------------------------------------------------------------
!  増分 w_dHDiv, w_dVor, w_dHsfc を計算
!  変数 IntScheme で切替え. デフォルトは semi-implicit
!  * Case 1 では w_dHsfc のみ計算, Intscheme は 'explicit' を指定
!----------------------------------------------------------------------
!
!!    call get_dVor(2.0*DelTime, xy_VelLonCos_n, xy_VelLatCos_n, xy_Vor_n, w_dVor)

    select case (IntScheme)
    case ('explicit') 
!!      call get_dHDiv(2.0*DelTime, xy_VelLonCos_n, xy_VelLatCos_n,      &
!!&                    xy_Hsfc_n, xy_Vor_n, w_Hsfc_n, w_dHDiv) 

      call get_dHsfc(2.0*DelTime, xy_VelLonCos_n, xy_VelLatCos_n,      &
&                    xy_Hsfc_n, xy_Vor_n, w_HDiv_n, w_dHsfc) 
    case default
!!      call get_dHDiv(2.0*DelTime, xy_VelLonCos_n, xy_VelLatCos_n,      &
!!&                    xy_Hsfc_n, xy_Vor_n, w_Hsfc_b, w_dHDiv)  

      call get_dHsfc(2.0D0*DelTime, xy_VelLonCos_n, xy_VelLatCos_n,      &
&                    xy_Hsfc_n, xy_Vor_n, w_HDiv_b, w_dHsfc) 
    end select

!----------------------------------------------------------------------
!  1 ステップ計算
!  * Case 1 では w_Hsfc のみ計算
!----------------------------------------------------------------------
!
!!    w_Vor_a = w_Vor_b + w_dVor + 2.0D0*DelTime*w_NumVis_w(w_Vor_b)
!!
!!    w_HDiv_a = ((1.0D0 + beta*rn(:,1))*w_HDiv_b +                  &
!!&                 w_dHDiv + 2.0D0*DelTime*w_NumVis_w(w_HDiv_b))   &
!!&                 /(1.0D0 - beta*rn(:,1))

    w_Hsfc_a = ((1.0D0 + beta*rn(:,1))*w_Hsfc_b +                  &
&                 w_dHsfc + 2.0D0*DelTime*w_NumVisScaler_w(w_Hsfc_b))& 
&                 /(1.0D0 - beta*rn(:,1))

!---------------------------------------------------------------------
!  Case 1 の解析解のための計算
!  * 山の中心座標の時間発展をオイラー法で計算する.
!  * 計算精度を上げるため時間刻を細かくとる
!---------------------------------------------------------------------

    do it = 1, Ndiv

    VelLon_ref =  2.0D0*pi/(12.0D0*86400.0D0) * &
&              (cos(alpha0) + tan(Lat_ref)*cos(Lon_ref)*sin(alpha0)) 

    VelLat_ref = - 2.0D0*pi/(12.0D0*86400.0D0) * &
&                   sin(Lon_ref)*sin(alpha0)

    Lon_ref = Lon_ref + DelTime/Ndiv*VelLon_ref 
    Lat_ref = Lat_ref + DelTime/Ndiv*VelLat_ref 

!-- 球面三角法を使って計算する場合 -- 
!-- 精度がよくないようなので使わない(2005/06/20 Odaka)
!
!    DelLat_ref = asin( sin(2.0D0*pi/(12.0D0*86400.0D0))*sin(alpha0) )
!    DelLon_ref = asin( sin(2.0D0*pi/(12.0D0*86400.0D0))*cos(alpha0)/cos(DelLat_ref) )
!    if ( cos(Lon_ref) < 0.0D0 ) then
!      DelLat_ref = - DelLat_ref
!    end if 
!
!    Lon_ref = Lon_ref + DelTime/Ndiv*DelLon_ref 
!    Lat_ref = Lat_ref + DelTime/Ndiv*DelLat_ref 
!
!-------------------------------------
    end do 

!----------------------------------------------------------------------
!  時間フィルターをかけて, 変数を 1 ステップ更新する.
!  * Case 1 では w_Vor と w_Hdiv は定常
!----------------------------------------------------------------------
!
!!    call timefilter(w_Vor_a , w_Vor_n , w_Vor_b )   
!!    call timefilter(w_HDiv_a, w_HDiv_n, w_HDiv_b)      
    call timefilter(w_Hsfc_a, w_Hsfc_n, w_Hsfc_b)   

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

!!    call uv(w_Vor_n, w_HDiv_n, xy_VelLon_n, xy_VelLat_n) 
!!    call uvcos(w_Vor_n, w_HDiv_n, xy_VelLonCos_n, xy_VelLatCos_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)

!!    call uv(w_Vor_b, w_HDiv_b, xy_VelLon_b, xy_VelLat_b) 

!!    xy_Vor_b  = xy_w(w_Vor_b)
!!    xy_HDiv_b = xy_w(w_HDiv_b)
    xy_Hsfc_b = xy_w(w_Hsfc_b)
              
!----------------------------------------------------------------------
!  出力
!  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 +         &
&                     Grav*(xy_Hsfc_n+HsfcAvr - xy_Topo))*0.5D0 )

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


!-- 解析解と誤差の計算

      gamma = Rplanet*acos(sin(Lat_ref)*sin(xy_Lat) + &
&                 cos(Lat_ref)*cos(xy_Lat)*cos(xy_Lon - Lon_ref))

      where ( gamma < Rplanet/3 ) 
        xy_HsfcRef = 0.5D0*1000D0*(1.0 + cos(3*pi*gamma/Rplanet))
      elsewhere
        xy_HsfcRef = 0.0D0
      end where

      xy_HsfcError = xy_Hsfc_n + HsfcAvr - xy_Topo - xy_HsfcRef

      L1norm = &
&       IntLonLat_xy(abs(xy_HsfcError))/IntLonLat_xy(abs(xy_HsfcRef))

      L2norm = &
&       sqrt(IntLonLat_xy(xy_HsfcError**2))/sqrt(IntLonLat_xy(xy_HsfcRef**2))

      Linfnorm = &
&       maxval(abs(xy_HsfcError))/maxval(abs(xy_HsfcRef))

!!      call MessageNotify('Message',      &
!!&                        'shallow_zd',   &
!!&                        'Time=%f', d=(/Time/))

!-- 出力

      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 OutputRestart_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_Coli    , xy_Topo      )

  deallocate(xy_VelLonCos_n, xy_VelLatCos_n)

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

  deallocate(xy_HsfcRef, xy_HsfcError,gamma)

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

!-- CPU 時間の出力 (Fujitsu frt 用)
!  call clock(CpuTime, 0, 2)
!
!      call MessageNotify('Message',      &
!&                        'shallow_zd_benc01',   &
!&                        'CpuTime=%f', d=(/CpuTime/))

  stop

contains

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

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

  subroutine get_dVor(dt, xy_VelLonCos, xy_VelLatCos, xy_Vor, w_dVor)

    real(8),intent(in) :: dt    
    real(8),intent(in) :: xy_VelLonCos(im,jm)
    real(8),intent(in) :: xy_VelLatCos(im,jm)
    real(8),intent(in) :: xy_Vor(im,jm)
    real(8),intent(out):: w_dVor((nm+1)*(nm+1))

    call BeginSub('get_dVor')

    w_dVor = dt*(                                               &
&     - w_DivLambdaMu_xy_xy((xy_Coli + xy_Vor)*xy_VelLonCos,    &
&                   (xy_Coli + xy_Vor)*xy_VelLatCos) / Rplanet)

    call EndSub('get_dVor')

  end subroutine get_dVor

!----------------------------------------------------------------------
!  get_dHDiv
!  発散の式の増分 w_dHDiv を計算する
!----------------------------------------------------------------------

  subroutine get_dHDiv(dt, xy_VelLonCos, xy_VelLatCos, xy_Hsfc, &
&                        xy_Vor, w_Hsfc, w_dHDiv)

    real(8),intent(in) :: dt
    real(8),intent(in) :: xy_VelLonCos(im,jm)
    real(8),intent(in) :: xy_VelLatCos(im,jm)
    real(8),intent(in) :: xy_Hsfc(im,jm)
    real(8),intent(in) :: xy_Vor(im,jm)
    real(8),intent(in) :: w_Hsfc((nm+1)*(nm+1))
    real(8),intent(out):: w_dHDiv((nm+1)*(nm+1))

    call BeginSub('get_dHDiv')

    w_dHDiv = dt*(                                                    &
&     w_DivLambdaMu_xy_xy(  (xy_Coli + xy_Vor)*xy_VelLatCos,          &
&                         - (xy_Coli + xy_Vor)*xy_VelLonCos)/Rplanet  &
&     - w_Lapla_w(Grav*w_Hsfc                                         &
&        + 0.5D0*w_xy((xy_VelLonCos**2 + xy_VelLatCos**2)/cos(xy_Lat)**2)) &
&          /Rplanet**2                                                &
&     + alpha*DelTime*Grav*w_Lapla_w(                                 &
&       w_DivLambdaMu_xy_xy((xy_Hsfc - xy_Topo)*xy_VelLonCos,         &
&                   (xy_Hsfc - xy_Topo)*xy_VelLatCos)/Rplanet)/Rplanet**2)

    call EndSub('get_dHDiv')

  end subroutine get_dHDiv

!----------------------------------------------------------------------
!  get_dHsfc
!  質量保存の式の増分 w_dHsfc を計算する
!----------------------------------------------------------------------

  subroutine get_dHsfc(dt, xy_VelLonCos, xy_VelLatCos, xy_Hsfc, &
&                        xy_Vor, w_HDiv, w_dHsfc)

    real(8),intent(in) :: dt
    real(8),intent(in) :: xy_VelLonCos(im,jm)
    real(8),intent(in) :: xy_VelLatCos(im,jm)
    real(8),intent(in) :: xy_Hsfc(im,jm)
    real(8),intent(in) :: xy_Vor(im,jm)
    real(8),intent(in) :: w_HDiv((nm+1)*(nm+1))
    real(8),intent(out):: w_dHsfc((nm+1)*(nm+1))

    call BeginSub('get_dHsfc')

    w_dHsfc = dt*(                                                    &
&     - w_DivLambdaMu_xy_xy((xy_Hsfc - xy_Topo)*xy_VelLonCos,         &
&                           (xy_Hsfc - xy_Topo)*xy_VelLatCos)/Rplanet &       
&     - HsfcAvr*w_HDiv                                                &
&     - alpha*DelTime*HsfcAvr*(                                       &
&         w_DivLambdaMu_xy_xy( (xy_Coli + xy_Vor)*xy_VelLatCos,       &
&                     - (xy_Coli + xy_Vor)*xy_VelLonCos)/Rplanet      &
&         - w_Lapla_w(                                                &
&             w_xy((xy_VelLonCos**2 + xy_VelLatCos**2)/cos(xy_Lat)**2)*0.5D0)&
&             /Rplanet**2))

    call EndSub('get_dHsfc')
     
  end subroutine get_dHsfc

!----------------------------------------------------------------------
!  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

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

  subroutine uvcos(w_Vor, w_HDiv, xy_VelLonCos, xy_VelLatCos)

    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_VelLonCos(im,jm)
    real(8),intent(out) :: xy_VelLatCOs(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_VelLonCos = (delta*xy_GradLambda_w(w_VelPot)         &
&                  - xy_GradMu_w(w_Stream))/Rplanet

    xy_VelLatCos = (xy_GradLambda_w(w_Stream)              &
&                  + delta*xy_GradMu_w(w_VelPot))/Rplanet

    call EndSub('uv')

  end subroutine uvcos

!----------------------------------------------------------------------
!  timefilter
!  Asselin (1972) の時間フィルターを適用し, 変数を 1 ステップ更新する.
!----------------------------------------------------------------------

  subroutine timefilter(w_data_a, w_data_n, w_data_b)

    real(8),intent(inout) :: 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('timefilter')

  !----- 時間フィルター -----

    w_data_n = w_data_n &
&     + 0.5D0*TfilCoef*(w_data_b - 2.0D0*w_data_n + w_data_a)

  !----- 1 ステップ進める -----
    w_data_b = w_data_n
    w_data_n = w_data_a 

    call EndSub('timefilter')

  end subroutine timefilter

!----------------------------------------------------------------------
!  w_DivLambdaMu_xy_xy
!  \lamda, \mu 座標 を用いた発散の計算
!----------------------------------------------------------------------
  function w_DivLambdaMu_xy_xy(xy_dataLonCos, xy_dataLatCos)

    real(8)             :: w_DivLambdaMu_xy_xy((nm+1)*(nm+1))
    real(8), intent(in) :: xy_dataLonCos(im,jm)
    real(8), intent(in) :: xy_dataLatCos(im,jm)

    call BeginSub('w_DivLambdaMu_xy_xy')
  
    w_DivLambdaMu_xy_xy = &
&      w_DivLambda_xy(xy_dataLonCos) + w_DivMu_xy(xy_dataLatCos)

    call EndSub('w_DivLambdaMu_xy_xy')

  end function w_DivLambdaMu_xy_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 ) 

    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) ! 発散

    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)

  end subroutine InputRestart_gtool4

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

  subroutine Output_gtool4_init
    call HistoryCreate( &                               ! ヒストリー作成
&          file=trim(OutputFile2D),                                 & 
&          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),&
&          history=hst_2d)

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

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

    call HistoryPut('topo',xy_Topo, hst_2d)

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

    call HistoryAddVariable( &                          ! 変数定義
&          varname='h_ref', dims=(/'lon','lat','t  '/), & 
&          longname='height(reference)',       &
&          units='m', xtype='double', history=hst_2d)
    call HistoryAddVariable( &                          ! 変数定義
&          varname='h_error', dims=(/'lon','lat','t  '/), & 
&          longname='height error',         &
&          units='m', xtype='double', history=hst_2d)

    call HistoryCreate( &                               ! ヒストリー作成
&          file=trim(OutputFile1D),                               & 
&          title=trim(ExpTitle),                                  &
&          source=trim(ExpSrc),                                   &
&          institution=trim(ExpInst),                             &
&          dims=(/'t  '/), dimsizes=(/0/),                        &
&          longnames=(/'time     '/),     &
&          units=(/'sec.'/),                                      &
&          origin=real(InitTime), interval=real(NStepOut*DelTime),&
&          history=hst_1d)

    call HistoryAddVariable( &                          ! 変数定義
&          varname='mass', dims=(/'t  '/), & 
&          longname='total mass',          &
&          units='m^3', xtype='double', history=hst_1d)
    call HistoryAddVariable( &                          ! 変数定義
&          varname='energy_all', dims=(/'t  '/), & 
&          longname='total energy',          &
&          units='m^5/s^2', xtype='double', history=hst_1d)
    call HistoryAddVariable( &                          ! 変数定義
&          varname='angmoment_all', dims=(/'t  '/), & 
&          longname='total angular momentum', &
&          units='m^5/s', xtype='double', history=hst_1d)

    call HistoryAddVariable( &                          ! 変数定義
&          varname='l1norm', dims=(/'t  '/), & 
&          longname='L1 mass norm',          &
&          units=' ', xtype='double', history=hst_1d)
    call HistoryAddVariable( &                          ! 変数定義
&          varname='l2norm', dims=(/'t  '/), & 
&          longname='L_2 mass norm', units=' ', xtype='double', &
&          history=hst_1d)
    call HistoryAddVariable( &                          ! 変数定義
&          varname='linfnorm', dims=(/'t  '/), & 
&          longname='L_inf mass norm',         &
&          units=' ', xtype='double', history=hst_1d)


  end subroutine Output_gtool4_init

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

  subroutine Output_gtool4
    call HistoryPut('t'      ,Time             ,hst_2d)
    call HistoryPut('u'      ,xy_VelLon_n      ,hst_2d)
    call HistoryPut('v'      ,xy_VelLat_n      ,hst_2d)
    call HistoryPut('h'      ,xy_Hsfc_n+HsfcAvr,hst_2d)
    call HistoryPut('vor'    ,xy_Vor_n         ,hst_2d)
    call HistoryPut('div'    ,xy_HDiv_n        ,hst_2d)
    call HistoryPut('h_ref'  ,xy_HsfcRef      ,hst_2d)
    call HistoryPut('h_error',xy_HsfcError     ,hst_2d)

    call HistoryPut('t'       ,Time            ,hst_1d)
    call HistoryPut('l1norm'  ,L1norm          ,hst_1d)
    call HistoryPut('l2norm'  ,L2norm          ,hst_1d)
    call HistoryPut('linfnorm',Linfnorm        ,hst_1d)
    call HistoryPut('mass'    ,MassAll         ,hst_1d)
    call HistoryPut('energy_all'  ,EnergyAll   ,hst_1d)
    call HistoryPut('angmoment_all' ,AngMomentAll,hst_1d)

  end subroutine Output_gtool4

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

  subroutine Output_gtool4_close
    call HistoryClose(history=hst_1d)
    call HistoryClose(history=hst_2d)
  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),                 &
&          history=hst_rst)

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

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

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

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

  end subroutine OutputRestart_gtool4_init

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

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

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

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

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

  subroutine OutputRestart_gtool4_close
    call HistoryClose(history=hst_rst)
  end subroutine OutputRestart_gtool4_close


end program shallow_zd_bench01

