!----------------------------------------------------------------------
!     Copyright (c) 2001 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
! Sample program for gtool_history/gtool4 and ISPACK   2001/12/26 S.Takehiro
! 
! Solving 2-D shallow water system on a sphere
!     Du/dt  - 2\Omega\sin\phi v  =  -g/a\cos\phi dh/d\lambda,
!     Dv/dt  + 2\Omega\sin\phi u  =  -g/a  dh/d\phi,
!     dh/dt + ( 1/\cos\phi d((H+h)u)/d\lambda 
!               + 1/\cos\phi d((H+h)v\cos\phi)/d\phi ) = 0.
!
  program spshall1

    use sppack_single
    use gtool_history
    implicit none

  !---- 空間解像度設定 ----
    integer, parameter :: im=64, jm=32             ! 格子点の設定(X,Y)
    integer, parameter :: nm=21

  !---- 変数 ----
    real(8)            :: gg_u(im,jm)              ! 格子点データ(速度経度成分)
    real(8)            :: gg_v(im,jm)              ! 格子点データ(速度緯度成分)
    real(8)            :: gg_h(im,jm)              ! 格子点データ(変位)

    real(8)            :: ep_u((nm+1)*(nm+1))      ! スペクトル(速度経度成分)
    real(8)            :: ep_v((nm+1)*(nm+1))      ! スペクトル(速度緯度成分)
    real(8)            :: ep_h((nm+1)*(nm+1))      ! スペクトル(変位)

  !---- 時間積分パラメター ----
    real(8), parameter :: dt = 1e-3                ! 時間ステップ間隔
    integer, parameter :: nt = 10000, ndisp = 500  ! 時間積分数, 表示ステップ
 
  !---- 物理パラメター ----
    real(8), parameter :: pi = 3.141592653589793   ! 円周率
    real(8), parameter :: r0 = 1.0                 ! 球の半径
    real(8), parameter :: h0 = 1.0                 ! 粘性係数
    real(8), parameter :: omega = 1.0              ! 回転角速度
    real(8), parameter :: grav = 1.0               ! 粘性係数

    real(8), parameter :: a0 = 1.0                 ! 初期変位の振幅
    real(8), parameter :: x0 = pi, y0 = 0.0        ! 初期変位の位置
    real(8), parameter :: ax = pi/10, ay = pi/10   ! 初期変位の大きさ

  !------ 作業変数 ------
    integer             :: i, j, it=0

  !---------------- 座標値の設定 ---------------------
    call spinitial_single(nm,im,jm)                ! ISPACK初期化

  !------------------- 初期値設定 ----------------------
    gg_u  = 0.0 ; gg_v = 0.0
    gg_h  = a0 * exp( -((gg_lon-x0)/ax)**2 -((gg_lat-y0)/ay)**2 )

    ep_u = ep_gg(gg_u) ; ep_v = ep_gg(gg_v) 
    ep_h = ep_gg(gg_h)

    call output_gtool4_init                        ! ヒストリー初期化
    call output_gtool4                             ! 初期値出力

  !------------------- 時間積分 ----------------------
    do it=1,nt
       ep_u = ep_u &
            + dt * ep_gg( - gg_u * gg_gradlon_ep(ep_u) / r0    &
                          - gg_v * gg_gradlat_ep(ep_u) / r0    &
                          + gg_u * gg_v * tan(gg_lat) / r0   &
                          + 2 * omega * sin(gg_lat) * gg_v   &
                          - grav * gg_gradlon_ep(ep_h)/ r0   )
       gg_u = gg_ep(ep_u)

       ep_v = ep_v &
            + dt * ep_gg( - gg_u * gg_gradlon_ep(ep_v) / r0    &
                          - gg_v * gg_gradlat_ep(ep_v) / r0    &
                          - gg_u * gg_u * tan(gg_lat) / r0   &
                          - 2 * omega * sin(gg_lat) * gg_u   &
                          - grav * gg_gradlat_ep(ep_h) / r0 )
       gg_v = gg_ep(ep_v)

       ep_h = ep_h &
            + dt * ( - ep_div_gg_gg( gg_u*(H0+gg_h), gg_v*(H0+gg_h) ) / r0 )
       gg_h = gg_ep(ep_h)

       if(mod(it,ndisp) .eq. 0)then                        ! 出力
          call output_gtool4
       endif
    enddo

    call output_gtool4_close
    stop

  contains
    subroutine output_gtool4_init
      call HistoryCreate( &                               ! ヒストリー作成
           file='spshall2.nc', title='Shallow water equation on a sphere', &
           source='Sample program of gtool_history/gtool4',       &
           institution='GFD_Dennou Club davis/spmodel project',   &
           dims=(/'lon','lat','t  '/), dimsizes=(/im,jm,0/),            &
           longnames=(/'longitude','latitude ','time     '/),     &
           units=(/'1','1','1'/),                        &
           origin=0.0, interval=real(ndisp*dt) )

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

      call HistoryAddVariable( &                          ! 変数定義
           varname='h', dims=(/'lon','lat','t  '/), & 
           longname='surface displacement ', units='1', xtype='double')
      call HistoryAddVariable( &                          ! 変数定義
           varname='u', dims=(/'lon','lat','t  '/), & 
           longname='velocity(longitude) ', units='1', xtype='double')
      call HistoryAddVariable( &                          ! 変数定義
           varname='v', dims=(/'lon','lat','t  '/), & 
           longname='velocity(latitude) ', units='1', xtype='double')
    end subroutine output_gtool4_init

    subroutine output_gtool4
      write(6,*) 'it = ',it
      call HistoryPut('u',gg_u)
      call HistoryPut('v',gg_v)
      call HistoryPut('h',gg_h)
    end subroutine output_gtool4

    subroutine output_gtool4_close
      call HistoryClose
    end subroutine output_gtool4_close

  end program spshall1
