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