! Sample program for gtool_history/gtool4 and ISPACK   2001/11/29 S.Takehiro
! 
! Solving 2-D barotropic system on a sphere
!     d\zeta/dt + J(\psi,\zeta) = \nabla^2\zeta
!     \nabla^2\psi = \zeta
!
  program stbaro1

    use stpack
    use gtool_history
    implicit none

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

  !---- 変数 ----
    real(8)            :: gzeta(jm,im)          ! 格子点データ(渦度)
    real(8)            :: gpsi(jm,im)           ! 格子点データ(流線)
    real(8)            :: szeta((nm+1)*(nm+1))  ! スペクトル(渦度)
    real(8)            :: spsi((nm+1)*(nm+1))   ! スペクトル(流線)

  !---- 座標変数など ----
    real(8)            :: gx(jm,im)             ! 格子点データ(経度)
    real(8)            :: gy(jm,im)             ! 格子点データ(緯度)
    real(8)            :: sq0((nm+1)*(nm+1))    ! スペクトル(惑星渦度)
    real(8), parameter :: radius=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/6, ay = pi/6     ! 初期渦の大きさ

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


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

    do j=1,jm
       gx(j,:)=xgrid()
    enddo
    do i=1,im 
       gy(:,i)=ygrid()
    enddo

  !------------------- 初期値設定 ----------------------
    sq0 = s_g(2 * omega * sin(gy))

  !------------------- 初期値設定 ----------------------
    gzeta = zeta0 * exp( -((gx-x0)/ax)**2 -((gy-y0)/ay)**2 )
    szeta = s_g(gzeta)
    szeta(l_nm(0,0)) = 0.0
    spsi  = s_laplainv_s(szeta)
    gpsi  = g_s(spsi)

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

 !------------------- 時間積分 ----------------------
    do it=1,nt
       szeta = szeta + dt * ( - s_jac_ss(spsi,sq0) + nu*s_lapla_s(szeta) )
       spsi = s_laplainv_s(szeta)

       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='stbaro1.nc', title='Barotropic flow on a sphere', &
           source='Sample program of gtool_history/gtool4',       &
           institution='GFD_Dennou Club davis/spmodel project',   &
           dims=(/'x','y','t'/), dimsizes=(/im,jm,0/),            &
           longnames=(/'longitude','latitude ','time     '/),     &
           units=(/'1','1','1'/),                        &
           origin=0.0, interval=real(ndisp*dt) )

      call HistoryPut('x',gx(1,1:im)*180/pi)              ! 変数出力
      call HistoryPut('y',gy(1:jm,1)*180/pi)              ! 変数出力

      call HistoryAddVariable( &                          ! 変数定義
           varname='psi', dims=(/'x','y','t'/), & 
           longname='stream function', units='1', xtype='double')
      call HistoryAddVariable( &                          ! 変数定義
           varname='zeta', dims=(/'x','y','t'/), & 
           longname='voticity', units='1', xtype='double')

    end subroutine output_gtool4_init

    subroutine output_gtool4
      gpsi  = g_s(spsi)
      gzeta = g_s(szeta)

      write(6,*) 'it = ',it
      call HistoryPut('psi',transpose(gpsi))
      call HistoryPut('zeta',transpose(gzeta))
    end subroutine output_gtool4

    subroutine output_gtool4_close
      call HistoryClose
    end subroutine output_gtool4_close

  end program stbaro1
