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 ) write(*,*) x, initfunc end function initfunc