program main
  use dcl
  use gms
  use gtool_history
  implicit none
  integer, parameter :: nx = 100
  integer, parameter :: nloop = 1000
  real, parameter :: fdump = 1000.0
  real(8), parameter :: delta_x = 2000.0D0
  real(8),parameter :: xmax  = delta_x * nx, xmin = 0.0D0
  integer, parameter :: margin = 1
  integer, parameter :: lbx = -margin, ubx = nx+margin
  real(8), dimension(lbx:ubx) :: u0
  real(8), dimension(lbx:ubx) :: h0
  real(8), dimension(lbx:ubx) :: draw_h
  integer, dimension(3),parameter :: u_grid = (/gms_on_grid, gms_none_grid, gms_none_grid/)
  integer, dimension(3),parameter :: h_grid = (/gms_off_grid, gms_none_grid, gms_none_grid/)
  type(var_x) :: u, u_a, u_b
  type(var_x) :: 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
  real :: time=0.0
  external initfunc
  real(8)::initfunc

!  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_margin_x(margin)     !糊代の長さ
  call set_real_min_x(xmin)           !x軸の最小座標値
  call gms_set_interval_x(delta_x)    !dx
  call dump_gms_modelparm 

!---initialize memory manager
  call allocate_work_area_x(20)

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

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


  call gt4init(h, "test.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(initfunc, h_b)
  u_b = 0.0D0


!---norishiro fix
  call cyclic_boundary_x(u_b)
  call cyclic_boundary_x(h_b)


!---initial integrate by Euler method


  h = h_b -   depth * d_x(u_b) * dt 
  u = u_b  -   grav * d_x(h_b) * dt 


  call cyclic_boundary_x(u)
  call cyclic_boundary_x(h)


   do i = 0, nloop
      time = real(dt) * i

     h_a = h_b - depth * d_x(u) * 2.0D0 * dt
     u_a = u_b -  grav * d_x(h) * 2.0D0 * dt 

     call cyclic_boundary_x(u_a)
     call cyclic_boundary_x(h_a)

 !-- time filter

     u = u + 0.1D0 * (u_a - 2.0D0 * u + u_b) 
     h = h + 0.1D0 * (h_a - 2.0D0 * h + h_b) 

     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


!  FOR NEXT STEP
     u_b = u
     h_b = h
     h   = h_a
     u   = u_a

  end do

!  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

  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_x), 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)
    real(8), intent(in) :: x 
    real(8), parameter ::pi = 3.14159265358979D0
!   initfunc = sin( 2.0D0 * pi * x / 2.0D5 )
!	initfunc = sin( 2.0D0 * pi * x / 4.0D5 )

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

  end function initfunc
