!----------------------------------------------------------------------
!     Copyright (c) 2003 Masatsugu Odaka. All rights reserved.
!----------------------------------------------------------------------
! Sample program for gt4f90io and ISPACK  2005/06/24 M. Odaka
! 
! Making initial condition for 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 initial conditions are similar to those of Williamson et al. (1992) 
! [J. Comp. Phys., 102, 211-224]: 
!   case 1: advection of cosine bell over Pole
!   case 2: global steady state nonlinear zonal geostrophic flow
!   case 3: steady state nonlinear zonal geostrophic flow with conpact support
!   case 4: Forced nonlinear system with a translating low
!   case 5: zonal flow over an isolated mountain
!   case 6: Rossby-Haurwitz wave
!
program shallow_zd_init

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

  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_init.nml'

!----------------------------------------------------------------------
!  変数
!----------------------------------------------------------------------

  real(8), allocatable :: xy_VelLon(:,:)    ! 速度経度成分
  real(8), allocatable :: xy_VelLat(:,:)    ! 速度緯度成分
  real(8), allocatable :: xy_Hsfc(:,:)      ! 表面変位
  real(8), allocatable :: xy_Vor(:,:)       ! 渦度
  real(8), allocatable :: xy_HDiv(:,:)      ! 発散
  real(8), allocatable :: xy_Stream(:,:)    ! 流線関数
  real(8), allocatable :: xy_Topo(:,:)      ! 表面地形

  real(8), allocatable :: xy_VelLonOrg(:,:) ! 座標変換前の速度経度成分(case3)
  real(8), allocatable :: xy_VelLatOrg(:,:) ! 座標変換前の速度緯度成分(case3)
  real(8), allocatable :: xy_LonOrg(:,:)    ! 座標変換前の経度(case3)
  real(8), allocatable :: xy_LatOrg(:,:)    ! 座標変換前の緯度(case3)
  real(8), allocatable :: xy_VelLonMean(:,:)! 平均速度経度成分(case4)


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

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

  real(8), parameter  :: Tinit = 0.0D0            ! 初期時刻

!----------------------------------------------------------------------
!  初期条件パラメタ 
!----------------------------------------------------------------------

!---- 初期条件パラメター (case1,2,3,5) ---- 
  real(8), parameter :: alpha = 0.0D0          ! 流れの軸の傾き
!  real(8), parameter :: alpha = 0.05D0         ! 流れの軸の傾き
!  real(8), parameter :: alpha = pi*0.25        ! 流れの軸の傾き
!  real(8), parameter :: alpha = pi/3.0D0       ! 流れの軸の傾き
!  real(8), parameter :: alpha = pi*0.5 - 0.05D0 ! 流れの軸の傾き
!  real(8), parameter :: alpha = pi*0.5         ! 流れの軸の傾き

!---- 初期条件パラメター (case1) ---- 
  real(8), parameter :: Lon_ref = pi*1.5
  real(8), parameter :: Lat_ref = 0.0D0

!---- 初期条件パラメター (case3) ---- 
  real(8), parameter :: Lat_b = - pi/6.0    
  real(8), parameter :: Lat_e = pi*0.5
  real(8), parameter :: Xe0   = 0.3D0

!---- 初期条件パラメター (case4) ---- 
  real(8), parameter :: Lon_0 = 0.0D0
  real(8), parameter :: Lat_0 = pi/4.0
  real(8), parameter :: sigma = (12.74244D0)**2

!---- 初期条件パラメター (case5) ---- 

  real(8), parameter :: TopoMax  = 2.0D3          ! 山頂高度
  real(8), parameter :: Lat_topo = pi/6.0         ! 山頂緯度
  real(8), parameter :: Lon_topo = pi*0.5         ! 山頂経度
  real(8), parameter :: R_topo   = pi/9.0         ! 地形の半径

!---- 初期条件パラメター (case6) ----
  real(8), parameter :: Omega0 = 7.848D-6        ! 初期角速度
  real(8), parameter :: Wnum = 4.0D0             ! 初期波数

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

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

  namelist /expset/  ExpTitle, ExpInst, ExpSrc, ExpModel, ExpCase

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

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

  namelist /fileset/ OutputRstfile

