!----------------------------------------------------------------------
!  Copyright (C) 2001--2010 SPMODEL Development Group. All rights reserved.
!----------------------------------------------------------------------
! Sample program for gt4_history/gt4f90io and ISPACK   2002/08/21 S.Takehiro
!                                                      2004/01/26 M.Odaka
!                                                      2010/03/09 S.Takehiro
!
! Solving a non-linear 2-D shallow water system
!     with an isolated mountain
!     du/dt + u du/dx + v du/dy
!           - 2 omega v = -g dh/dx - \nu\lapla^4 u,
!     dv/dt + u dv/dx + v dv/dy
!           + 2 omega u = -gdh/dy - \nu\lapla^4 v,
!     dh/dt + ( d( (h-ht)u )/dx
!             + d( (h-ht)v )/dy ) = - \nu\lapla^4 h.
!
program spshallow_xy_nonlinear

  use ee_module
  use gt4_history
  implicit none

  !---- parameters for horizontal resolution ----
!!$  integer, parameter :: im=256            ! number of grid (X,Y)
!!$  integer, parameter :: im=128            ! number of grid (X,Y)
  integer, parameter :: im=64             ! number of grid (X,Y)
  integer, parameter :: jm=im             ! number of grid (X,Y)
  integer, parameter :: nm=(im-1)/3

  real(8), parameter :: XLength = 10000.0d3 ! unit is m
  real(8), parameter :: YLength = XLength

  !---- variables ----
  real(8)            :: yx_U(0:jm-1,0:im-1)          ! zonal flow (grid)
  real(8)            :: yx_V(0:jm-1,0:im-1)          ! meridional flow (grid)
  real(8)            :: yx_H(0:jm-1,0:im-1)          ! water surface height (grid)

  real(8)            :: ee_U(-nm:nm,-nm:nm)       ! zonal flow (spectrum)
  real(8)            :: ee_V(-nm:nm,-nm:nm)       ! meridional flow (spectrum)
  real(8)            :: ee_H(-nm:nm,-nm:nm)       ! water surface height (spectrum)

  real(8)            :: yx_H0(0:jm-1,0:im-1)         ! basic surface height (grid)
  real(8)            :: ee_H0(-nm:nm,-nm:nm)      ! basic surface height (spectrum)

  real(8)            :: yx_Zeta(0:jm-1,0:im-1)       ! vorticity (grid)
  real(8)            :: yx_Htopo(0:jm-1,0:im-1)      ! topography (grid)

  !---- physical parameters ----
  real(8), parameter :: pi = 3.141592653589793   ! 
  real(8), parameter :: H0 = 10D3                ! basic water depth
!!$  real(8), parameter :: Omega = 0.0d0            ! rotation rate
!!$  real(8), parameter :: Omega = 2*pi/86400.0d0   ! rotation rate
  real(8), parameter :: Omega = 2*pi/8640.0d0    ! rotation rate
  real(8), parameter :: Grav = 9.8               ! gravitational constant

  real(8), parameter :: Lon0 = pi, Lat0=pi/6     ! central location of mountain
  real(8), parameter :: W0 = 2000.0d3            ! inverse of initial bump width
  real(8), parameter :: Htopo0= 2500             ! mountain height

  integer, parameter :: ndiff=4                  ! order of hyper-diffusion
  real(8), parameter :: Nu=2.338D16              ! hyper-diffusion coefficient

  !---- parameters for time integration ----
!!$  real(8), parameter :: dt = 100.0/(im/64)          ! time increment
!!$  integer, parameter :: nt = 1500*im/64             ! integration step
!!$  integer, parameter :: ndisp = 5*im/64             ! display step

  real(8), parameter :: dt = XLength / im / sqrt( Grav * H0 ) * 0.1 ! time increment
  integer, parameter :: nt = XLength / sqrt( Grav * H0 ) * 0.1      ! integration step
  integer, parameter :: ndisp = XLength / sqrt( Grav * H0 ) / 1000  ! display step

  !------ work variable ------
  integer            :: i, j, it=0
  real(8)            :: Delta

  write( 6, * ) dt, nt, ndisp
  
  !---------------- setting coordinates ---------------------
  call ee_Initial(im,jm,nm,nm,0.0D0,XLength,0.0D0,YLength)

  !------------------- setting initial condition ----------------------
  yx_U  = 0.0d0
  yx_V  = 0.0d0

  do j=0,jm-1
    do i=0,im-1
      delta = sqrt( (XLength/2-x_X(j))**2 + (YLength/2.0-y_Y(i))**2 )
      if ( abs(delta) <= W0 ) then 
        yx_H(j,i)= H0 + 1.0d0 * (1+cos(pi/2*Delta/W0))/2
      else
        yx_H(j,i)= H0
      endif
    enddo
  enddo


  ee_U = ee_yx(yx_U)
  ee_V = ee_yx(yx_V)
  ee_H = ee_yx(yx_H)

  !------------------- setting topography ----------------------

