!----------------------------------------------------------------------
!     Copyright (c) 2001 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
! Sample program for gtool_history/gtool4 and ISPACK   2001/12/26 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)            :: gg_u(0:jm,0:im-1)      ! 格子データ(速度 X 成分)
    real(8)            :: gg_v(0:jm,0:im-1)      ! 格子データ(速度 Y 成分)
    real(8)            :: gg_h(0:jm,0:im-1)      ! 格子データ(表面変位)
    real(8)            :: ec_u(-km:km,0:lm)      ! スペクトル(速度 X 成分)
    real(8)            :: es_v(-km:km,lm)        ! スペクトル(速度 Y 成分)
    real(8)            :: ec_h(-km:km,0:lm)      ! スペクトル(表面変位)

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

 !---- 座標変数など ----
    real(8)            :: gg_x(0:jm,0:im-1)      ! 格子データ(X)
    real(8)            :: gg_y(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
          gg_x(j,i)= (xmax-xmin)/im*i + xmin
          gg_y(j,i)= (ymax-ymin)/jm*j + ymin
       enddo
    enddo

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

 !-------- 初期値設定 -----------
    ! 初期値なにもなし
    gg_h = 0.0
    gg_u = 0.0
    gg_v = 0.0

    ! 初期値を式で与える場合
!    gg_h = sin(4*pi/(xmax-xmin)*gg_x)*cos(pi/(ymax-ymin)*gg_y)
!    gg_u = sin(4*pi/(xmax-xmin)*gg_x)*cos(pi/(ymax-ymin)*gg_y)
!    gg_h = exp( -(gg_x**2+gg_y**2)/a0**2 )

    ! リスタートファイル読み込みの場合
!    call LoadVariable('shallow_eqbeta_init.nc@uvel,t=50', gg_data, im*(jm+1))
!    gg_u=transpose(gg_data) ; ec_u = ec_gg(gg_u)
!    call LoadVariable('shallow_eqbeta_init.nc@vvel,t=50', gg_data, im*(jm+1))
!    gg_v=transpose(gg_data) ; es_v = es_gg(gg_v)
!    call LoadVariable('shallow_eqbeta_init.nc@height,t=50', gg_data, im*(jm+1))
!    gg_h=transpose(gg_data) ; ec_h = ec_gg(gg_h)

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

    ec_qh=ec_gg(gg_qh)

    call output_gtool4_init
    call output_gtool4

 !-------- 時間積分 -----------
    do it=1,nt
       ec_u = ec_u &
              + dt* (   ec_gg( beta * gg_y * (gg_es(es_v))) &
                      - grav * ec_dx_ec(ec_h) )
       es_v = es_v &
              + dt* (  -es_gg( beta * gg_y * (gg_ec(ec_u))) &
                      - grav * es_dy_ec(ec_h) )
       ec_h = ec_h &
              + dt*( -H0*(ec_dx_ec(ec_u) + ec_dy_es(es_v)) -ec_qh )

       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',gg_x(1,0:im-1))                          ! 変数出力
      call HistoryPut('y',gg_y(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/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
      gg_u   = gg_ec(ec_u)
      gg_v   = gg_es(es_v)
      gg_h   = gg_ec(ec_h)

      call HistoryPut('t',real(ndisp*dt))
      call HistoryPut('uvel',transpose(gg_u))
      call HistoryPut('vvel',transpose(gg_v))
      call HistoryPut('height',transpose(gg_h))
      call HistoryPut('wvel', &
           transpose(gg_ec(-H0*(ec_dx_ec(ec_u)+ec_dy_es(es_v))-ec_qh)))
      call HistoryPut('hdiv', &
           transpose(gg_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

    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