!----------------------------------------------------------------------
!  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)            :: HsfcAvr                  ! 平均変位

  namelist /paramset/ Rplanet, Omega, Grav

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

  logical            :: DebugOn = .false.        ! デバッグ出力スイッチ
  namelist /debugset/ DebugOn

!----------------------------------------------------------------------
!  作業変数 
!----------------------------------------------------------------------
  integer             :: jy                      ! ループ変数
  real(8),allocatable :: gamma(:,:)              ! case 1,3,4 初期値計算用座標
  real(8),allocatable :: Xe(:,:)                 ! case 3 初期値計算用座標
  real(8)             :: VelLon0                 ! 初期速度


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

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

  open(10,file=trim(NameListFile))
  read(10,nml=expset)  
  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(im,jm), xy_VelLat(im,jm), xy_Hsfc(im,jm)  , &
&          xy_Vor(im,jm)   , xy_HDiv(im,jm)  , xy_Stream(im,jm), &
&          xy_Topo(im,jm))

  allocate(xy_LonOrg(im,jm), xy_LatOrg(im,jm), &
&          xy_VelLonOrg(im,jm), xy_VelLatOrg(im,jm), &
&          gamma(im,jm), Xe(im,jm), xy_VelLonMean(im,jm))

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

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


!== 初期値の計算 ======================================================
!  パラメータ ExpCase に与えた文字列ごとに初期値を用意する.
!
  call DbgMessage(fmt='call %c', c1='w_initial') 
  select case (ExpCase)

!----------------------------------------------------------------------
!  case 1: cos 型 bell の移流
!----------------------------------------------------------------------

  case ('case1')                                 

    call MessageNotify('Message',         &
&                      'shallow_zd_init', &
&                      'Case 1: advection of cosine bell over Pole')

    HsfcAvr = 1.0D3                                   ! 平均変移
    VelLon0 = 2.0D0*pi*Rplanet/(12.0D0*86400.0D0)     ! 初期速度

    xy_Stream = - Rplanet*VelLon0*(sin(xy_Lat)*cos(alpha)                &
&                                  - cos(xy_Lon)*cos(xy_Lat)*sin(alpha))  

    xy_VelLon = VelLon0*(cos(xy_Lat)*cos(alpha)  &
&                          + cos(xy_Lon)*sin(xy_Lat)*sin(alpha))

    xy_VelLat = - VelLon0*sin(xy_Lon)*sin(alpha)

    xy_Vor  = xy_w(w_Lapla_w(w_xy(xy_Stream)))/Rplanet**2 

    xy_HDiv = 0.0D0

