!----------------------------------------------------------------------
!  Copyright (c) 2002--2017 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
! Sample program for SPML(based on gt4f90io and ISPACK)
!
!    2002/08/19 S.Takehiro
!    2004/01/26 M.Odaka
!    2004/06/10 M.Odaka
!    2005/03/16 S.Takehiro
!    2005/10/12 H.Taniguchi -> S.Takehiro
!    2017/05/26 S.Takehiro
! 
! Solving equatorial beta plane shallow water equations (Gill 1980)
!     du/dt - beta y v = -g dh/dx - ru
!     dv/dt + beta y u = -g dh/dy - rv
!     dh/dt + H0( du/dx + dv/dy ) = -qh - nh
!     du/dy = v = dh/dy = 0 at y=ymin,ymax
!
  program equatorial_beta_shallow_water
    use esc_module
    use gtool_history
    implicit none

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

 !---- 変数 ----
    real(8)            :: yx_U(0:jm,0:im-1)      ! 格子データ(速度 X 成分)
    real(8)            :: yx_V(0:jm,0:im-1)      ! 格子データ(速度 Y 成分)
    real(8)            :: yx_H(0:jm,0:im-1)      ! 格子データ(表面変位)
    real(8)            :: ec_U(-km:km,0:lm)      ! スペクトル(速度 X 成分)
    real(8)            :: ec_U0(-km:km,0:lm)      ! スペクトル(速度 X 成分)
    real(8)            :: es_V(-km:km,lm)        ! スペクトル(速度 Y 成分)
    real(8)            :: ec_H(-km:km,0:lm)      ! スペクトル(表面変位)

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

 !---- 座標変数など ----
    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 :: Rd = 0.0              ! レイリー摩擦係数
    real(8), parameter :: Nc = 0.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=0

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

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

 !-------- 初期値設定 -----------
    ! 初期値なにもなし
    yx_H = 0.0
    yx_U = 0.0
    yx_V = 0.0

!!$    ! 初期値を式で与える場合
!!$    yx_H = sin(4*pi/(xmax-xmin)*yx_X)*cos(pi/(ymax-ymin)*yx_Y)
!!$    yx_U = sin(4*pi/(xmax-xmin)*yx_X)*cos(pi/(ymax-ymin)*yx_Y)
!!$    yx_H = exp( -(yx_X**2+yx_Y**2)/A0**2 )

!!$    ! リスタートファイル読み込みの場合
!!$    call HistoryGet('shallow_eqbeta_init.nc', 'uvel', xy_data, 200.0 )
!!$    yx_U = transpose(xy_data)
!!$    call HistoryGet('shallow_eqbeta_init.nc', 'vvel', xy_data, 200.0 )
!!$    yx_V = transpose(xy_data)    
!!$    call HistoryGet('shallow_eqbeta_init.nc', 'height', xy_data, 200.0 )
!!$    yx_H = transpose(xy_data)
    
    ec_U = ec_yx(yx_U)
    es_V = es_yx(yx_V)
    ec_H = ec_yx(yx_H)

 !----------------------- 熱強制分布 --------------------------
    do i=0,im-1
       do j=0,jm
          if ( abs(yx_X(j,i)) .le. Xrange ) then
             yx_QH(j,i) = Q0 * exp( -yx_Y(j,i)**2/A0 )&
                             * cos(pi*yx_X(j,i)/(2*Xrange))
          else
             yx_QH(j,i) = 0.0
          endif
       enddo
    enddo

    ec_QH=ec_yx(yx_QH)

    call output_gtool4_init
    call output_gtool4

 !-------- 時間積分 -----------
    do it=1,nt
       ec_U = ec_U &
              + dt* (   ec_yx( Beta * yx_Y * (yx_es(es_V))) &
                      - Grav * ec_dx_ec(ec_H)               &
                      - Rd * ec_U )
       es_V = es_V &
              + dt* (  -es_yx( Beta * yx_Y * (yx_ec(ec_U))) &
                      - Grav * es_dy_ec(ec_H)               &
                      - Rd * es_V )
       ec_H = ec_H &
              + dt*( -H0*(ec_dx_ec(ec_U) + ec_dy_es(es_V)) -ec_QH & 
                      - Nc * ec_H )

       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_eqbeta1.nc', &
           title='Equatorial beta-plane shallow water model ', &
           source='A sample program of ISPACK and gtool_history/gtool4', &
           institution='GFD_Dennou Club spmodel 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',x_X)                                   ! 変数出力
      call HistoryAddattr('x','topology','circular')             ! 周期属性
      call HistoryAddattr('x','modulo',xmax-xmin)                ! 周期属性
      call HistoryPut('y',y_Y)                                   ! 変数出力

      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/s', xtype='double')
      call HistoryAddVariable( &                                 ! 変数定義
           varname='hdiv', dims=(/'x','y','t'/), & 
           longname='Horizontal divergence', units='1/s', xtype='double')
    end subroutine output_gtool4_init

    subroutine output_gtool4
      write(6,*) 'it = ',it
      yx_U   = yx_ec(ec_U)
      yx_V   = yx_es(es_V)
      yx_H   = yx_ec(ec_H)

      call HistoryPut('t',real(it*dt))
      call HistoryPut('uvel',transpose(yx_U))
      call HistoryPut('vvel',transpose(yx_V))
      call HistoryPut('height',transpose(yx_H))
      call HistoryPut('wvel', &
           transpose(yx_ec(-H0*(ec_dx_ec(ec_U)+ec_dy_es(es_V))-ec_QH)))
      call HistoryPut('hdiv', &
           transpose(yx_ec(ec_dx_ec(ec_U)+ec_dy_es(es_V))))
    end subroutine output_gtool4

    subroutine output_gtool4_close
      call HistoryClose
    end subroutine output_gtool4_close

  end program equatorial_beta_shallow_water
