character*4 cdnm
parameter ( cdnm = 'step' )
logical ltrace
parameter ( ltrace = .false. )
integer i, Nx
real*8 set_floor
real*8 delta_t, delta_tb2
logical independent_residual
parameter ( independent_residual = .true. )
integer ret, calc_prim
c============================================================
c============================================================
c-------------------------------------------------------------
c Use a better variable for the number of grid points.
c Notice the RNPL variable is g1_nx.
c-------------------------------------------------------------
nx = g1_nx
delta_t = lambda * (x(2)-x(1))
delta_tb2 = 0.5d0 * delta_t
if ( ltrace ) then
write(0,*) cdnm,': Nx=',Nx
write(0,*) cdnm,': Ng=',Ng
write(0,*) cdnm,': dt=',delta_t
endif
c-------------------------------------------------------------
c First step. Computes quantities at half-advanced time step.
c On return from this section, tau_np1, etc. contain an
c estimate of the dynamical variables at the half step.
c-------------------------------------------------------------
c-------------------------------------------------------------
c Tracing output to unit 80, which, by default, will go to
c an ASCII file called 'fort.80'
c-------------------------------------------------------------
write(80,*) '***************************************'
write(80,*) '*** 1. SPATIAL COORDINATES *******'
write(80,*) '***************************************'
call dvdump(x,Nx,'r',80)
write(80,*)
write(80,*) '***************************************'
write(80,*) '*** 1. BEFORE CALL TO CALC_FLUX *******'
write(80,*) '***************************************'
call dvdump(rho_n,Nx,'rho_n',80)
call dvdump(v_n,Nx,'v_n',80)
call calc_flux(rho_n, v_n,
& x, FF1_n, FF2_n, Nx, Ng, Gamma, floor)
write(80,*)
write(80,*) '***************************************'
write(80,*) '*** 1. BEFORE CALL TO UPDATE_Q *******'
write(80,*) '***************************************'
call dvdump(FF1_n,Nx,'FF1_n',80)
call dvdump(FF2_n,Nx,'FF2_n',80)
call dvdump(tau_np1,Nx,'tau_np1',80)
call dvdump(tau_n,Nx,'tau_n',80)
call dvdump(S_np1,Nx,'S_np1',80)
call dvdump(S_n,Nx,'S_n',80)
call update_q (tau_np1, S_np1, tau_n, S_n, x,
& FF1_n, FF2_n, delta_tb2, Nx, Ng)
write(80,*)
write(80,*) '***************************************'
write(80,*) '*** 1. BEFORE CALL TO UPDATE_BOUNDARY *'
write(80,*) '***************************************'
call dvdump(tau_np1,Nx,'tau_np1',80)
call dvdump(tau_n,Nx,'tau_n',80)
call dvdump(S_np1,Nx,'S_np1',80)
call dvdump(S_n,Nx,'S_n',80)
call update_boundary(tau_np1, S_np1, Nx, Ng)
write(80,*)
write(80,*) '***************************************'
write(80,*) '*** 1. AFTER CALL TO UPDATE_BOUNDARY *'
write(80,*) '***************************************'
call dvdump(tau_np1,Nx,'tau_np1',80)
call dvdump(S_np1,Nx,'S_np1',80)
write(80,*)
do i = 1, Nx
tau_np1(i) = max(tau_np1(i),floor+abs(S_np1(i)))
ret= calc_prim(tau_np1(i), S_np1(i),
& rho_np1(i), v_np1(i), P_np1(i), Gamma)
if (ret .lt. 0) then
write(0,*) cdnm, ': Stopped at the 1st R-K iteration'
write(0,*) cdnm, ': position x(',i,')=', x(i)
stop
endif
enddo
c-------------------------------------------------------------
c Second step. Computes quantities at full-advanced time step.
c Numerical fluxes are computed using half-step values.
c------------------------------------------------------------
write(80,*) '+++++++++++++++++++++++++++++++++++++++'
write(80,*) '+++ 2. BEFORE CALL TO CALC_FLUX +++++++'
write(80,*) '+++++++++++++++++++++++++++++++++++++++'
call dvdump(tau_np1,Nx,'tau_np1',80)
call dvdump(S_np1,Nx,'S_np1',80)
call dvdump(rho_np1,Nx,'rho_np1',80)
call dvdump(v_np1,Nx,'v_np1',80)
call dvdump(P_np1,Nx,'P_np1',80)
write(80,*)
call calc_flux(rho_np1, v_np1,
& x, FF1_n, FF2_n, Nx, Ng, Gamma, floor)
write(80,*) '+++++++++++++++++++++++++++++++++++++++'
write(80,*) '+++ 2. BEFORE CALL TO UPDATE_Q +++++++'
write(80,*) '+++++++++++++++++++++++++++++++++++++++'
call dvdump(FF1_n,Nx,'FF1_n',80)
call dvdump(FF2_n,Nx,'FF2_n',80)
call dvdump(tau_np1,Nx,'tau_np1',80)
call dvdump(tau_n,Nx,'tau_n',80)
call dvdump(S_np1,Nx,'S_np1',80)
call dvdump(S_n,Nx,'S_n',80)
write(80,*)
call update_q (tau_np1, S_np1, tau_n, S_n, x,
& FF1_n, FF2_n, delta_t, Nx, Ng)
write(80,*) '+++++++++++++++++++++++++++++++++++++++'
write(80,*) '+++ 2. BEFORE CALL TO UPDATE_BOUNDARY +'
write(80,*) '+++++++++++++++++++++++++++++++++++++++'
write(80,*)
call dvdump(tau_np1,Nx,'tau_np1',80)
call dvdump(S_np1,Nx,'S_np1',80)
write(80,*)
call update_boundary(tau_np1, S_np1, Nx, Ng)
write(80,*) '***************************************'
write(80,*) '*** 2. AFTER CALL TO UPDATE_BOUNDARY *'
write(80,*) '***************************************'
write(80,*)
call dvdump(tau_np1,Nx,'tau_np1',80)
call dvdump(S_np1,Nx,'S_np1',80)
write(80,*)
do i = 1, Nx
tau_np1(i) = max(tau_np1(i),floor+abs(S_np1(i)))
ret= calc_prim(tau_np1(i), S_np1(i),
& rho_np1(i), v_np1(i), P_np1(i), Gamma)
if (ret .lt. 0) then
write(0,*) cdnm, ': Stopped at the 2nd R-K iteration'
write(0,*) cdnm, ': position x(',i,')=', x(i)
stop
endif
enddo
write(80,*) '+++++++++++++++++++++++++++++++++++++++'
write(80,*) '+++ 2. AFTER PRIMIITVE COMPUTATION +'
write(80,*) '+++++++++++++++++++++++++++++++++++++++'
write(80,*)
call dvdump(tau_np1,Nx,'tau_np1',80)
call dvdump(S_np1,Nx,'S_np1',80)
call dvdump(rho_np1,Nx,'rho_np1',80)
call dvdump(v_np1,Nx,'v_np1',80)
call dvdump(P_np1,Nx,'P_np1',80)
c if ( independent_residual ) then
c call indep_res( q_np1, q_n, x, res_n, Nx, Ng, lambda)
c endif
return