! Sample program for gtool_history/gtool4 and ISPACK   2001/11/29 S.Takehiro
! 
! Solving equatorial beta plane shallow water equations (Gill 1980)
!     du/dt - beta y v = -g dh/dx
!     dv/dt + beta y u = -g dh/dy
!     dh/dt + H0( du/dx + dv/dy ) = -qh
!     du/dy = v = dh/dy = 0 at y=ymin,ymax
!
  program equatorial_beta_shallow_water

    use c2pack
    use gtool_history
    use gtool
    implicit none

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

 !---- 変数 ----
    real(8)            :: guvel(0:jm,0:im-1)     ! 格子データ(速度 X 成分)
    real(8)            :: gvvel(0:jm,0:im-1)     ! 格子データ(速度 Y 成分)
    real(8)            :: gheight(0:jm,0:im-1)   ! 格子データ(表面変位)
    real(8)            :: cuvel(-km:km,0:lm)     ! スペクトル(速度 X 成分)
    real(8)            :: svvel(-km:km,lm)       ! スペクトル(速度 Y 成分)
    real(8)            :: cheight(-km:km,0:lm)   ! スペクトル(表面変位)

    real(8)            :: gqh(0:jm,0:im-1)       ! 格子データ(熱源分布)
    real(8)            :: cqh(-km:km,0:lm)       ! スペクトル(熱源分布)
    !real(8)            :: gdata(0:im-1,0:jm)                                 

 !---- 座標変数など ----
    real(8)            :: gx(0:jm,0:im-1)          ! 格子データ(X)
    real(8)            :: gy(0:jm,0:im-1)          ! 格子データ(Y)
    real(8), parameter :: xmin=-200.0, xmax=200.0  ! 領域範囲(X 方向)
    real(8), parameter :: ymin=-5.0, ymax=5.0      ! 領域範囲(Y 方向) 

 !---- 物理パラメター ----
    real(8), parameter :: grav=1.0              ! 重力加速度
    real(8), parameter :: beta=1.0              ! 赤道βパラメタ
    real(8), parameter :: H0 = 1.0              ! 平均の深さ
    real(8), parameter :: q0=1.0                ! 熱源分布の振幅
    real(8), parameter :: a0=0.2                ! 熱源分布の Y 方向の幅
    real(8), parameter :: xrange=20.0           ! 熱源分布の X 方向の幅

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

    integer :: i, j, it

    real(8) :: pi
    pi = atan(1.0D0)*4.0D0

    do i=0,im-1
       do j=0,jm
          gx(j,i)= (xmax-xmin)/im*i + xmin
          gy(j,i)= (ymax-ymin)/jm*j + ymin
       enddo
    enddo

    call c2initial(im,jm,km,lm,(xmax-xmin),(ymax-ymin))

 !-------- 初期値設定 -----------
    ! 初期値なにもなし
    gheight = 0.0
    guvel = 0.0
    gvvel = 0.0

    ! 初期値を式で与える場合
!    gheight = sin(4*pi/(xmax-xmin)*gx)*cos(pi/(ymax-ymin)*gy)
!    guvel = sin(4*pi/(xmax-xmin)*gx)*cos(pi/(ymax-ymin)*gy)
!    gheight = exp( -(gx**2+gy**2)/a0**2 )

    ! リスタートファイル読み込みの場合
!    call LoadVariable('shallow_eqbeta_init.nc@uvel,t=50', gdata, im*(jm+1))
!    guvel=transpose(gdata) ; cuvel = c_g(guvel)
!    call LoadVariable('shallow_eqbeta_init.nc@vvel,t=50', gdata, im*(jm+1))
!    gvvel=transpose(gdata) ; svvel = s_g(gvvel)
!    call LoadVariable('shallow_eqbeta_init.nc@height,t=50', gdata, im*(jm+1))
!    gheight=transpose(gdata) ; cheight = c_g(gheight)

 !----------------------- 熱強制分布 --------------------------
    do i=0,im-1
       do j=0,jm
          if ( abs(gx(j,i)) .le. xrange ) then
             gqh(j,i) = q0*exp( -gy(j,i)**2/a0 )*cos(pi*gx(j,i)/(2*xrange))
          else
             gqh(j,i) = 0.0
          endif
       enddo
    enddo

    cqh=c_g(gqh)

    call output_gtool4_init
    call output_gtool4

 !-------- 時間積分 -----------
    do it=1,nt
       cuvel = cuvel &
              + dt* (   c_g(beta*gy*(g_s(svvel))) &
                      - grav*c_dx_c(cheight) )
       svvel = svvel &
              + dt* (  -s_g(beta*gy*(g_c(cuvel))) &
                      - grav*s_dy_c(cheight) )
       cheight = cheight &
              + dt*( -H0*(c_dx_c(cuvel) + c_dy_s(svvel)) -cqh )

       if(mod(it,ndisp) .eq. 0)then                      ! 出力
          call output_gtool4
       endif
    enddo

    call output_gtool4_close
    stop

  contains

   !----------------------- gtool4 設定 ----------------------------
    subroutine output_gtool4_init
      call HistoryCreate( &                             ! ヒストリー作成
           file='shallow_eqbeta.nc', title='Equatorial beta-plane shallow water model ', &
           source='A sample program of ISPACK and gtool_history/gtool4', &
           institution='GFD_Dennou Club davis project',     &
           dims=(/'x','y','t'/), dimsizes=(/im,jm+1,0/),    &
           longnames=(/'X-coordinate','Y-coordinate','time        '/),&
           units=(/'m','m','s'/),                           &
           origin=0.0, interval=real(ndisp*dt) )

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

      call HistoryAddVariable( &                                 ! 変数定義
           varname='uvel', dims=(/'x','y','t'/), & 
           longname='X-velocity', units='m/s', xtype='double')
      call HistoryAddVariable( &                                 ! 変数定義
           varname='vvel', dims=(/'x','y','t'/), & 
           longname='Y-velocity', units='m/s', xtype='double')
      call HistoryAddVariable( &                                 ! 変数定義
           varname='height', dims=(/'x','y','t'/), & 
           longname='height', units='m', xtype='double')
      call HistoryAddVariable( &                                 ! 変数定義
           varname='wvel', dims=(/'x','y','t'/), & 
           longname='Z-velocity', units='m', xtype='double')
    end subroutine output_gtool4_init

    subroutine output_gtool4
      write(6,*) 'it = ',it
      guvel   = g_c(cuvel)
      gvvel   = g_s(svvel)
      gheight = g_c(cheight)

      call HistoryPut('uvel',transpose(guvel))
      call HistoryPut('vvel',transpose(gvvel))
      call HistoryPut('height',transpose(gheight))
      call HistoryPut('wvel', &
           transpose(g_c(-c_dx_c(cuvel)-c_dy_s(svvel)-cqh)))
    end subroutine output_gtool4

    subroutine output_gtool4_close
      call HistoryClose
    end subroutine output_gtool4_close

    subroutine LoadVariable(var_uri, array, length)
        ! 変数を指示する URI
        character(len = *), intent(in):: var_uri
        ! 配列とその長さ
        integer, intent(in):: length
        real(8), intent(out):: array(length)
        type(GT_VARIABLE):: var
    continue
        call open(var, var_uri)
        call get(var, array, length)
        call close(var)
      end subroutine LoadVariable

  end program equatorial_beta_shallow_water
