! Sample program for gtool_history/gtool4 and ISPACK   2001/11/12 S.Takehiro
! 
! Solving 2-D shallow water fluid system
!     du/dt  - fv  =  -g dh/dx,
!     dv/dt  + fu  =  -g dh/dy,
!     dh/dt + H( du/dx + dv/dy ) = 0.
!
program shall_l

  use n2pack
  use gtool_history
  implicit none

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

 !---- 変数 ----
  real(8)            :: gu(0:jm-1,0:im-1)       ! 格子データ(X 速度)
  real(8)            :: gv(0:jm-1,0:im-1)       ! 格子データ(Y 速度)
  real(8)            :: gh(0:jm-1,0:im-1)       ! 格子データ(高度)

  real(8)            :: su(-lm:lm,-km:km)       ! スペクトルデータ
  real(8)            :: sv(-lm:lm,-km:km)       ! スペクトルデータ
  real(8)            :: sh(-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=1d7
  real(8), parameter :: ymin=0.0, ymax=1d7

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

 !---- 物理パラメター ----
  real(8), parameter :: grav=9.8                ! 重力
  real(8), parameter :: H0=4000.0               ! 平均深さ
  real(8), parameter :: f=1e-4                  ! コリオリパラメター

  real(8), parameter :: x1=(xmin+xmax)/2.0      ! 初期分布 X 座標
  real(8), parameter :: y1=(xmin+xmax)/2.0      ! 初期分布 Y 座標
  real(8), parameter :: sigma=1.0d6             ! 初期分布の大きさ

  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)    ! スペクトル初期化

 !------------------- 初期値設定 ----------------------
  gh = exp(-((gx-x1)**2+(gy-y1)**2)/(2*sigma**2))
  gu = 0.0
  gv = 0.0

  su = s_g(gu) ; sv = s_g(gv) ; sh = s_g(gh)

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

 !------------------- 時間積分 ----------------------
  do it=1,nt
     su = su + dt * (   f * sv - grav * s_dx_s(sh) )      ! 時間積分
     sv = sv + dt * ( - f * su - grav * s_dy_s(sh) )
     sh = sh + dt * H0 * ( -s_dx_s(su) - s_dy_s(sv) ) 

     if(mod(it,ndisp) .eq. 0)then                    ! 出力
        gu = g_s(su) ; gv = g_s(sv) ; gh = g_s(sh)
        call output_gtool4
     endif
  enddo

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

contains

  subroutine output_gtool4_init
    call HistoryCreate( &                                  ! ヒストリー作成
           file='shall_l.nc', title='2D shallow water 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=(/'m','m','s'/),                           &
           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='u', dims=(/'x','y','t'/), & 
           longname='U-velocity', units='m/s', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='v', dims=(/'x','y','t'/), & 
           longname='V-velocity', units='m/s', xtype='double')
    call HistoryAddVariable( &                                ! 変数定義
           varname='h', dims=(/'x','y','t'/), & 
           longname='height', units='m', xtype='double')
  end subroutine output_gtool4_init

  subroutine output_gtool4
    write(6,*) 'it = ',it
    call HistoryPut('u',transpose(gu))
    call HistoryPut('v',transpose(gv))
    call HistoryPut('h',transpose(gh))
  end subroutine output_gtool4

  subroutine output_gtool4_close
    call HistoryClose
  end subroutine output_gtool4_close

end program shall_l
