! Sample program for gtool_history/gtool4 and ISPACK   2001/09/23 S.Takehiro
! 
! Solving 2-D non-divergent (barotropic) fluid system
!     d\zeta/dt + J(\psi,\zeta) + beta d\psi/dx = \nabla^2\zeta,
!     \nabla\psi = \zeta
!
program baro

  use n2pack
  use gtool_history
  implicit none

!  integer,parameter      :: im=32, jm=32     ! 格子点の設定(X,Y)
!  integer,parameter      :: km=10, lm=10     ! 切断波数の設定(X,Y)
!  integer,parameter      :: km=15, lm=15     ! 切断波数の設定(X,Y)
  integer,parameter      :: im=128, jm=128   ! 格子点の設定(X,Y)
  integer,parameter      :: km=63,  lm=63    ! 切断波数の設定(X,Y)

 !---- 変数 ----
  real(8) :: gpsi(0:jm-1,0:im-1)              ! 格子データ(流線)
  real(8) :: gzeta(0:jm-1,0:im-1)             ! 格子データ(渦度)

  real(8) :: spsi(-lm:lm,-km:km)             ! スペクトルデータ
  real(8) :: szeta(-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=1.0
  real(8), parameter :: ymin=0.0, ymax=1.0

 !---- 粘性設定 ----
  real(8), parameter :: nu=1e-17             ! 粘性係数
  integer, parameter :: nv=5                 ! 粘性の階数
  real(8)            :: snu(-lm:lm,-km:km)   ! スペクトルデータ

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

 !---- 物理パラメター ----
  real(8), parameter  :: x1=(xmax-xmin)*0.4  ! 初期渦1位置 X 座標
  real(8), parameter  :: y1=(ymax-ymin)*0.4  ! 初期渦1位置 Y 座標
  real(8), parameter  :: x2=(xmax-xmin)*0.6  ! 初期渦2位置 X 座標
  real(8), parameter  :: y2=(ymax-ymin)*0.6  ! 初期渦2位置 Y 座標
  real(8), parameter  :: sigma=0.05          ! 初期渦の大きさ

  real(8), parameter :: pi=3.141592653589793
  integer :: i,j, k,l, it

  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)

 !------------------- 粘性設定 ----------------------
  do k=-km,km
     do l=-lm,lm
        snu(l,k)=nu*((2*pi*k/(xmax-xmin))**2+(2*pi*l/(ymax-ymin))**2)**nv
     enddo
  enddo

 !------------------- 初期値設定 ----------------------
  gzeta=  exp(-((gx-x1)**2+(gy-y1)**2)/(2*sigma**2)) &
       + exp(-((gx-x2)**2+(gy-y2)**2)/(2*sigma**2))
  szeta = s_g(gzeta)
  szeta(0,0)=0.0
  spsi=s_laplainv_s(szeta)

  call output_gtool4_init                            ! ヒストリー初期化

  do it=1,nt
    ! szeta = szeta + &
    !        dt*(  - s_jac_ss(spsi,szeta) &
    !             - beta*s_dx_s(spsi) - snu*szeta  )

     szeta = (szeta + dt*(-s_jac_ss(spsi,szeta)))/(1+snu*dt)

     spsi = s_laplainv_s(szeta)

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

  call output_gtool4_close
  stop

contains
  subroutine output_gtool4_init
    call HistoryCreate( &                             ! ヒストリー作成
           file='baro.nc', title='2D barotrpic 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(1,0:im-1))                         ! 変数出力
    call HistoryPut('y',gy(0:jm-1,1))                         ! 変数出力

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

  subroutine output_gtool4
    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 baro

