! Sample program for gtool_history/gtool4 and ISPACK   2001/11/12 S.Takehiro
! 
! 2002/01/13 TAKAHASHI Koko
! Rossby Wave :
! Solving 2-D non-divergent (barotropic) fluid system
!     d\zeta/dt  + d\zeta/dx + d\psi/dx = q
!     \lapla psi = \zeta
!     q = exp( - (x^2 + y^2)/(0.05 * /labmda) )
!
program n2_rossbywave

  use n2pack
  use gtool_history
  implicit none

 !---- 空間解像度設定 ----
  integer, parameter :: im=64, jm=64            ! 格子点の設定(X,Y)
  integer, parameter :: km=10, lm=10            ! 切断波数の設定(X,Y)

 !---- 変数 ----
  real(8)            :: gzeta(0:jm-1,0:im-1)    ! 格子データ(渦度)
  real(8)            :: szeta0(-lm:lm,-km:km)   ! スペクトルデータ(t-dt)
  real(8)            :: szeta1(-lm:lm,-km:km)   ! スペクトルデータ(t)
  real(8)            :: szeta2(-lm:lm,-km:km)   ! スペクトルデータ(t+dt)

  real(8)            :: gpsi(0:jm-1,0:im-1)     ! 格子データ(流線)
  real(8)            :: spsi(-lm:lm,-km:km)     ! スペクトルデータ

  real(8)            :: gq(0:jm-1,0:im-1)       ! 格子データ(地形)
  real(8)            :: sq(-lm:lm,-km:km)       ! スペクトルデータ

 !---- 座標変数など ----
  real(8)            :: gx(0:jm-1,0:im-1), gy(0:jm-1,0:im-1)
  real(8), parameter :: xmin=0.0, xmax=100.0
  real(8), parameter :: ymin=0.0, ymax=100.0

 !---- 時間積分パラメター ----
  real(8), parameter :: dt=1e-5                  ! 時間ステップ間隔
  integer, parameter :: nt=1000000, ndisp=50000 ! 時間積分数, 表示ステップ

 !---- 物理パラメター ----
  real(8), parameter :: rx=0.5                     ! 地形の効果
  real(8), parameter :: ry=0.5                     ! 地形の効果

  integer :: i, j, it                              ! DO 変数

 !---------------- 座標値の設定 ---------------------
  do i=0,im-1
     do j=0,jm-1
        gx(j,i)= (xmax-xmin)/im*i
        gy(j,i)= (ymax-ymin)/jm*j
     enddo
  enddo

  call n2initial(im,jm,km,lm,xmax-xmin,ymax-ymin)    ! スペクトル初期化

 !------------------- 初期値設定 ----------------------
  gpsi = 0.0
  gq = exp( - (gx - xmax/5)**2/(rx**2) - (gy - ymax/2)**2/(ry**2) )
  ! 地形の位置

  spsi = s_g(gpsi)
  szeta1 = s_lapla_s(spsi)
  szeta0 = szeta1
  sq = s_g(gq)

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

 !------------------- 時間積分 ----------------------
  do it=1,nt                                         ! Leap frog による時間積分
      szeta2 = szeta0 &
            + 2 * dt * ( - s_dx_s(sq) &
            - s_dx_s(szeta0) - s_dx_s(spsi) )
      spsi = s_laplainv_s(szeta2)

     szeta0 = szeta1 ; szeta1 = szeta2

     if(mod(it,ndisp) .eq. 0)then                    ! 出力
        gzeta = g_s(szeta1)
         gpsi  = g_s(spsi)
        call output_gtool4
     endif
  enddo

  call output_gtool4_close                           ! ヒストリー後処理
  stop

contains

  subroutine output_gtool4_init
    call HistoryCreate( &                                  ! ヒストリー作成
           file='rossby-n2.nc', title=' model',  &
           source='Sample program of gtool_history/gtool4', &
           institution='GFD_Dennou Club davis project',     &
           dims=(/'x','y','t'/), dimsizes=(/im,jm,0/),      &
           longnames=(/'X-coordinate','Y-coordinate','time        '/),&
           units=(/'1','1','1'/),                           &
           origin=0.0, interval=real(ndisp*dt) )

    call HistoryPut('x',gx(0,0:im-1))                         ! 変数出力
    call HistoryPut('y',gy(0:jm-1,0))                         ! 変数出力

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

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

  subroutine output_gtool4_close
    call HistoryClose
  end subroutine output_gtool4_close

end program n2_rossbywave