! 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

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

  double precision, dimension(0:jm-1,0:im-1) :: psi, zeta   ! 格子データ
  double precision, dimension(-lm:lm,-km:km) :: spsi,szeta  ! スペクトルデータ
  double precision, dimension(0:jm-1,0:im-1) :: x, y        ! 座標変数

  double precision, parameter :: xmin=0.0, xmax=1.0  ! 計算範囲(X)
  double precision, parameter :: ymin=0.0, ymax=1.0  ! 計算範囲(Y)

  double precision, parameter :: beta=1.0, nu=1e-8   ! パラメター
  double precision, parameter :: dt=1e-1             ! 時間積分刻み
  integer, parameter          :: nt=1000, ndisp=10   ! 時間ステップ数

  double precision, parameter  :: x1=(xmax-xmin)*0.5    ! 初期渦用パラメター
  double precision, parameter  :: y1=(ymax-ymin)*0.5    ! 初期渦用パラメター
  double precision, parameter  :: sigma=0.05, zeta0=0.2 ! 初期渦用パラメター

  integer :: i,j, k,l

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

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

 !------------------- 初期値設定 ----------------------
  zeta=  zeta0*exp(-((x-x1)**2+(y-y1)**2)/(2*sigma**2))   ! 初期渦分布
  szeta = s_g(zeta)                           ! スペクトルデータへ変換
  szeta(0,0)=0.0                              ! 渦度の平均値は 0. 
  psi=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) + nu*s_lapla_s(szeta)  )

     spsi = s_laplainv_s(szeta)

     if(mod(it,ndisp) .eq. 0)then             ! 出力処理
        psi = g_s(spsi)
        zeta = 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',x(1,0:im-1))                          ! 変数出力
    call HistoryPut('y',y(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(psi))
    call HistoryPut('zeta',transpose(zeta))
  end subroutine output_gtool4

  subroutine output_gtool4_close
    call HistoryClose
  end subroutine output_gtool4_close

end program baro

