Update remaining source

Author: mandli

src/1d_classic/shallow/claw1.f

@@ -6,7 +6,7 @@ subroutine claw1(meqn,mwaves,maux,mbc,mx,
 c     ==============================================================
 c
 c  Solves a hyperbolic system of conservation laws in one space dimension
-c  of the general form 
+c  of the general form
 c
 c     capa * q_t + A q_x = psi
 c
@@ -24,7 +24,7 @@ subroutine claw1(meqn,mwaves,maux,mbc,mx,
 c
 c  The user must supply the following subroutines:
 c
-c    bc1, rp1        subroutines specifying the boundary conditions and 
+c    bc1, rp1        subroutines specifying the boundary conditions and
 c                    Riemann solver.
 c                    These are described in greater detail below.
 c
@@ -44,7 +44,7 @@ subroutine claw1(meqn,mwaves,maux,mbc,mx,
 c  These routines must be declared EXTERNAL in the main program.
 c  For description of the calling sequences, see below.
 c
-c  Dummy routines b4step1.f and src1.f are available in 
+c  Dummy routines b4step1.f and src1.f are available in
 c       claw/clawpack/1d/lib
 c
 c
@@ -74,12 +74,12 @@ subroutine claw1(meqn,mwaves,maux,mbc,mx,
 c       physical domain.  In addition there are mbc grid cells
 c       along each edge of the grid that are used for boundary
 c       conditions.
-c 
-c    q(meqn, 1-mbc:mx+mbc) 
+c
+c    q(meqn, 1-mbc:mx+mbc)
 c        On input:  initial data at time tstart.
 c        On output: final solution at time tend.
 c        q(m,i) = value of mth component in the i'th cell.
-c        Values within the physical domain are in q(m,i) 
+c        Values within the physical domain are in q(m,i)
 c                for i = 1,2,...,mx
 c        mbc extra cells on each end are needed for boundary conditions
 c        as specified in the routine bc1.
@@ -102,15 +102,15 @@ subroutine claw1(meqn,mwaves,maux,mbc,mx,
 c            stored in aux(mcapa,i).
 c            In this case we require method(7).ge.mcapa.
 c
-c    dx = grid spacing in x.  
+c    dx = grid spacing in x.
 c         (for a computation in ax <= x <= bx,  set dx = (bx-ax)/mx.)
 c
 c    tstart = initial time.
 c
 c    tend = Desired final time (on input).
 c              If tend<tstart, then claw1 returns after a single successful
 c                 time step has been taken (single-step mode).
-c              Otherwise, as many steps are taken as needed to reach tend, 
+c              Otherwise, as many steps are taken as needed to reach tend,
 c                 up to a maximum of nv(1).
 c         = Actual time reached (on output).
 c
@@ -152,13 +152,13 @@ subroutine claw1(meqn,mwaves,maux,mbc,mx,
 c                       comes from the previous step, the Courant number will
 c                       not in general be exactly equal to the desired value
 c                       If the actual Courant number in the next step is
-c                       greater than 1, then this step is redone with a 
+c                       greater than 1, then this step is redone with a
 c                       smaller dt.
 c
 c         method(2) = 1 if Godunov's method is to be used, with no 2nd order
 c                       corrections.
 c                   = 2 if second order correction terms are to be added, with
-c                       a flux limiter as specified by mthlim.  
+c                       a flux limiter as specified by mthlim.
 c
 c         method(3)  is not used in one-dimension.
 c
@@ -223,7 +223,7 @@ subroutine claw1(meqn,mwaves,maux,mbc,mx,
 c            If mwork is too small then the program returns with info = 4
 c            and prints the necessary value of mwork to unit 6.
 c
-c            
+c
 c    info = output value yielding error information:
 c         = 0 if normal return.
 c         = 1 if mbc.lt.2
@@ -241,7 +241,7 @@ subroutine claw1(meqn,mwaves,maux,mbc,mx,
 c    User-supplied subroutines
 c    -------------------------
 c
-c    bc1 = subroutine that specifies the boundary conditions.  
+c    bc1 = subroutine that specifies the boundary conditions.
 c         This subroutine should extend the values of q from cells
 c         1:mx to the mbc ghost cells along each edge of the domain.
 c
@@ -290,7 +290,7 @@ subroutine claw1(meqn,mwaves,maux,mbc,mx,
 c         are passed in using auxl and auxr.  Again, in the standard routines
 c         auxl=auxr=aux in the call to rp1.
 c
-c          On output, 
+c          On output,
 c              wave(m,mw,i) is the m'th component of the jump across
 c                              wave number mw in the ith Riemann problem.
 c              s(mw,i) is the wave speed of wave number mw in the
@@ -307,7 +307,7 @@ subroutine claw1(meqn,mwaves,maux,mbc,mx,
 c           into the specification of amdq and apdq.
 c
 c    src1 = subroutine for the source terms that solves the equation
-c               capa * q_t = psi 
+c               capa * q_t = psi
 c           over time dt.
 c
 c           If method(5)=0 then the equation does not contain a source
@@ -338,7 +338,7 @@ subroutine claw1(meqn,mwaves,maux,mbc,mx,
 c                other tasks which must be done every time step.
 c
 c          The form of this subroutine is
-c      
+c
 c  -------------------------------------------------
 c      subroutine b4step1(mbc,mx,meqn,q,xlower,dx,time,dt,maux,aux)
 c      implicit double precision (a-h,o-z)
@@ -353,12 +353,12 @@ subroutine claw1(meqn,mwaves,maux,mbc,mx,
 c
 c  Copyright 1994 -- 2002 R. J. LeVeque
 c
-c  This software is made available for research and instructional use only. 
+c  This software is made available for research and instructional use only.
 c  You may copy and use this software without charge for these non-commercial
 c  purposes, provided that the copyright notice and associated text is
 c  reproduced on all copies.  For all other uses (including distribution of
-c  modified versions), please contact the author at the address given below. 
-c  
+c  modified versions), please contact the author at the address given below.
+c
 c  *** This software is made available "as is" without any assurance that it
 c  *** will work for your purposes.  The software may in fact have defects, so
 c  *** use the software at your own risk.
@@ -370,18 +370,18 @@ subroutine claw1(meqn,mwaves,maux,mbc,mx,
 c    Author:  Randall J. LeVeque
 c             Applied Mathematics
 c             Box 352420
-c             University of Washington, 
+c             University of Washington,
 c             Seattle, WA 98195-2420
 c             [email protected]
 c =========================================================================
 c
 c
 c
-c            
+c
 c    ======================================================================
 c    Beginning of claw1 code
 c    ======================================================================
-c 
+c
 
       use gauges_module, only: num_gauges, update_gauges,
      &                         print_gauges_and_reset_nextLoc
@@ -435,7 +435,7 @@ subroutine claw1(meqn,mwaves,maux,mbc,mx,
       i0wave = i0f + (mx + 2*mbc) * meqn
       i0s = i0wave + (mx + 2*mbc) * meqn * mwaves
       i0dtdx = i0s + (mx + 2*mbc) * mwaves
-      i0qwork = i0dtdx + (mx + 2*mbc) 
+      i0qwork = i0dtdx + (mx + 2*mbc)
       i0amdq = i0qwork + (mx + 2*mbc) * meqn
       i0apdq = i0amdq + (mx + 2*mbc) * meqn
       i0dtdx = i0apdq + (mx + 2*mbc) * meqn
@@ -463,7 +463,7 @@ subroutine claw1(meqn,mwaves,maux,mbc,mx,
 c           # save old q in case we need to retake step with smaller dt:
             call copyq1(meqn,mbc,mx,q,work(i0qwork))
             endif
-c           
+c
    40    continue
          dt2 = dt / 2.d0
          thalf = t + dt2  !# midpoint in time for Strang splitting
@@ -532,8 +532,8 @@ subroutine claw1(meqn,mwaves,maux,mbc,mx,
                 endif
                 dtmin = dmin1(dt,dtmin)
                 dtmax = dmax1(dt,dtmax)
-                
-                    
+
+
               else
                 dt = dtv(2)
               endif
@@ -563,13 +563,13 @@ subroutine claw1(meqn,mwaves,maux,mbc,mx,
                     go to 900
                   endif
                endif
-               
+
          if (.not. topo_finalized) then
              ! modify topo using dtopo arrays:
              call topo_update(t)
              call setaux(mbc,mx,xlower,dx,maux,aux)
          endif
-      
+
 c
 c        # see if we are done:
          nv(2) = nv(2) + 1
@@ -578,7 +578,7 @@ subroutine claw1(meqn,mwaves,maux,mbc,mx,
   100    continue
 c
   900  continue
-c 
+c
 c      # return information
 c
        if (method(1).eq.1 .and. t.lt.tend .and. nv(2) .eq. maxn) then
@@ -602,5 +602,5 @@ subroutine claw1(meqn,mwaves,maux,mbc,mx,
           call print_gauges_and_reset_nextLoc(ii,meqn)
        end do
 
-       return 
+       return
        end

src/2d/shallow/b4step2.f90

@@ -1,5 +1,6 @@
 ! ============================================
-subroutine b4step2(mbc,mx,my,meqn,q,xlower,ylower,dx,dy,t,dt,maux,aux,actualstep)
+subroutine b4step2(mbc, mx, my, meqn, q, xlower, ylower, dx, dy, t, dt,        &
+                   maux, aux, actualstep)
 ! ============================================
 !
 ! # called before each call to step
@@ -36,7 +37,7 @@ subroutine b4step2(mbc,mx,my,meqn,q,xlower,ylower,dx,dy,t,dt,maux,aux,actualstep
     real(kind=8), intent(inout) :: xlower, ylower, dx, dy, t, dt
     real(kind=8), intent(inout) :: q(meqn,1-mbc:mx+mbc,1-mbc:my+mbc)
     real(kind=8), intent(inout) :: aux(maux,1-mbc:mx+mbc,1-mbc:my+mbc)
-    logical, intent (in) :: actualstep
+    logical, intent(in) :: actualstep
 
     ! Local storage
     integer :: index,i,j,k,dummy
@@ -95,6 +96,8 @@ subroutine b4step2(mbc,mx,my,meqn,q,xlower,ylower,dx,dy,t,dt,maux,aux,actualstep
     endif
 
     ! Set wind and pressure aux variables for this grid
-    call set_storm_fields(maux,mbc,mx,my,xlower,ylower,dx,dy,t,aux)
+    if (actualstep) then
+        call set_storm_fields(maux,mbc,mx,my,xlower,ylower,dx,dy,t,aux)
+    end if
 
 end subroutine b4step2

src/2d/shallow/multilayer/b4step2.f90

@@ -1,14 +1,15 @@
-subroutine b4step2(mbc,mx,my,meqn,q,xlower,ylower,dx,dy,t,dt,maux,aux)
+subroutine b4step2(mbc, mx, my, meqn, q, xlower, ylower, dx, dy, t, dt,        &
+                   maux, aux, actualstep)
 ! ============================================
-! 
+!
 ! # called before each call to step
 ! # use to set time-dependent aux arrays or perform other tasks.
-! 
+!
 ! This particular routine sets negative values of q(1,i,j) to zero,
 ! as well as the corresponding q(m,i,j) for m=1,meqn.
 ! This is for problems where q(1,i,j) is a depth.
 ! This should occur only because of rounding error.
-! 
+!
 ! Also calls movetopo if topography might be moving.
 
     use amr_module, only: xlowdomain => xlower
@@ -27,14 +28,15 @@ subroutine b4step2(mbc,mx,my,meqn,q,xlower,ylower,dx,dy,t,dt,maux,aux)
 
     use multilayer_module, only: num_layers, KAPPA_UNIT, dry_tolerance
     use multilayer_module, only: check_richardson, richardson_tolerance
-    
+
     implicit none
-    
+
     ! Subroutine arguments
     integer, intent(in) :: mbc,mx,my,meqn,maux
     real(kind=8), intent(in) :: xlower, ylower, dx, dy, t, dt
     real(kind=8), intent(inout) :: q(meqn,1-mbc:mx+mbc,1-mbc:my+mbc)
     real(kind=8), intent(inout) :: aux(maux,1-mbc:mx+mbc,1-mbc:my+mbc)
+    logical, intent(in) :: actualstep
 
     ! Local storage
     integer :: index,i,j,k
@@ -58,14 +60,17 @@ subroutine b4step2(mbc,mx,my,meqn,q,xlower,ylower,dx,dy,t,dt,maux,aux)
     end forall
 
     ! Move the topography if needed
-    if (aux_finalized < 2) then
+    if (aux_finalized < 2 .and. actualstep) then
         call setaux(mbc, mx, my, xlower, ylower, dx, dy, maux, aux)
-    endif 
+    endif
 
-    call set_storm_fields(maux,mbc,mx,my,xlower,ylower,dx,dy,t,aux)
+    ! Set wind and pressure aux variables for this grid
+    if (actualstep) then
+        call set_storm_fields(maux,mbc,mx,my,xlower,ylower,dx,dy,t,aux)
+    end if
 
     ! Check Richardson number -- Only implemented for 2 layers
-    if (num_layers == 2 .and. check_richardson) then
+    if (num_layers == 2 .and. check_richardson .and. actualstep) then
         do i=1,mx
             do j=1,my
                 dry_state = .false.

src/2d/shallow/multilayer/stepgrid.f

@@ -28,7 +28,7 @@ subroutine stepgrid(q,fm,fp,gm,gp,mitot,mjtot,mbc,dt,dtnew,dx,dy,
 
       use geoclaw_module, only: grav, rho
       use amr_module
-      use fgout_module, only: FGOUT_num_grids, FGOUT_fgrids, 
+      use fgout_module, only: FGOUT_num_grids, FGOUT_fgrids,
      &                        FGOUT_tcfmax, fgout_interp, fgout_grid,
      &                        FGOUT_ttol
       use multilayer_module, only: num_layers, dry_tolerance
@@ -54,7 +54,7 @@ subroutine stepgrid(q,fm,fp,gm,gp,mitot,mjtot,mbc,dt,dtnew,dx,dy,
 c
       FGOUT_tcfmax = -rinfinity
       level = node(nestlevel,mptr)
-      
+
 
       if (dump) then
          write(outunit,*)" at start of stepgrid: dumping grid ",mptr
@@ -115,7 +115,7 @@ subroutine stepgrid(q,fm,fp,gm,gp,mitot,mjtot,mbc,dt,dtnew,dx,dy,
 c
       tc0=time !# start of computational step
       tcf=tc0+dt !# end of computational step
-      
+
 c     Check if fgout interpolation needed before and after step:
 
       allocate(fgout_interp_needed(FGOUT_num_grids))
@@ -128,7 +128,7 @@ subroutine stepgrid(q,fm,fp,gm,gp,mitot,mjtot,mbc,dt,dtnew,dx,dy,
               fgout_interp_needed(ng) = .false.
           else
               fgout_tnext = fgout%output_times(fgout%next_output_index)
-              fgout_interp_needed(ng) =  
+              fgout_interp_needed(ng) =
      &          ((fgout%x_low < xlowmbc + mx * dx) .and.
      &           (fgout%x_hi  > xlowmbc) .and.
      &           (fgout%y_low < ylowmbc + my * dy) .and.
@@ -146,13 +146,13 @@ subroutine stepgrid(q,fm,fp,gm,gp,mitot,mjtot,mbc,dt,dtnew,dx,dy,
 
 
 c::::::::::::::::::::::::fgout Output:::::::::::::::::::::::::::::::::
-c     This has been moved to tick.f, after advancing all patches on 
+c     This has been moved to tick.f, after advancing all patches on
 c     finest level.  No need to check on each patch separately.
 c::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
 
-       call b4step2(mbc,mx,my,nvar,q,
-     &             xlowmbc,ylowmbc,dx,dy,time,dt,maux,aux,actualstep)
-      
+!        call b4step2(mbc,mx,my,nvar,q,
+!      &             xlowmbc,ylowmbc,dx,dy,time,dt,maux,aux,actualstep)
+
 c::::::::::::::::::::::::FGOUT DATA before step:::::::::::::::::::::::
 c     # fill in values at fgout points affected at time tc0
       do ng=1,FGOUT_num_grids
@@ -183,7 +183,7 @@ subroutine stepgrid(q,fm,fp,gm,gp,mitot,mjtot,mbc,dt,dtnew,dx,dy,
      &              q,aux,dx,dy,dt,cflgrid,
      &              fm,fp,gm,gp,rpn2,rpt2)
 c
-c            
+c
         mptr_level = node(nestlevel,mptr)
 
 c       write(outunit,811) mptr, mptr_level, cflgrid
@@ -319,5 +319,3 @@ subroutine stepgrid(q,fm,fp,gm,gp,mitot,mjtot,mbc,dt,dtnew,dx,dy,
 c
       return
       end
-
-