|
filpatch_geo.f.html |
 |
|
Source file: filpatch_geo.f
|
|
Directory: /Users/rjl/git/rjleveque/clawpack-4.6.3/geoclaw/2d/lib
|
|
Converted: Mon Jan 21 2013 at 20:15:56
using clawcode2html
|
|
This documentation file will
not reflect any later changes in the source file.
|
c
c ---------------------------------------------------------------
c
recursive subroutine filrecur(level,nvar,valbig,aux,naux,
1 time,mitot,mjtot,
2 nrowst,ncolst,ilo,ihi,jlo,jhi)
c :::::::::::::::::::::::::::: FILPATCH :::::::::::::::::::::::::;
c
c fill the portion of valbig from rows nrowst
c and cols ncolst
c the patch can also be described by the corners (xlp,ybp) by (xrp,ytp).
c vals are needed at time time, and level level,
c
c first fill with values obtainable from the level level
c grids. if any left unfilled, then enlarge remaining rectangle of
c unfilled values by 1 (for later linear interp), and recusively
c obtain the remaining values from coarser levels.
c
c :::::::::::::::::::::::::::::::::::::::;:::::::::::::::::::::::;
use geoclaw_module
implicit double precision (a-h,o-z)
include "call.i"
logical set, sticksout
dimension valbig(mitot,mjtot,nvar), aux(mitot,mjtot,naux)
c use stack-based scratch arrays instead of alloc, since dont really
c need to save beyond these routines, and to allow dynamic memory resizing
c
c use 1d scratch arrays that are potentially the same size as
c current grid, since may not coarsen.
c need to make it 1d instead of 2 and do own indexing, since
c when pass it in to subroutines they treat it as having different
c dimensions than the max size need to allocate here
c
dimension valcrse((ihi-ilo+2)*(jhi-jlo+2)*nvar) ! NB this is a 1D array
dimension auxcrse((ihi-ilo+2)*(jhi-jlo+2)*naux) ! the +2 is to expand on coarse grid to enclose fine
c
dimension flaguse(ihi-ilo+1,jhi-jlo+1)
logical reloop
logical fineflag((ihi-ilo+2)*(jhi-jlo+2)*nvar)
double precision finemass((ihi-ilo+2)*(jhi-jlo+2))
double precision etacrse((ihi-ilo+2)*(jhi-jlo+2))
double precision velmax((ihi-ilo+2)*(jhi-jlo+2))
double precision velmin((ihi-ilo+2)*(jhi-jlo+2))
double precision slopex((ihi-ilo+2)*(jhi-jlo+2))
double precision slopey((ihi-ilo+2)*(jhi-jlo+2))
integer icount((ihi-ilo+2)*(jhi-jlo+2))
ivalc(i,j,ivar) = i + nrowc*(j - 1)
& + nrowc*ncolc*(ivar-1)
icrse(i,j) = i + nrowc*(j-1)
c
c # index into first component of aux = topo:
iauxc(i,j) = i + nrowc*(j-1)
sticksout(iplo,iphi,jplo,jphi) =
& (iplo .lt. 0 .or. jplo .lt. 0 .or.
& iphi .ge. iregsz(levc) .or. jphi .ge. jregsz(levc))
!-- write(*,*)" entering filrecur with level ",level
!-- write(*,*)" and patch indices ilo,ihi,jlo,jhi ",
!-- & ilo,ihi,jlo,jhi
c write(*,*)" in filrecur for level ",level,mitot,mjtot
c
c We begin by filling values for grids at level level. If all values can be
c filled in this way, we return;
nrowp = ihi - ilo + 1
ncolp = jhi - jlo + 1
c locuse = igetsp(nrowp*ncolp)
hxf = hxposs(level)
hyf = hyposs(level)
xlp = xlower + ilo*hxf
xrp = xlower + (ihi+1)*hxf
ybp = ylower + jlo*hyf
ytp = ylower + (jhi+1)*hyf
call intfil
& (valbig,mitot,mjtot,time,flaguse,nrowst,ncolst,
& ilo,ihi,jlo,jhi,level,nvar,naux)
c & (valbig,mitot,mjtot,time,locuse,nrowst,ncolst,
c
c Trimbd returns set = true if all of the entries are filled (=1.).
c set = false, otherwise. If set = true, then no other levels are
c are required to interpolate, and we return.
c
c Note that the used array is filled entirely in intfil, i.e. the
c marking done there also takes into account the points filled by
c the boundary conditions. bc2amr will be called later, after all 4
c boundary pieces filled.
c call trimbd(alloc(locuse),nrowp,ncolp,set,il,ir,jb,jt)
call trimbd(flaguse,nrowp,ncolp,set,il,ir,jb,jt)
if (set) go to 90 ! all done except for bcs
c
c otherwise make recursive calls to coarser levels to fill remaining unset points
c
if (level .eq. 1) then
write(outunit,*)" error in filrecur - level 1 not set"
write(outunit,900) nrowst,ncolst
write(*,*)" error in filrecur - level 1 not set"
write(*,*)" should not need more recursion "
write(*,*)" to set patch boundaries"
write(*,900) nrowst,ncolst
900 format("start at row: ",i4," col ",i4)
stop
endif
c set = false. we will have to interpolate some values from coarser
c levels. We begin by initializing the level level arrays, so that we can use
c purely recursive formulation for interpolating.
levc = level - 1
hxc = hxposs(levc)
hyc = hyposs(levc)
isl = il + ilo - 1
isr = ir + ilo - 1
jsb = jb + jlo - 1
jst = jt + jlo - 1
c
c coarsen
lratiox = intratx(levc)
lratioy = intraty(levc)
iplo = (isl-lratiox +nghost*lratiox)/lratiox - nghost
jplo = (jsb-lratioy +nghost*lratioy)/lratioy - nghost
iphi = (isr+lratiox )/lratiox
jphi = (jst+lratioy )/lratioy
xlc = xlower + iplo*hxc
ybc = ylower + jplo*hyc
xrc = xlower + (iphi+1)*hxc
ytc = ylower + (jphi+1)*hyc
nrowc = iphi - iplo + 1
ncolc = jphi - jplo + 1
ntot = nrowc*ncolc*(nvar+naux)
c write(*,876) nrowc,ncolc, ihi-ilo+2,jhi-jlo+2
876 format(" needed coarse grid size ",2i5," allocated ",2i5)
if (nrowc .gt. ihi-ilo+2 .or. ncolc .gt. jhi-jlo+2) then
write(*,*)" did not make big enough work space in filrecur "
write(*,*)" need coarse space with nrowc,ncolc ",nrowc,ncolc
write(*,*)" coarser level is ",levc," nghost ",nghost
write(*,*)" with ratios lratiox lratioy ",lratiox,lratioy
write(*,*)" made space for ilo,ihi,jlo,jhi ",ilo,ihi,jlo,jhi
write(*,*)" isl,isr,jsb,jst ",isl,isr,jsb,jst
write(*,*)"iplo,iphi,jplo,jphi ",iplo,iphi,jplo,jphi
write(*,*)" orig il,ir,jb,jt ",il,ir,jb,jt
write(*,*)"filpatch called with mitot,mjtot ",mitot,mjtot
call outtre(1,.false.,nvar,naux)
stop
endif
c loccrse = igetsp(ntot)
c locauxc = loccrse + nrowc*ncolc*nvar
if (naux.gt.0) then
maxmx = nrowc - 2*nghost
mx = maxmx
maxmy = ncolc - 2*nghost
my = maxmy
xl = xlc + nghost*hxc
yb = ybc + nghost*hyc
call setaux(maxmx,maxmy,nghost,mx,my,xl,yb,hxc,hyc,
& naux,auxcrse)
c & naux,alloc(locauxc))
endif
if ((xperdom .or. (yperdom .or. spheredom)) .and.
& sticksout(iplo,iphi,jplo,jphi)) then
call prefilrecur(levc,nvar,valcrse,auxcrse,
1 naux,time,nrowc,ncolc,1,1,
2 iplo,iphi,jplo,jphi)
else
c call filpatch2(levc,nvar,alloc(loccrse),alloc(locauxc),naux,
call filrecur(levc,nvar,valcrse,auxcrse,naux,
1 time,nrowc,ncolc,1,1,
2 iplo,iphi,jplo,jphi)
endif
c interpolate back up
20 continue
* !loop through coarse cells determining intepolation slopes
* !these will be saved for fine grid loop
* !prevents finding the same slope possibly lratiox*lratioy times
* !all fine gid depths will be found before any momentum
reloop = .false.
toldry= drytolerance
do ic = 2, nrowc-1
do jc = 2, ncolc-1
icount(icrse(ic,jc)) = 0
finemass(icrse(ic,jc)) = 0.d0
do ivar=1,nvar
fineflag(ivalc(ic,jc,ivar)) = .false.
enddo
* !find interpolation slope for eta = q(:,1)+ aux(:,1)
do i=-1,1
etacrse(icrse(ic+i,jc)) = valcrse(ivalc(ic+i,jc,1))
& + auxcrse(iauxc(ic+i,jc))
if (valcrse(ivalc(ic+i,jc,1)).lt.toldry) then
etacrse(icrse(ic+i,jc)) = sealevel
endif
enddo
s1 = etacrse(icrse(ic,jc))- etacrse(icrse(ic-1,jc))
s2 = etacrse(icrse(ic+1,jc))- etacrse(icrse(ic,jc))
if (s1*s2.le.0) then
slopex(icrse(ic,jc))= 0.d0
else
slopex(icrse(ic,jc))=dmin1(dabs(s1),dabs(s2))*dsign(1.d0,
& etacrse(icrse(ic+1,jc))- etacrse(icrse(ic-1,jc)))
endif
do j=-1,1
etacrse(icrse(ic,jc+j)) = valcrse(ivalc(ic,jc+j,1))
& + auxcrse(iauxc(ic,jc+j))
if (valcrse(ivalc(ic,jc+j,1)).lt.toldry) then
etacrse(icrse(ic,jc+j)) = sealevel
endif
enddo
s1 = etacrse(icrse(ic,jc))- etacrse(icrse(ic,jc-1))
s2 = etacrse(icrse(ic,jc+1))- etacrse(icrse(ic,jc))
if (s1*s2.le.0) then
slopey(icrse(ic,jc))= 0.d0
else
slopey(icrse(ic,jc))=dmin1(dabs(s1),dabs(s2))*dsign(1.d0,
& etacrse(icrse(ic,jc+1))- etacrse(icrse(ic,jc-1)))
endif
end do
end do
* !loop through patch: note this includes multiple coarse cells
do iff = 1,nrowp
ic = 2 + (iff - (isl - ilo) - 1)/lratiox
eta1 = (-0.5d0+dble(mod(iff-1,lratiox)))/dble(lratiox)
do jf = 1,ncolp
jc = 2 + (jf - (jsb - jlo) - 1)/lratioy
eta2 = (-0.5d0+dble(mod(jf -1,lratioy)))/dble(lratioy)
cc = icrse(ic,jc)
c flag = alloc(iadflag(iff,jf))
flag = flaguse(iff,jf)
if (flag .eq. 0.0) then
* !interp. from coarse cells to fine grid to find eta
icount(icrse(ic,jc)) = icount(icrse(ic,jc)) + 1
etafine = etacrse(icrse(ic,jc))
& + eta1*slopex(icrse(ic,jc))+eta2*slopey(icrse(ic,jc))
hfine = max(etafine - aux(iff+nrowst-1,jf+ncolst-1,1),
& 0.d0)
valbig(iff+nrowst-1,jf+ncolst-1,1) = hfine
finemass(icrse(ic,jc)) = finemass(icrse(ic,jc)) +
& valbig(iff+nrowst-1,jf+ncolst-1,1)
if (valbig(iff+nrowst-1,jf+ncolst-1,1).lt.toldry) then
fineflag(ivalc(ic,jc,1)) = .true.
reloop = .true.
endif
endif
enddo
enddo
c ! determine momentum
do ivar = 2,nvar
* !find interpolation slope for momentum = q(:,ivar)
do ic = 2, nrowc-1
do jc = 2, ncolc-1
s1 = valcrse(ivalc(ic,jc,ivar))
& - valcrse(ivalc(ic-1,jc,ivar))
s2 = valcrse(ivalc(ic+1,jc,ivar))
& - valcrse(ivalc(ic,jc,ivar))
if (s1*s2.le.0) then
slopex(icrse(ic,jc))= 0.d0
else
slopex(icrse(ic,jc))=dmin1(dabs(s1),dabs(s2))
& *dsign(1.d0, valcrse(ivalc(ic+1,jc,ivar))
& - valcrse(ivalc(ic-1,jc,ivar)))
endif
s1 = valcrse(ivalc(ic,jc,ivar))
& - valcrse(ivalc(ic,jc-1,ivar))
s2 = valcrse(ivalc(ic,jc+1,ivar))
& - valcrse(ivalc(ic,jc,ivar))
if (s1*s2.le.0) then
slopey(icrse(ic,jc))= 0.d0
else
slopey(icrse(ic,jc))=dmin1(dabs(s1),dabs(s2))
& *dsign(1.d0, valcrse(ivalc(ic,jc+1,ivar))
& - valcrse(ivalc(ic,jc-1,ivar)))
endif
if (valcrse(ivalc(ic,jc,1)).gt.toldry) then
velmax(icrse(ic,jc)) = valcrse(ivalc(ic,jc,ivar))
& /valcrse(ivalc(ic,jc,1))
velmin(icrse(ic,jc)) = valcrse(ivalc(ic,jc,ivar))
& /valcrse(ivalc(ic,jc,1))
else
velmax(icrse(ic,jc)) = 0.d0
velmin(icrse(ic,jc)) = 0.d0
endif
* !look for bounds on velocity to avoid generating new extrema
* !necessary since interpolating momentum linearly
* !yet depth is not interpolated linearly
do ii = -1,1,2
if (valcrse(ivalc(ic+ii,jc,1)).gt.toldry) then
velmax(icrse(ic,jc)) = max(velmax(icrse(ic,jc))
& ,valcrse(ivalc(ic+ii,jc,ivar))
& /valcrse(ivalc(ic+ii,jc,1)))
velmin(icrse(ic,jc)) = min(velmin(icrse(ic,jc))
& ,valcrse(ivalc(ic+ii,jc,ivar))
& /valcrse(ivalc(ic+ii,jc,1)))
endif
if (valcrse(ivalc(ic,jc+ii,1)).gt.toldry) then
velmax(icrse(ic,jc)) = max(velmax(icrse(ic,jc))
& ,valcrse(ivalc(ic,jc+ii,ivar))
& /valcrse(ivalc(ic,jc+ii,1)))
velmin(icrse(ic,jc)) = min(velmin(icrse(ic,jc))
& ,valcrse(ivalc(ic,jc+ii,ivar))
& /valcrse(ivalc(ic,jc+ii,1)))
endif
enddo
end do
end do
* !determine momentum in fine cells
do iff = 1,nrowp
ic = 2 + (iff - (isl - ilo) - 1)/lratiox
eta1 = (-0.5d0+dble(mod(iff-1,lratiox)))/dble(lratiox)
do jf = 1,ncolp
jc = 2 + (jf - (jsb - jlo) - 1)/lratioy
eta2 = (-0.5d0+dble(mod(jf -1,lratioy)))/dble(lratioy)
flag = flaguse(iff,jf)
if (flag .eq. 0.0) then
if (.not.(fineflag(ivalc(ic,jc,1)))) then
* !this is a normal wet cell. intepolate normally
hvf = valcrse(ivalc(ic,jc,ivar))
& + eta1*slopex(icrse(ic,jc))
& + eta2*slopey(icrse(ic,jc))
vf = hvf/valbig(iff+nrowst-1,jf+ncolst-1,1)
if (vf.lt.velmin(icrse(ic,jc)).or.
& vf.gt.velmax(icrse(ic,jc))) then
fineflag(ivalc(ic,jc,ivar))=.true.
reloop = .true.
else
valbig(iff+nrowst-1,jf+ncolst-1,ivar) = hvf
endif
endif
endif
enddo
enddo
* !loop again and reset momentum to conserve momentum
* !in the case of gained momentum, or if velocity bounds violated
if (reloop) then
do iff = 1,nrowp
ic = 2 + (iff - (isl - ilo) - 1)/lratiox
eta1 =(-0.5d0+dble(mod(iff-1,lratiox)))/dble(lratiox)
do jf = 1,ncolp
jc = 2 + (jf - (jsb - jlo) - 1)/lratioy
eta2 = (-0.5d0+dble(mod(jf -1,lratioy)))/dble(lratioy)
flag = flaguse(iff,jf)
if (flag.eq.0.0) then
if (fineflag(ivalc(ic,jc,1))
& .or.fineflag(ivalc(ic,jc,ivar))) then
if (finemass(icrse(ic,jc)).gt.toldry) then
hcrse = valcrse(ivalc(ic,jc,1))
hcnt = dble(icount(icrse(ic,jc)))
hfineave = finemass(icrse(ic,jc))/hcnt
dividemass = max(hcrse,hfineave)
hfine = valbig(iff+nrowst-1,jf+ncolst-1,1)
Vnew = valcrse(ivalc(ic,jc,ivar))/dividemass
valbig(iff+nrowst-1,jf+ncolst-1,ivar) =
& Vnew*valbig(iff+nrowst-1,jf+ncolst-1,1)
else
valbig(iff+nrowst-1,jf+ncolst-1,ivar)=0.d0
endif
endif
endif
enddo
enddo
endif
enddo
c call reclam(loccrse,ntot)
90 continue
c
c set bcs, whether or not recursive calls needed. set any part of patch that stuck out
c
call bc2amr(valbig,aux,mitot,mjtot,nvar,naux,
1 hxf,hyf,level,time,
2 xlp,xrp,ybp,ytp,
3 xlower,ylower,xupper,yupper,
4 xperdom,yperdom,spheredom)
c call reclam(locuse,nrowp*ncolp)
return
end