program main
  use dcl
  use gms
  use gtool_history
  implicit none
  integer, parameter :: nx = 100
  integer, parameter :: nloop = 1000
  integer, parameter :: ny = 100
  real, parameter :: fdump = 1000.0
  real(8), parameter :: delta_x = 2000.0D0
  real(8), parameter :: delta_y = 2000.0D0
  real(8),parameter :: xmax  = delta_x * nx, xmin = 0.0D0
  real(8),parameter :: ymax  = delta_y * ny, ymin = 0.0D0
  integer, parameter :: margin = 1
  integer, parameter :: lbx = -margin, ubx = nx+margin
  integer, parameter :: lby = -margin, uby = ny+margin
  real(8), dimension(lbx:ubx,lby:uby) :: u0
  real(8), dimension(lbx:ubx,lby:uby) :: v0
  real(8), dimension(lbx:ubx,lby:uby) :: h0
  real(8), dimension(lbx:ubx,lby:uby) :: draw_h
  integer, dimension(3),parameter :: u_grid = (/gms_on_grid, gms_off_grid, gms_none_grid/)
  integer, dimension(3),parameter :: v_grid = (/gms_off_grid, gms_on_grid, gms_none_grid/)
  integer, dimension(3),parameter :: h_grid = (/gms_off_grid, gms_off_grid, gms_none_grid/)
  type(var_xy) :: u, u_a, u_b
  type(var_xy) :: v, v_a, v_b
  type(var_xy) :: h, h_a, h_b
  real(8), parameter :: depth = 1.0D0, grav=9.8D0, dt = 100.0D0
  real(8), parameter :: pi = 3.14159265358979D0
  integer :: i, j
  real :: time=0.0
  external initfunc
  real(8)::initfunc
	integer :: dd1, dd2

!  interface
!     real(8) function initfunc(x)
!       real(8),intent(in) :: x
!     end function initfunc
!  end interface

!  call init_graph

!---set model_parameter
  call set_grid_num_x(nx) !格子点数
  call set_grid_num_y(ny) !格子点数
  call set_margin_x(margin)     !糊代の長さ
  call set_margin_y(margin)     !糊代の長さ
  call set_real_min_x(xmin)           !x軸の最小座標値
  call set_real_min_y(ymin)           !y軸の最小座標値
  call gms_set_interval_x(delta_x)    !dx
  call gms_set_interval_y(delta_y)    !dx
  call dump_gms_modelparm 

!---initialize memory manager
  call allocate_work_area_xy(19)
!  call allocate_work_area_y(20)

!---allocate and initialize variable
  call def_var_xy(u,   u_grid)
  call def_var_xy(u_a, u_grid)
  call def_var_xy(u_b, u_grid)

  call def_var_xy(v,   v_grid)
  call def_var_xy(v_a, v_grid)
  call def_var_xy(v_b, v_grid)

  call def_var_xy(h,   h_grid)
  call def_var_xy(h_a, h_grid)
  call def_var_xy(h_b, h_grid)


  call gt4init(h, "2d.nc", fdump, "h", "displacement", "m")

 !---make initial wave

!  do i = lbx, ubx

!     h0(i) = sin( 2.0D0 * pi * ( real(i) - 0.5 ) / (real(nx)) )

!  end do

!  u0 = 0.0D0
!---------------------------------------------------------------



!  call input_var(u0, u_b)
!  call input_var(h0, h_b)

  call input_func_var_xy(initfunc, h_b) 
  u_b = 0.0D0
  v_b = 0.0D0
!  write(*,*) h_b
!---norishiro fix
  call cyclic_boundary_x_xy(u_b)
!write(*,*) u_b
  call cyclic_boundary_x_xy(v_b)
!write(*,*) v_b
  call cyclic_boundary_x_xy(h_b)
!write(*,*) h_b
  call cyclic_boundary_y_xy(u_b)
!write(*,*) u_b
  call cyclic_boundary_y_xy(v_b)
!write(*,*) v_b
  call cyclic_boundary_y_xy(h_b)
!write(*,*) h_b


!---initial integrate by Euler method


  h = h_b -   depth * d_x(u_b) * dt - depth * d_y(v_b) * dt 
  u = u_b  -   grav * d_x(h_b) * dt 
  v = v_b  -   grav * d_y(h_b) * dt
	
	dd1 = size_x(v)
	dd2 = size_y(v)

	u0 = get(u)
do i = lbx+1, ubx-1
	do j = lby+1, uby-2	
!		write(*,*) i, j	
!		write(*,*) u0(i, j)
	end do
!	write(*,*)
end do	

!	write(*,*) dd1, dd2  

!write(*,*) v


  call cyclic_boundary_x_xy(u)
!write(*,*) u
  call cyclic_boundary_x_xy(v)
!write(*,*) v
  call cyclic_boundary_x_xy(h)
!write(*,*) h
  call cyclic_boundary_y_xy(u)
!write(*,*) u
  call cyclic_boundary_y_xy(v)
!write(*,*) v
  call cyclic_boundary_y_xy(h)
!write(*,*) h

!
   do i = 0, nloop
      time = real(dt) * i
