!----------------------------------------------------------------------
!     Copyright (c) 2001 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
! Sample program for gtool_history/gtool4 and ISPACK   2001/12/26 S.Takehiro
! 
! Solving nonlinear 2-D shallow water system on a sphere
!     Du/dt  - uv\tan\phi/a - 2\Omega\sin\phi v  =  -g/a\cos\phi dh/d\lambda,
!     Dv/dt  + u^2\tan\phi/a + 2\Omega\sin\phi u  =  -g/a dh/d\phi,
!     dh/dt + (   1/\cos\phi d((H+h)u/d\lambda 
!               + 1/\cos\phi d((H+h)v\cos\phi)/d\phi ) = 0.
!
program shall_n

  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)            :: gg_u(0:jm-1,0:im-1)     ! 格子データ(X 速度)
  real(8)            :: gg_v(0:jm-1,0:im-1)     ! 格子データ(Y 速度)
  real(8)            :: gg_h(0:jm-1,0:im-1)     ! 格子データ(高度)

  real(8)            :: ee_u(-lm:lm,-km:km)     ! スペクトルデータ
  real(8)            :: ee_v(-lm:lm,-km:km)     ! スペクトルデータ
  real(8)            :: ee_h(-lm:lm,-km:km)     ! スペクトルデータ

  real(8)            :: gg_zeta(0:jm-1,0:im-1)  ! 格子データ(渦度)
  real(8)            :: ee_zeta(-lm:lm,-km:km)  ! スペクトルデータ

 !---- 座標変数など ----
  real(8)            :: gg_x(0:jm-1,0:im-1), gg_y(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)            :: ee_nu(-lm:lm,-km:km)    ! スペクトルデータ

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

 !---- 物理パラメター ----
  real(8), parameter :: grav=1.0                ! 重力
  real(8), parameter :: H0=0.05                 ! 平均深さ
  real(8), parameter :: f=0.0                   ! コリオリパラメター

  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 変数

 !---------------- 座標値の設定 ---------------------
  do i=0,im-1
     do j=0,jm-1
        gg_x(j,i)= (xmax-xmin)/im*i
        gg_y(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
        ee_nu(l,k)=nu*((2*pi*k/(xmax-xmin))**2+(2*pi*l/(ymax-ymin))**2)**nv
     enddo
  enddo

 !------------------- 初期値設定 ----------------------
  gg_zeta=  exp(-((gg_x-x1)**2+(gg_y-y1)**2)/(2*sigma**2)) &
       + exp(-((gg_x-x2)**2+(gg_y-y2)**2)/(2*sigma**2))
  ee_zeta =  ee_gg(gg_zeta) ; ee_zeta(0,0)=0.0 ; gg_zeta= gg_ee(ee_zeta)

  gg_h = 0.0 ; ee_h = ee_gg(gg_h)
  ee_u = -ee_dy_ee(ee_laplainv_ee(ee_gg(gg_zeta)))
  ee_v =  ee_dx_ee(ee_laplainv_ee(ee_gg(gg_zeta)))
  gg_u = gg_ee(ee_u) ; gg_v = gg_ee(ee_v)
  gg_zeta = gg_ee(ee_dx_ee(ee_v)-ee_dy_ee(ee_u))

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

 !------------------- 時間積分 ----------------------
  do it=1,nt
     ee_u = ( ee_u + dt * ( &                              ! 時間積分
                   - ee_gg(gg_u*gg_ee(ee_dx_ee(ee_gg(gg_u)))) &
                   - ee_gg(gg_v*gg_ee(ee_dy_ee(ee_gg(gg_u)))) &
                   + f * ee_v - grav * ee_dx_ee(ee_h) ) ) &
          /(1+ee_nu*dt)

     ee_v = ( ee_v + dt * ( &
                   - ee_gg(gg_u*gg_ee(ee_dx_ee(ee_gg(gg_v)))) &
                   - ee_gg(gg_v*gg_ee(ee_dy_ee(ee_gg(gg_v)))) &
                   - f * ee_u - grav * ee_dy_ee(ee_h) ) ) &
          /(1+ee_nu*dt)

     ee_h = ee_h + dt * ( & 
          H0 *( -ee_dx_ee(ee_u) - ee_dy_ee(ee_v) ) &
          - ee_dx_ee(ee_gg(gg_h*gg_u)) - ee_dy_ee(ee_gg(gg_h*gg_v)) )

     gg_u = gg_ee(ee_u) ; gg_v = gg_ee(ee_v) ; gg_h = gg_ee(ee_h)
     gg_zeta = gg_ee(ee_dx_ee(ee_v)-ee_dy_ee(ee_u))

     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='shall_n.nc', title='2D nonlinear 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',gg_x(0,0:im-1))                         ! 変数出力
    call HistoryPut('y',gg_y(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')
    call HistoryAddVariable( &                                ! 変数定義
           varname='zeta', dims=(/'x','y','t'/), & 
           longname='vorticity', units='1/s', xtype='double')
  end subroutine output_gtool4_init

  subroutine output_gtool4
    write(6,*) 'it = ',it
    call HistoryPut('u',transpose(gg_u))
    call HistoryPut('v',transpose(gg_v))
    call HistoryPut('h',transpose(gg_h))
    call HistoryPut('zeta',transpose(gg_zeta))
  end subroutine output_gtool4

  subroutine output_gtool4_close
    call HistoryClose
  end subroutine output_gtool4_close

end program shall_n
