!----------------------------------------------------------------------
!     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 barotropic system on a sphere
!     d\zeta/dt + J(\psi,\zeta) = \nabla^2\zeta
!     \nabla^2\psi = \zeta
!
  program spbaro1

    use sppack
    use gtool_history
    implicit none

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

  !---- 変数 ----
    real(8)            :: gg_zeta(im,jm)           ! 格子点データ(渦度)
    real(8)            :: gg_psi(im,jm)            ! 格子点データ(流線)
    real(8)            :: ep_zeta((nm+1)*(nm+1))   ! スペクトル(渦度)
    real(8)            :: ep_psi((nm+1)*(nm+1))    ! スペクトル(流線)

  !---- 座標変数など ----
    real(8), parameter :: r0=1.0

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

    real(8), parameter :: zeta0 = 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(nm,im,jm)                       ! ISPACK初期化

  !------------------- 初期値設定 ----------------------
    gg_zeta = zeta0 * exp( -((gg_lon-x0)/ax)**2 -((gg_lat-y0)/ay)**2 )
    ep_zeta = ep_gg(gg_zeta)
    ep_zeta(l_nm(0,0)) = 0.0
    ep_psi = ep_laplainv_ep(ep_zeta)
    gg_psi = gg_ep(ep_psi)

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

 !------------------- 時間積分 ----------------------
    do it=1,nt
       ep_zeta = ep_zeta &
            + dt * ( - ep_jacobian_ep_ep(ep_psi,ep_zeta) &
                     - omega * ep_dlon_ep(ep_psi) &
                     + nu * ep_lapla_ep(ep_zeta) )
       ep_psi = ep_laplainv_ep(ep_zeta)

       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='spbaro1.nc', title='Barotropic flow 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='psi', dims=(/'lon','lat','t  '/), & 
           longname='stream function', units='1', xtype='double')
      call HistoryAddVariable( &                          ! 変数定義
           varname='zeta', dims=(/'lon','lat','t  '/), & 
           longname='voticity', units='1', xtype='double')

    end subroutine output_gtool4_init

    subroutine output_gtool4
      gg_psi  = gg_ep(ep_psi)
      gg_zeta = gg_ep(ep_zeta)

      write(6,*) 'it = ',it
      call HistoryPut('psi',gg_psi)
      call HistoryPut('zeta',gg_zeta)
    end subroutine output_gtool4

    subroutine output_gtool4_close
      call HistoryClose
    end subroutine output_gtool4_close

  end program spbaro1