!     h_a = h_b - depth * d_x(u) * 2.0D0 * dt
     h_a = h_b - depth * d_x(u) * 2.0D0 * dt - depth * d_y(v) * 2.0D0 * dt
     u_a = u_b -  grav * d_x(h) * 2.0D0 * dt 
     v_a = v_b -  grav * d_y(h) * 2.0D0 * dt 

!write(*,*) u_a

     call cyclic_boundary_x_xy(u_a)
     call cyclic_boundary_x_xy(v_a)
     call cyclic_boundary_x_xy(h_a)
     call cyclic_boundary_y_xy(u_a)
     call cyclic_boundary_y_xy(v_a)
     call cyclic_boundary_y_xy(h_a)

 !-- time filter

     u = u + 0.1D0 * (u_a - 2.0D0 * u + u_b) 
!write(*,*) u
     v = v + 0.1D0 * (v_a - 2.0D0 * v + v_b) 
!write(*,*) v
     h = h + 0.1D0 * (h_a - 2.0D0 * h + h_b) 
!write(*,*) h

     if ( mod(time, fdump) == 0.0 ) then 
!        call draw_graph(real(position_x(h)), real(get(h)), real(xmax))
        call gt4_timeout("t", time)
        call gt4_varout("h", h)

     end if
!write(*,*) h

!  FOR NEXT STEP
     u_b = u
     v_b = v
     h_b = h
     h   = h_a
     u   = u_a
     v   = v_a

  end do

!write(*,*) size_x_xy(h),size_y_xy(h),h

!  call end_graph
  call gt4end




contains
  subroutine draw_graph(x, y, xmax)
    real :: x(:), y(:), xmax
    call DclNewFrame                ! 新しい描画領域を作成する
    call DclSetWindow( 0.,xmax, -1., 1. )
    call DclSetViewport( 0.2, 0.8, 0.2, 0.8 )
    call DclSetTransFunction 
    
    call DclDrawScaledGraph( x, y ) ! おまかせでグラフを描く
  
  end subroutine draw_graph

  subroutine end_graph
    call DclCloseGraphics    
  end subroutine end_graph
  
  subroutine init_graph
    call DclOpenGraphics()          ! 出力装置のオープン
  end subroutine init_graph

!call gt4init(h, "2d.nc", fdump=1000, "h", "displacement", "m")

  subroutine gt4init(var, fname, ndump, vname, lname, vunit)

    type(var_xy), intent(in) :: var
    character(*), intent(in):: fname, vname, lname, vunit
    real        , intent(in)::ndump
!    real :: dt=100.0D0
!    integer :: i
    
    call HistoryCreate("2d.nc", "gms_test", "gms_shallow1d", "miki",   		&
	         & (/'x', 'y','t'/), (/size_x_xy(var),size_y_xy(var),0/), 		 	&
          	 &(/"x-coordinate","y-coordinate","time   "/), (/"m","m","s"/), 	&
				 & 0.0, ndump)

    call HistoryPut('x', pos_x(var))
    call HistoryPut('y', pos_y(var))

    call HistoryAddVariable( &                                ! 変数定義
           varname=vname, dims=(/'x','y','t'/), & 
           longname=lname, units=vunit, xtype='double')

!    call HistoryPut(vname, pos_x_xy(var))	
!    call HistoryPut(vname, pos_y_xy(var))

!		do i=1,ndump*real(dt),dt
!			read(*,*) h
!			call HistoryPut(vname, pos_x_xy(var))
!	     	call HistoryPut(vname, pos_y_xy(var))
!		end do

  end subroutine gt4init

!  subroutine gt4init(var, fname, ndump, vname, lname, vunit)

!    type(var_x), intent(in) :: var
!    character(*), intent(in):: fname, vname, lname, vunit
!    real        , intent(in)::ndump
    
!    call HistoryCreate("test.nc", "gms_test", "gms_shallow1d", "masuo", &
!	         & (/'x', 't'/), (/size_x(var),0/), &
!                 (/"x-coordinate","time        "/), (/"m","s"/), 0.0, ndump)
!    call HistoryPut('x', pos_x(var))
!    call HistoryAddVariable( &                                ! 変数定義
!           varname=vname, dims=(/'x','t'/), & 
!           longname=lname, units=vunit, xtype='double')
!  end subroutine gt4init

  subroutine gt4_varout(vname, var)
    character(*), intent(in) :: vname
    type(var_xy), intent(in)  ::var

    call HistoryPut(vname, get(var) )

  end subroutine gt4_varout

  subroutine gt4_timeout(vname, time)
    character(*), intent(in) :: vname
    real, intent(in)  ::time

    call HistoryPut(vname, time )

  end subroutine gt4_timeout

  subroutine gt4end
    call HistoryClose
  end subroutine gt4end

end program main

  real(8) function initfunc(x,y)
    real(8), intent(in) :: x 
    real(8), intent(in) :: y
    real(8), parameter ::pi = 3.14159265358979D0
    initfunc = sin( 2.0D0 * pi * x / 4.0D5 )

    initfunc = 1./ exp(((x-100000.)/10000.) ** 2. + ((y-100000.)/10000.)**2.)
!write(*,*) x, y,initfunc

  end function initfunc