!!$  do j=0,jm-1
!!$    do i=0,im-1
!!$      delta = sqrt( (XLength/2-x_X(j))**2 + (YLength-y_Y(i))**2 )
!!$      if ( abs(delta) <= W0 ) then 
!!$        yx_Htopo(j,i) = htopo0 * (1+cos(pi/2*Delta/W0))/2
!!$      else
!!$        yx_Htopo(j,i) = 0
!!$      endif
!!$    enddo
!!$  enddo
  yx_Htopo = 0.0d0

  call output_gtool4_init                        ! initialization of output
  call output_gtool4                             ! output initial condition

  !------------------- time integration ----------------------
  do it=1,nt
    ee_U = ( ee_U &
      &     + dt * ee_yx( - yx_U * yx_ee( ee_Dx_ee(ee_U) )     &
      &                   - yx_V * yx_ee( ee_Dy_ee(ee_U) )     &
      &                   + 2 * Omega * yx_V                   &
      &                   - Grav * yx_ee( ee_Dx_ee(ee_H) )   ) &
      &   )

    yx_U = yx_ee(ee_U)

    ee_V = (   ee_V &
      &     + dt * ee_yx( - yx_U * yx_ee( ee_Dx_ee(ee_V) )    &
      &                   - yx_V * yx_ee( ee_Dy_ee(ee_V) )    &
      &                   - 2 * Omega * yx_U                  &
      &                   - Grav * yx_ee( ee_Dy_ee(ee_H) ) )  &
      &   )

    yx_V = yx_ee(ee_V)

    ee_H = (   ee_H                                                &
      &     + dt * ( - ee_Dx_ee( ee_yx( yx_U*(yx_H-yx_Htopo) ) )   &
      &              - ee_Dy_ee( ee_yx( yx_V*(yx_H-yx_Htopo) ) ) ) &
      &   )

    yx_H = yx_ee(ee_H)

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

  call output_gtool4_close
  stop

contains
  subroutine output_gtool4_init
    call HistoryCreate( &                               ! create output file
      file='spshallow_xy.nc',                         &
      title='Shallow water equation on the Cartesian coordinate', &
      source='Sample program of gtool_History/gtool4',       &
      institution='GFD_Dennou Club davis/spmodel project',   &
      dims=(/'x','y','t'/), dimsizes=(/im,jm,0/),            &
      longnames=(/'x   ','y   ','time'/),                    &
      units=(/'m   ','m   ','sec.'/),                        &
      origin=0.0, interval=real(ndisp*dt) )

    call HistoryPut('x',x_X)                         ! output longitude
    call HistoryAddattr('x','topology','circular')   ! add attribute
    call HistoryPut('y',y_Y)                         ! output latitude
    call HistoryAddattr('y','topology','circular')   ! add attribute

    call HistoryAddVariable( &                          ! define variable
      varname='h', dims=(/'y','x','t'/), & 
      longname='surface displacement ', units='m', xtype='double')
    call HistoryAddVariable( &                          ! define variable
      varname='u', dims=(/'y','x','t'/), & 
      longname='velocity(longitude) ', units='m/s', xtype='double')
    call HistoryAddVariable( &                          ! define variable
      varname='v', dims=(/'y','x','t'/), & 
      longname='velocity(latitude) ', units='m/s', xtype='double')

    call HistoryAddVariable( &                          ! define variable
      varname='zeta', dims=(/'y','x','t'/), & 
      longname='vorticity', units='1/s', xtype='double')

    call HistoryAddVariable( &                          ! define variable
      varname='htopo', dims=(/'y','x'/), & 
      longname='orography', units='m', xtype='double')
    call HistoryPut('htopo',yx_Htopo)
  end subroutine output_gtool4_init

  subroutine output_gtool4
    write(6,*) 'it = ',it
    call HistoryPut('t',real(it*dt))
    call HistoryPut('u',yx_U)
    call HistoryPut('v',yx_V)
    call HistoryPut('h',yx_H)
    yx_Zeta = yx_ee(ee_Dx_ee(ee_yx(yx_V)) - ee_Dy_ee(ee_yx(yx_U)))
    call HistoryPut('zeta',yx_Zeta)
  end subroutine output_gtool4

  subroutine output_gtool4_close
    call HistoryClose
  end subroutine output_gtool4_close

end program spshallow_xy_nonlinear