!    xy_HDiv  = xy_w(w_DivLambdaMu_xy_xy  &
!&      (xy_VelLon*cos(xy_Lat), xy_VelLat*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_Hsfc = 0.5D0*HsfcAvr*(1.0 + cos(3*pi*gamma/Rplanet))
    elsewhere
      xy_Hsfc = 0.0D0
    end where

    xy_Topo = 0.0D0

!----------------------------------------------------------------------
!  case 2: 非線形帯状地衡流 1
!----------------------------------------------------------------------

  case ('case2')                                 
    call MessageNotify('Message',         &
&                      'shallow_zd_init', &
&        'Case 2: global steady state nonlinear zonal geostrophic flow')

    HsfcAvr = 2.94D4/Grav                           ! 平均変移
    VelLon0 = 2.0D0*pi*Rplanet/(12.0*86400.0)       ! 初期速度

    xy_Stream = - Rplanet*VelLon0*(sin(xy_Lat)*cos(alpha) &
&                                    - cos(xy_Lon)*cos(xy_Lat)*sin(alpha))

    xy_VelLon = VelLon0*(cos(xy_Lat)*cos(alpha)  &
&                          + cos(xy_Lon)*sin(xy_Lat)*sin(alpha))

    xy_VelLat = - VelLon0*sin(xy_Lon)*sin(alpha)

    xy_Vor  = xy_w(w_Lapla_w(w_xy(xy_Stream)))/Rplanet**2 

    xy_Hdiv = 0.0D0

    xy_Hsfc = HsfcAvr - (Rplanet*Omega*VelLon0 + 0.5D0*VelLon0**2)* &
&              (- cos(xy_Lon)*cos(xy_Lat)*sin(alpha)                &
&                   + sin(xy_Lat)*cos(alpha))**2/Grav

    xy_Topo = 0.0D0

!----------------------------------------------------------------------
!  case 3: 非線形帯状地衡流 2
!----------------------------------------------------------------------

  case ('case3')                                 
    call MessageNotify('Message',      &
&                      'shallow_zd_init', &
& 'Case 3: steady state nonlinear zonal geostrophic flow with conpact support')

    HsfcAvr = 2.94D4 / Grav                           ! 平均変移
    VelLon0 = 2.0D0*pi*Rplanet/(12.0*86400.0)         ! 初期速度

    xy_VelLonOrg = 0.0D0
    xy_VelLatOrg = 0.0D0

    xy_VelLon = 0.0D0
    xy_VelLat = 0.0D0

  !-- 精度の都合上, α=0 (流れの傾きのない場合)とその他の場合で
  !   計算方法を変更する.

    if (alpha .eq. 0.0D0) then  

  !-- 座標変換前の緯度経度を計算する
  !-- α=0 の場合, 座標変換前と変換度の緯度経度を一致させる

      xy_LatOrg = xy_Lat
      xy_LonOrg = xy_Lon

    else

  !-- 座標変換前の緯度経度を計算する

      xy_LatOrg = asin(sin(xy_Lat)*cos(alpha) &
&                        - cos(xy_Lat)*cos(xy_Lon)*sin(alpha))

      gamma = cos(alpha)*cos(xy_Lon)*cos(xy_Lat) + sin(alpha)*sin(xy_Lat)

      where (gamma .ge. 0.0D0) 
        xy_LonOrg = asin(sin(xy_Lon)*cos(xy_Lat)/cos(xy_LatOrg))
      elsewhere
        xy_LonOrg = pi - asin(sin(xy_Lon)*cos(xy_Lat)/cos(xy_LatOrg))
      end where
 
    end if

  !-- 速度場の計算

    Xe = Xe0*(xy_LatOrg - Lat_b)/(Lat_e - Lat_b)

    where (Xe .gt. 0.0D0) 
      xy_VelLonOrg = VelLon0*exp(-1.0/Xe)*exp(-1.0/(Xe0 - Xe))*exp(4.0/Xe0)
    end where

    xy_VelLat = - xy_VelLonOrg*sin(alpha)*sin(xy_LonOrg)/cos(xy_Lat)
   
    where (cos(xy_Lon) .ne. 0.0D0)
      xy_VelLon = xy_VelLat*sin(xy_Lat)*sin(xy_Lon)/cos(xy_Lon) + &
&                   xy_VelLonOrg*cos(xy_LonOrg)/cos(xy_Lon)
    end where

  !-- 台形公式で h を計算

    xy_Hsfc(:,1) = HsfcAvr
    do jy = 2, jm
      xy_Hsfc(:,jy) = xy_Hsfc(:,jy-1) - Rplanet/Grav*                     &
&                   ((2.0*Omega*sin(xy_LatOrg(:,jy))                      &
&                       + xy_VelLonOrg(:,jy)*tan(xy_LatOrg(:,jy))/Rplanet)&
&                    *xy_VelLonOrg(:,jy) +                                &
&                    (2.0*Omega*sin(xy_LatOrg(:,jy-1))                    & 
&                       + xy_VelLonOrg(:,jy-1)*tan(xy_LatOrg(:,jy-1))/Rplanet)&
&                    *xy_VelLonOrg(:,jy-1))                               &
&                   * (xy_LatOrg(:,jy) - xy_LatOrg(:,jy-1))*0.5
    end do

    xy_Vor = (xy_GradLon_w(w_xy(xy_VelLat)) - &
&              xy_GradLat_w(w_xy(xy_VelLon)))/Rplanet

    xy_HDiv  = xy_w(w_DivLambdaMu_xy_xy  &
&      (xy_VelLon*cos(xy_Lat), xy_VelLat*cos(xy_Lat)))/Rplanet

    xy_Topo = 0.0D0

!----------------------------------------------------------------------
!  case 4: 移動性低気圧を持つ強制された非線形の流れ
!----------------------------------------------------------------------

  case ('case4') 

    call MessageNotify('Message',         &
&                      'shallow_zd_init', &
&                    'Case 4: Forced nonlinear system with a translating low')


    VelLon0 = 40.0D0     ! 平均速度
    HsfcAvr = 1.0D5/Grav ! 平均変移

    gamma = sin(Lat_0)*sin(xy_Lat) + cos(Lat_0)*cos(xy_Lat)*cos( &
&             xy_Lon - Lon_0)
    xy_Stream = - 0.03D0*Grav*HsfcAvr/(2.0D0*Omega*sin(Lat_0)) &
&                   *exp( - sigma*(1 - gamma)/(1 + gamma) )
  
    xy_VelLonMean = VelLon0*(sin(2.0D0*xy_Lat))**14 

    xy_VelLon = xy_VelLonMean - xy_GradLat_w(w_xy(xy_Stream))/Rplanet
    xy_VelLat = xy_GradLon_w(w_xy(xy_Stream))/Rplanet

  !-- 台形公式で h を計算

    xy_Hsfc(:,1) = HsfcAvr
    do jy = 2, jm
      xy_Hsfc(:,jy) = xy_Hsfc(:,jy-1) - Rplanet/Grav*                   &
&                   ((2.0*Omega*sin(xy_Lat(:,jy))                       &
&                       + xy_VelLonMean(:,jy)*tan(xy_Lat(:,jy))/Rplanet)&
&                    *xy_VelLonMean(:,jy) +                             &
&                    (2.0*Omega*sin(xy_Lat(:,jy-1))                     & 
&                       + xy_VelLonMean(:,jy-1)*tan(xy_Lat(:,jy-1))/Rplanet)&
&                    *xy_VelLonMean(:,jy-1))                            &
&                   * (xy_Lat(:,jy) - xy_Lat(:,jy-1))*0.5
    end do

    xy_Hsfc = xy_Hsfc + 2.0D0*Omega*sin(xy_Lat)*xy_Stream/Grav

    xy_Vor = (xy_GradLon_w(w_xy(xy_VelLat)) - &
&              xy_GradLat_w(w_xy(xy_VelLon)))/Rplanet

    xy_HDiv  = xy_w(w_DivLambdaMu_xy_xy  &
&      (xy_VelLon*cos(xy_Lat), xy_VelLat*cos(xy_Lat)))/Rplanet

    xy_Topo = 0.0D0

!----------------------------------------------------------------------
!  case 5: 山岳応答問題 
!----------------------------------------------------------------------

  case ('case5')                          
    call MessageNotify('Message',      &
&                      'shallow_zd_init', &
&                      'Case 5: zonal flow over an isolated mountain')

    VelLon0 = 20.0D0                                  ! 初期速度
    HsfcAvr = 5.96D3                                  ! 平均変位

    xy_Stream = - VelLon0*Rplanet*(sin(xy_Lat)*cos(alpha) -             &
&                                  cos(xy_Lon)*cos(xy_Lat)*sin(alpha))  

    xy_VelLon = - xy_GradLat_w(w_xy(xy_Stream))/Rplanet
    xy_VelLat =   xy_GradLon_w(w_xy(xy_Stream))/Rplanet

    xy_Vor  = xy_w(w_Lapla_w(w_xy(xy_Stream)))/Rplanet**2 

    xy_HDiv  = xy_w(w_DivLambdaMu_xy_xy  &
&      (xy_VelLon*cos(xy_Lat), xy_VelLat*cos(xy_Lat)))/Rplanet

    xy_Hsfc  = HsfcAvr - (Rplanet*Omega*VelLon0 + VelLon0**2*0.5)/Grav* &
&                        (- cos(xy_Lon)*cos(xy_Lat)*sin(alpha)          &
&                         + sin(xy_Lat)*cos(alpha))**2 

    xy_Topo = TopoMax * (1.0D0 - sqrt(min(R_topo**2,              &
&                                  (xy_Lon-Lon_topo)**2 +         &
&                                  (xy_Lat-Lat_topo)**2))/R_topo)


!----------------------------------------------------------------------
!  case 6: Rossby-Haurwitz 波 
!----------------------------------------------------------------------

  case ('case6')                                 
    call MessageNotify('Message',      &
&                      'shallow_zd_init', &
&                      'Case 6: Rossby-Haurwitz wave')

    VelLon0 = 0.0D0                                   ! 初期風速
    HsfcAvr = 8.0D3                                   ! 平均変位

    xy_Stream = - Rplanet**2*Omega0*sin(xy_Lat)*(                 &
&               1.0D0 - cos(xy_Lat)**Wnum*cos(xy_Lon*Wnum))     

    xy_VelLon = - xy_GradLat_w(w_xy(xy_Stream))/Rplanet
    xy_VelLat =   xy_GradLon_w(w_xy(xy_Stream))/Rplanet

    xy_Vor  = xy_w(w_Lapla_w(w_xy(xy_Stream)))/Rplanet**2 

    xy_HDiv  = xy_w(w_DivLambdaMu_xy_xy  &
&      (xy_VelLon*cos(xy_Lat), xy_VelLat*cos(xy_Lat)))/Rplanet

    xy_Hsfc = HsfcAvr + Rplanet**2/Grav*                          &
&        (0.5D0*Omega0*(2.0D0*Omega + Omega0)*cos(xy_Lat)**2      &
&           +                                                     &
&         0.25D0*Omega0**2*cos(xy_Lat)**(2.0D0*Wnum)*             &
&            ((Wnum + 1.0D0)*cos(xy_Lat)**2                       &
&               + (2.0D0*Wnum**2 - Wnum - 2.0D0)                  &
&               - 2.0D0*Wnum**2*cos(xy_Lat)**(-2)                 &
&            )                                                    &
&           +                                                     &
&         (2.0D0*(Omega + Omega0)*Omega0                          &   
&            /(Wnum + 1.0D0)/(Wnum + 2.0D0)*cos(xy_Lat)**Wnum     & 
&            *((Wnum**2 + 2.0D0*Wnum + 2.0)                       &
&                 - (Wnum + 1.0)**2*cos(xy_Lat)**2))              &
&             * cos(xy_lon*Wnum)                                  &
&           +                                                     &
&          0.25D0*Omega0**2*cos(xy_Lat)**(2*Wnum)                 &
&             * ((Wnum + 1.0D0)*cos(xy_Lat)**2 - (Wnum + 2.0D0))  &
&             * cos(2.0D0*xy_lon*Wnum))

    xy_Topo = 0.0D0

  case default
    call MessageNotify('Error',      &
&                      'shallow_zd_init', &
&                      'invalid case number; program is forcedly terminated.')
  end select

!== 出力と後処理 ======================================================

!----------------------------------------------------------------------
!  出力
!----------------------------------------------------------------------

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

!----------------------------------------------------------------------
!  配列解放
!----------------------------------------------------------------------

  deallocate(xy_VelLon, xy_VelLat, xy_Hsfc  , &
&            xy_Vor   , xy_HDiv  , xy_Stream, &
&            xy_Topo    )

  deallocate(xy_LonOrg, xy_LatOrg, &
&            xy_VelLonOrg, xy_VelLatOrg, Xe, gamma, xy_VelLonMean)


contains

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

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

!----------------------------------------------------------------------
!  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(Tinit), 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='height', 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='height', 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='height', 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
!  初期値の出力を行う. 
!  - Adams-Bashforth スキームのため 3 ステップ分のデータを出力する.
!----------------------------------------------------------------------

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

    call HistoryPut('u_n'  ,xy_VelLon)
    call HistoryPut('v_n'  ,xy_VelLat)
    call HistoryPut('h_n'  ,xy_Hsfc)
    call HistoryPut('vor_n',xy_Vor)
    call HistoryPut('div_n',xy_HDiv)

    call HistoryPut('u_b'  ,xy_VelLon)
    call HistoryPut('v_b'  ,xy_VelLat)
    call HistoryPut('h_b'  ,xy_Hsfc)
    call HistoryPut('vor_b',xy_Vor)
    call HistoryPut('div_b',xy_HDiv)

    call HistoryPut('u_bb'  ,xy_VelLon)
    call HistoryPut('v_bb'  ,xy_VelLat)
    call HistoryPut('h_bb'  ,xy_Hsfc)
    call HistoryPut('vor_bb',xy_Vor)
    call HistoryPut('div_bb',xy_HDiv)
  end subroutine OutputRestart_gtool4

!----------------------------------------------------------------------
!  output_gtool4_close
!  入出力後処理を行う
!----------------------------------------------------------------------

  subroutine OutputRestart_gtool4_close
    call HistoryClose
  end subroutine OutputRestart_gtool4_close

end program shallow_zd_init


