|
stepgrid_geo.f.html |
 |
|
Source file: stepgrid_geo.f
|
|
Directory: /home/rjl/git/rjleveque/clawpack-4.x/geoclaw/2d/lib
|
|
Converted: Tue Jul 26 2011 at 12:58:45
using clawcode2html
|
|
This documentation file will
not reflect any later changes in the source file.
|
c
c -------------------------------------------------------------
c
subroutine stepgrid(q,fm,fp,gm,gp,mitot,mjtot,mbc,dt,dtnew,dx,dy,
& nvar,xlow,ylow,time,mptr,maux,aux)
c
c
c ::::::::::::::::::: STEPGRID ::::::::::::::::::::::::::::::::::::
c take a time step on a single grid. overwrite solution array q.
c A modified version of the clawpack routine step2 is used.
c
c return fluxes in fm,fp and gm,gp.
c patch has room for ghost cells (mbc of them) around the grid.
c everything is the enlarged size (mitot by mjtot).
c
c mbc = number of ghost cells (= lwidth)
c mptr = grid number (for debugging)
c xlow,ylow = lower left corner of enlarged grid (including ghost cells).
c dt = incoming time step
c dx,dy = mesh widths for this grid
c dtnew = return suggested new time step for this grid's soln.
c
c
c
c This version of stepgrid, stepgrid_geo.f allows output on
c fixed grids specified in setfixedgrids.data
c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
use geoclaw_module
implicit double precision (a-h,o-z)
external rpn2,rpt2
include "call.i"
include "fixedgrids.i"
common /comxyt/ dtcom,dxcom,dycom,tcom,icom,jcom
parameter (msize=max1d+4)
parameter (mwork=msize*(maxvar*maxvar + 13*maxvar + 3*maxaux +2))
dimension q(mitot,mjtot,nvar)
dimension fp(mitot,mjtot,nvar),gp(mitot,mjtot,nvar)
dimension fm(mitot,mjtot,nvar),gm(mitot,mjtot,nvar)
dimension aux(mitot,mjtot,maux)
dimension work(mwork)
logical debug, dump
data debug/.false./, dump/.false./
c
c # set tcom = time. This is in the common block comxyt that could
c # be included in the Riemann solver, for example, if t is explicitly
c # needed there.
tcom = time
if (dump) then
write(*,*)" dumping grid ",mptr
do i = 1, mitot
do j = 1, mjtot
c write(outunit,545) i,j,(q(i,j,ivar),ivar=1,nvar)
write(*,545) i,j,(q(i,j,ivar),ivar=1,nvar)
545 format(2i4,4e15.7)
end do
end do
endif
c
meqn = nvar
mx = mitot - 2*mbc
my = mjtot - 2*mbc
maxm = max(mx,my) !# size for 1d scratch array
mbig = maxm
xlowmbc = xlow + mbc*dx
ylowmbc = ylow + mbc*dy
c # method(2:7) and mthlim
c # are set in the amr2ez file (read by amr)
c
method(1) = 0
c
c
c # fluxes initialized in step2
c
mwork0 = (maxm+2*mbc)*(12*meqn + mwaves + meqn*mwaves + 2)
c
if (mwork .lt. mwork0) then
write(outunit,*) 'CLAW2 ERROR... mwork must be increased to ',
& mwork0
write(* ,*) 'CLAW2 ERROR... mwork must be increased to ',
& mwork0
stop
endif
c
c # partition work array into pieces needed for local storage in
c # step2 routine. Find starting index of each piece:
c
i0faddm = 1
i0faddp = i0faddm + (maxm+2*mbc)*meqn
i0gaddm = i0faddp + (maxm+2*mbc)*meqn
i0gaddp = i0gaddm + 2*(maxm+2*mbc)*meqn
i0q1d = i0gaddp + 2*(maxm+2*mbc)*meqn
i0dtdx1 = i0q1d + (maxm+2*mbc)*meqn
i0dtdy1 = i0dtdx1 + (maxm+2*mbc)
i0aux1 = i0dtdy1 + (maxm+2*mbc)
i0aux2 = i0aux1 + (maxm+2*mbc)*maux
i0aux3 = i0aux2 + (maxm+2*mbc)*maux
c
c
i0next = i0aux3 + (maxm+2*mbc)*maux !# next free space
mused = i0next - 1 !# space already used
mwork1 = mwork - mused !# remaining space (passed to step2)
c
c::::::::::::::::::::::::Fixed Grid Output:::::::::::::::::::::::::::::::::
tc0=time !# start of computational step
tcf=tc0+dt !# end of computational step
c # see if any f-grids should be written out
do ng=1,mfgrids
if (tc0.gt.tstartfg(ng).and.ilastoutfg(ng).lt.noutfg(ng)) then
c # fgrid ng may need to be written out
c # find the first output number that has not been written out and
c # find the first output number on a fixed grid that is >= tc0
c # which will not be written out
if (dtfg(ng).gt.0.d0) then
ioutfgend= 1+max(0,nint((tc0-tstartfg(ng))/dtfg(ng)))
else
ioutfgend=1
endif
ioutfgend=min(ioutfgend,noutfg(ng))
ioutfgstart=ilastoutfg(ng)+1
c # write-out fgrid times that are less than tc0, and have not been written yet
c # these should be the most accurate values at any given point in the fgrid
c # since tc0> output time
do ioutfg=ioutfgstart,ioutfgend
toutfg=tstartfg(ng)+(ioutfg-1)*dtfg(ng)
if (toutfg.lt.tc0) then
c # write out the solution for fixed grid ng
i0=i0fg(ng)
i02=i0fg2(ng)
c # test if arrival times should be output
ioutflag = ioutarrivaltimes(ng)*
& (noutfg(ng)-ilastoutfg(ng))
call fgridout(fgridearly(i0),fgridlate(i0),
& fgridoften(i02),xlowfg(ng),xhifg(ng),ylowfg(ng),
& yhifg(ng),mxfg(ng),myfg(ng),
& mfgridvars(ng),mfgridvars2(ng),toutfg,
& ioutfg,ng,ioutarrivaltimes(ng),ioutflag)
tlastoutfg(ng)=toutfg
ilastoutfg(ng)=ilastoutfg(ng)+1
endif
enddo
endif
enddo
c::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
call b4step2(mx,my,mbc,mx,my,nvar,q,
& xlowmbc,ylowmbc,dx,dy,time,dt,maux,aux)
c::::::::::::::::::::::::FIXED GRID DATA before step:::::::::::::::::::::::
c # fill in values at fixed grid points effected at time tc0
do ng=1,mfgrids
if ((xlowfg(ng).lt.xlowmbc+mx*dx.and.xhifg(ng).gt.xlowmbc).and.
& (ylowfg(ng).lt.ylowmbc+my*dy.and.yhifg(ng).gt.ylowmbc).and.
& (ilastoutfg(ng).lt.noutfg(ng).and.tcf.ge.tstartfg(ng))) then
if (tlastoutfg(ng)+dtfg(ng).ge.tc0.and.
& tlastoutfg(ng)+dtfg(ng).le.tcf) then
c # fixedgrid ng has an output time within [tc0,tcf] interval
c # and it overlaps this computational grid spatially
i0=i0fg(ng) !# index into the ng grid in the work array
call fgridinterp(fgridearly(i0),xlowfg(ng),ylowfg(ng),
& xhifg(ng),yhifg(ng),dxfg(ng),dyfg(ng),mxfg(ng),myfg(ng),
& tc0,mfgridvars(ng),q,nvar,mx,my,mbc,dx,dy,nvar,xlowmbc,
& ylowmbc,maux,aux,ioutarrivaltimes(ng),ioutsurfacemax(ng),0)
c # routine to spatially interpolate computational solution
c # at tc0 to the fixed grid spatial points,
c #saving solution, variables and tc0 at every grid point
endif
c # set maxima or minima if this is a new coarse step
if (tc0.ge.tcfmax) then
if (ioutsurfacemax(ng)+ioutarrivaltimes(ng).gt.0) then
i0=i0fg2(ng)
call fgridinterp(fgridoften(i0),xlowfg(ng),ylowfg(ng),
& xhifg(ng),yhifg(ng),dxfg(ng),dyfg(ng),mxfg(ng),myfg(ng),
& tc0,mfgridvars2(ng),q,nvar,mx,my,mbc,dx,dy,nvar,xlowmbc,
& ylowmbc,maux,aux,ioutarrivaltimes(ng),ioutsurfacemax(ng),2)
endif
endif
endif
enddo
tcfmax=max(tcfmax,tcf)
c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
c
c # take one step on the conservation law:
c
call step2(mbig,mx,my,nvar,maux,
& mbc,mx,my,
& q,aux,dx,dy,dt,cflgrid,
& fm,fp,gm,gp,
& work(i0faddm),work(i0faddp),
& work(i0gaddm),work(i0gaddp),
& work(i0q1d),work(i0dtdx1),work(i0dtdy1),
& work(i0aux1),work(i0aux2),work(i0aux3),
& work(i0next),mwork1,rpn2,rpt2)
c
c
mptr_level = node(nestlevel,mptr)
write(outunit,811) mptr, mptr_level, cflgrid
811 format(" Courant # of grid ",i5," level",i3," is ",d12.4)
c
cflmax = dmax1(cflmax,cflgrid)
cfl_level = dmax1(cfl_level,cflgrid)
c
c # update q
dtdx = dt/dx
dtdy = dt/dy
do 50 m=1,nvar
do 50 i=mbc+1,mitot-mbc
do 50 j=mbc+1,mjtot-mbc
if (mcapa.eq.0) then
c
c # no capa array. Standard flux differencing:
q(i,j,m) = q(i,j,m)
& - dtdx * (fm(i+1,j,m) - fp(i,j,m))
& - dtdy * (gm(i,j+1,m) - gp(i,j,m))
else
c # with capa array.
q(i,j,m) = q(i,j,m)
& - (dtdx * (fm(i+1,j,m) - fp(i,j,m))
& + dtdy * (gm(i,j+1,m) - gp(i,j,m))) / aux(i,j,mcapa)
endif
50 continue
c
c # Copied here from b4step2 since need to do before saving to qc1d:
do i=1,mitot
do j=1,mjtot
if (q(i,j,1).lt.drytolerance) then
q(i,j,1) = max(q(i,j,1),0.d0)
do m=2,nvar
q(i,j,m)=0.d0
enddo
endif
enddo
enddo
c
if (method(5).eq.1) then
c # with source term: use Godunov splitting
call src2(mx,my,nvar,mbc,mx,my,xlowmbc,ylowmbc,dx,dy,
& q,maux,aux,time,dt)
endif
c ::::::::::::::::::::::::Fixed Grid data afterstep:::::::::::::::::::::::
c # fill in values at fixed grid points effected at time tcf
do ng=1,mfgrids
if ((xlowfg(ng).lt.xlowmbc+mx*dx.and.xhifg(ng).gt.xlowmbc).and.
& (ylowfg(ng).lt.ylowmbc+my*dy.and.yhifg(ng).gt.ylowmbc).and.
& (ilastoutfg(ng).lt.noutfg(ng).and.tcf.ge.tstartfg(ng))) then
if (tlastoutfg(ng)+dtfg(ng).ge.tc0
& .and.tlastoutfg(ng)+dtfg(ng).le.tcf) then
c # fixedgrid ng has an output time within [tc0,tcf] interval
c # and it overlaps this computational grid spatially
i0=i0fg(ng) !# index into the ng grid in the work array
call fgridinterp(fgridlate(i0),xlowfg(ng),ylowfg(ng),
& xhifg(ng),yhifg(ng),dxfg(ng),dyfg(ng),mxfg(ng),myfg(ng),
& tcf,mfgridvars(ng),q,nvar,mx,my,mbc,dx,dy,nvar,xlowmbc,
& ylowmbc,maux,aux,ioutarrivaltimes(ng),ioutsurfacemax(ng),0)
c # routine to interpolate solution
c # at tcf to the fixed grid storage array,
c #saving solution and tcf at every grid point
endif
c # fill in values for eta if they need to be saved for later checking max/mins
c # check for arrival times
if (ioutsurfacemax(ng)+ioutarrivaltimes(ng).gt.0) then
i0=i0fg2(ng)
call fgridinterp(fgridoften(i0),xlowfg(ng),ylowfg(ng),
& xhifg(ng),yhifg(ng),dxfg(ng),dyfg(ng),mxfg(ng),myfg(ng),
& tc0,mfgridvars2(ng),q,nvar,mx,my,mbc,dx,dy,nvar,xlowmbc,
& ylowmbc,maux,aux,ioutarrivaltimes(ng),ioutsurfacemax(ng),1)
endif
endif
enddo
c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
c # output fluxes for debugging purposes:
if (debug) then
write(dbugunit,*)" fluxes for grid ",mptr
do 830 i = mbc+1, mitot-1
do 830 j = mbc+1, mjtot-1
write(dbugunit,831) i,j,fm(i,j,1),fp(i,j,1),
. gm(i,j,1),gp(i,j,1)
do 830 m = 2, meqn
write(dbugunit,832) fm(i,j,m),fp(i,j,m),
. gm(i,j,m),gp(i,j,m)
831 format(2i4,4d16.6)
832 format(8x,4d16.6)
830 continue
endif
c
c
c For variable time stepping, use max speed seen on this grid to
c choose the allowable new time step dtnew. This will later be
c compared to values seen on other grids.
c
if (cflgrid .gt. 0.d0) then
dtnew = dt*cfl/cflgrid
else
c # velocities are all zero on this grid so there's no
c # time step restriction coming from this grid.
dtnew = rinfinity
endif
c # give a warning if Courant number too large...
c
if (cflgrid .gt. cflv1) then
write(*,810) cflgrid,cflv1,mptr,mptr_level
write(outunit,810) cflgrid, cflv1,mptr,mptr_level
810 format('*** WARNING *** Courant number =', d12.4,
& ' is larger than input cfl_max = ', d12.4,
& ' on grid ',i3, ' level ',i3)
endif
c
!-- if (dump) then
!-- write(*,*)" at end of stepgrid: dumping grid ",mptr
!-- do i = 1, mitot
!-- do j = 1, mjtot
!--c write(outunit,545) i,j,(q(i,j,ivar),ivar=1,nvar)
!-- write(*,545) i,j,(q(i,j,ivar),ivar=1,nvar)
!-- end do
!-- end do
!-- endif
c
return
end