c c c c ===================================================== subroutine qinit(maxmx,maxmy,meqn,mbc,mx,my,xlower,ylower, & dx,dy,q,maux,aux) c ===================================================== c c # Set initial conditions for q. c implicit double precision (a-h,o-z) dimension q(1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc, meqn) common /comic/ hl,hr,hul,hur c do 20 i=1,mx xc = xlower + (i-0.5d0)*dx xclow = xlower + (i-1.0d0)*dx do 20 j=1,my yc = ylower + (j-0.5d0)*dy yclow = ylower + (j-1.0d0)*dy c # map the center of this computational cell to physical c # coordinates before evaluating the initial value funcion: call mapc2p(xc,yc,xp,yp) c q(i,j,1) = 1.d0 + c & 3.d0*dexp(-50.d0*((xp-1.0d0)**2 + (yp-1.0d0)**2)) c & dexp(-100.d0*((xp-0.8d0)**2 + (yp-1.0d0)**2)) call cellave(xclow,yclow,dx,dy,win) q(i,j,1) = hl*win + hr*(1.d0-win) q(i,j,2) = hul*win + hur*(1.d0-win) q(i,j,3) = 0.d0 20 continue return end c c c ===================================================== subroutine rpn2(ixy,maxm,meqn,mwaves,mbc,mx,ql,qr,auxl,auxr, & wave,s,amdq,apdq) c ===================================================== c c # Roe-solver for the 2D shallow water equations c # on a quadrilateral grid. The velocity components of data are c # rotated to the normal direction and then the standard Roe c # solver and entropy fix are applied, and rotation back is applied c # at the end. c c # solve Riemann problems along one slice of data. c c # On input, ql contains the state vector at the left edge of each cell c # qr contains the state vector at the right edge of each cell c c # This data is along a slice in the x-direction if ixy=1 c # or the y-direction if ixy=2. c # On output, wave contains the waves, s the speeds, c # and amdq, apdq the decomposition of the flux difference c # f(qr(i-1)) - f(ql(i)) c # into leftgoing and rightgoing parts respectively. c # With the Roe solver we have c # amdq = A^- \Delta q and apdq = A^+ \Delta q c # where A is the Roe matrix. An entropy fix can also be incorporated c # into the flux differences. c c # Note that the i'th Riemann problem has left state qr(i-1,:) c # and right state ql(i,:) c # From the basic clawpack routines, this routine is called with ql = qr c c implicit double precision (a-h,o-z) c dimension wave(1-mbc:maxm+mbc, meqn, mwaves) dimension s(1-mbc:maxm+mbc, mwaves) dimension ql(1-mbc:maxm+mbc, meqn) dimension qr(1-mbc:maxm+mbc, meqn) dimension apdq(1-mbc:maxm+mbc, meqn) dimension amdq(1-mbc:maxm+mbc, meqn) dimension auxl(1-mbc:maxm+mbc, 7) dimension auxr(1-mbc:maxm+mbc, 7) c c local arrays -- common block comroe is passed to rpt2sh c ------------ parameter (maxm2 = 1002) !# assumes at most 1000x1000 grid with mbc=2 dimension delta(3) logical efix dimension unorl(-1:maxm2), unorr(-1:maxm2) dimension utanl(-1:maxm2), utanr(-1:maxm2) dimension ax(-1:maxm2) dimension ay(-1:maxm2) common /sw/ g common /comroe/ u(-1:maxm2),v(-1:maxm2),a(-1:maxm2),h(-1:maxm2) c data efix /.true./ !# use entropy fix for transonic rarefactions c if (-1.gt.1-mbc .or. maxm2 .lt. maxm+mbc) then write(6,*) 'need to increase maxm2 in rpA' stop endif c c c # rotate the velocities q(2) and q(3) so that it is aligned with grid c # normal. The normal vector for the face at the i'th Riemann problem c # is stored in the aux array in locations (1,2) if ixy=1 or c # (4,5) if ixy=2. This normal is called (ax,ay) below. c # The ratio of the length of the cell side to the length of the c # computational cell is stored in aux(3) or aux(6) respectively. c c if (ixy.eq.1) then inx = 1 iny = 2 ilenrat = 3 else inx = 4 iny = 5 ilenrat = 6 endif c c # determine rotation matrix c [ ax ay ] c [-ay ax ] c c # note that this reduces to identity on standard cartesian grid c do i=2-mbc,mx+mbc ax(i) = auxl(i,inx) ay(i) = auxl(i,iny) unorl(i) = ax(i)*ql(i,2) + ay(i)*ql(i,3) unorr(i-1) = ax(i)*qr(i-1,2) + ay(i)*qr(i-1,3) utanl(i) = -ay(i)*ql(i,2) + ax(i)*ql(i,3) utanr(i-1) = -ay(i)*qr(i-1,2) + ax(i)*qr(i-1,3) enddo c c c # compute the Roe-averaged variables needed in the Roe solver. c # These are stored in the common block comroe since they are c # later used in routine rpt2 to do the transverse wave splitting. c do 10 i = 2-mbc, mx+mbc h(i) = (qr(i-1,1)+ql(i,1))*0.50d0 hsqrtl = dsqrt(qr(i-1,1)) hsqrtr = dsqrt(ql(i,1)) hsq2 = hsqrtl + hsqrtr u(i) = (unorr(i-1)/hsqrtl + unorl(i)/hsqrtr) / hsq2 v(i) = (utanr(i-1)/hsqrtl + utanl(i)/hsqrtr) / hsq2 a(i) = dsqrt(g*h(i)) 10 continue c c c # now split the jump in q at each interface into waves c c # find a1 thru a3, the coefficients of the 3 eigenvectors: do 20 i = 2-mbc, mx+mbc delta(1) = ql(i,1) - qr(i-1,1) delta(2) = unorl(i) - unorr(i-1) delta(3) = utanl(i) - utanr(i-1) a1 = ((u(i)+a(i))*delta(1) - delta(2))*(0.50d0/a(i)) a2 = -v(i)*delta(1) + delta(3) a3 = (-(u(i)-a(i))*delta(1) + delta(2))*(0.50d0/a(i)) c c # Compute the waves. c wave(i,1,1) = a1 wave(i,2,1) = ax(i)*a1*(u(i)-a(i)) - ay(i)*a1*v(i) wave(i,3,1) = ay(i)*a1*(u(i)-a(i)) + ax(i)*a1*v(i) s(i,1) = (u(i)-a(i)) * auxl(i,ilenrat) c wave(i,1,2) = 0.0d0 wave(i,2,2) = -ay(i)*a2 wave(i,3,2) = ax(i)*a2 s(i,2) = u(i) * auxl(i,ilenrat) c wave(i,1,3) = a3 wave(i,2,3) = ax(i)*a3*(u(i)+a(i)) - ay(i)*a3*v(i) wave(i,3,3) = ay(i)*a3*(u(i)+a(i)) + ax(i)*a3*v(i) s(i,3) = (u(i)+a(i)) * auxl(i,ilenrat) 20 continue c c c # compute flux differences amdq and apdq. c --------------------------------------- c if (efix) go to 110 c c # no entropy fix c ---------------- c c # amdq = SUM s*wave over left-going waves c # apdq = SUM s*wave over right-going waves c do 100 m=1,3 do 100 i=2-mbc, mx+mbc amdq(i,m) = 0.d0 apdq(i,m) = 0.d0 do 90 mw=1,mwaves if (s(i,mw) .lt. 0.d0) then amdq(i,m) = amdq(i,m) + s(i,mw)*wave(i,m,mw) else apdq(i,m) = apdq(i,m) + s(i,mw)*wave(i,m,mw) endif 90 continue 100 continue go to 900 c c----------------------------------------------------- c 110 continue c c # With entropy fix c ------------------ c c # compute flux differences amdq and apdq. c # First compute amdq as sum of s*wave for left going waves. c # Incorporate entropy fix by adding a modified fraction of wave c # if s should change sign. c do 200 i=2-mbc,mx+mbc c check 1-wave him1 = qr(i-1,1) s0 = (unorr(i-1)/him1 - dsqrt(g*him1)) * auxl(i,ilenrat) c check for fully supersonic case : if (s0.gt.0.0d0.and.s(i,1).gt.0.0d0) then do 60 m=1,3 amdq(i,m)=0.0d0 60 continue goto 200 endif c h1 = qr(i-1,1)+wave(i,1,1) hu1= unorr(i-1)+ ax(i)*wave(i,2,1) + ay(i)*wave(i,3,1) s1 = (hu1/h1 - dsqrt(g*h1))* auxl(i,ilenrat) !speed just to right of 1-wave if (s0.lt.0.0d0.and.s1.gt.0.0d0) then c transonic rarefaction in 1-wave sfract = s0*((s1-s(i,1))/(s1-s0)) else if (s(i,1).lt.0.0d0) then c 1-wave is leftgoing sfract = s(i,1) else c 1-wave is rightgoing sfract = 0.0d0 endif do 120 m=1,3 amdq(i,m) = sfract*wave(i,m,1) 120 continue c check 2-wave if (s(i,2).gt.0.0d0) then c #2 and 3 waves are right-going go to 200 endif do 140 m=1,3 amdq(i,m) = amdq(i,m) + s(i,2)*wave(i,m,2) 140 continue c c check 3-wave c hi = ql(i,1) s03 = (unorl(i)/hi + dsqrt(g*hi)) * auxl(i,ilenrat) h3=ql(i,1)-wave(i,1,3) hu3=unorl(i)- (ax(i)*wave(i,2,3) + ay(i)*wave(i,3,3)) s3=(hu3/h3 + dsqrt(g*h3)) * auxl(i,ilenrat) if (s3.lt.0.0d0.and.s03.gt.0.0d0) then c transonic rarefaction in 3-wave sfract = s3*((s03-s(i,3))/(s03-s3)) else if (s(i,3).lt.0.0d0) then c 3-wave is leftgoing sfract = s(i,3) else c 3-wave is rightgoing goto 200 endif do 160 m=1,3 amdq(i,m) = amdq(i,m) + sfract*wave(i,m,3) 160 continue 200 continue c c compute rightgoing flux differences : c do 220 m=1,3 do 220 i = 2-mbc,mx+mbc df = 0.0d0 do 210 mw=1,mwaves df = df + s(i,mw)*wave(i,m,mw) 210 continue apdq(i,m)=df-amdq(i,m) 220 continue c c 900 continue return end c c c ===================================================== subroutine rpt2(ixy,maxm,meqn,mwaves,mbc,mx, & ql,qr,aux1,aux2,aux3, & imp,asdq,bmasdq,bpasdq) c ===================================================== implicit double precision (a-h,o-z) c c # Riemann solver in the transverse direction for the shallow water c # equations on a quadrilateral grid. c c # Split asdq (= A^* \Delta q, where * = + or -) c # into down-going flux difference bmasdq (= B^- A^* \Delta q) c # and up-going flux difference bpasdq (= B^+ A^* \Delta q) c c # Uses Roe averages and other quantities which were c # computed in rpn2sh and stored in the common block comroe. c dimension ql(1-mbc:maxm+mbc, meqn) dimension qr(1-mbc:maxm+mbc, meqn) dimension asdq(1-mbc:maxm+mbc, meqn) dimension bmasdq(1-mbc:maxm+mbc, meqn) dimension bpasdq(1-mbc:maxm+mbc, meqn) dimension aux1(1-mbc:maxm+mbc, 7) dimension aux2(1-mbc:maxm+mbc, 7) dimension aux3(1-mbc:maxm+mbc, 7) c parameter (maxm2 = 1002) !# assumes at most 1000x1000 grid with mbc=2 dimension u(-1:maxm2),v(-1:maxm2),a(-1:maxm2),h(-1:maxm2) dimension ax(-1:maxm2) dimension ay(-1:maxm2) dimension wave(-1:maxm2, 3, 3) dimension s(-1:maxm2, 3) dimension delta(3) common /sw/ g c if (-1.gt.1-mbc .or. maxm2 .lt. maxm+mbc) then write(6,*) 'need to increase maxm2 in rpB' stop endif c c if (ixy.eq.1) then inx = 4 iny = 5 ilenrat = 6 else inx = 1 iny = 2 ilenrat = 3 endif c c # imp is used to flag whether wave is going to left or right, c # since states and grid are different on each side c if (imp.eq.1) then c # asdq = amdq, moving to left ix1 = 2-mbc ixm1 = mx+mbc else c # asdq = apdq, moving to right ix1 = 1-mbc ixm1 = mx+mbc endif c c -------------- c # up-going: c -------------- c c # determine rotation matrix for interface above cell, using aux3 c [ ax ay ] c [-ay ax ] c do i=ix1,ixm1 c if (imp.eq.1) then i1 = i-1 else i1 = i endif c ax(i) = aux3(i1,inx) ay(i) = aux3(i1,iny) h(i) = ql(i1,1) u(i) = (ax(i)*ql(i1,2) + ay(i)*ql(i1,3)) / h(i) v(i) = (-ay(i)*ql(i1,2) + ax(i)*ql(i1,3)) / h(i) a(i) = dsqrt(g*h(i)) enddo c c c c # now split asdq into waves: c c # find a1 thru a3, the coefficients of the 3 eigenvectors: do 20 i = ix1,ixm1 delta(1) = asdq(i,1) delta(2) = ax(i)*asdq(i,2) + ay(i)*asdq(i,3) delta(3) = -ay(i)*asdq(i,2) + ax(i)*asdq(i,3) a1 = ((u(i)+a(i))*delta(1) - delta(2))*(0.50d0/a(i)) a2 = -v(i)*delta(1) + delta(3) a3 = (-(u(i)-a(i))*delta(1) + delta(2))*(0.50d0/a(i)) c c # Compute the waves. c wave(i,1,1) = a1 wave(i,2,1) = ax(i)*a1*(u(i)-a(i)) - ay(i)*a1*v(i) wave(i,3,1) = ay(i)*a1*(u(i)-a(i)) + ax(i)*a1*v(i) s(i,1) = (u(i)-a(i)) * aux3(i1,ilenrat) c wave(i,1,2) = 0.0d0 wave(i,2,2) = -ay(i)*a2 wave(i,3,2) = ax(i)*a2 s(i,2) = u(i) * aux3(i1,ilenrat) c wave(i,1,3) = a3 wave(i,2,3) = ax(i)*a3*(u(i)+a(i)) - ay(i)*a3*v(i) wave(i,3,3) = ay(i)*a3*(u(i)+a(i)) + ax(i)*a3*v(i) s(i,3) = (u(i)+a(i)) * aux3(i1,ilenrat) 20 continue c c c # compute flux difference bpasdq c -------------------------------- c do 40 m=1,3 do 40 i=ix1,ixm1 bpasdq(i,m) = 0.d0 do 30 mw=1,mwaves bpasdq(i,m) = bpasdq(i,m) & + dmax1(s(i,mw),0.d0)*wave(i,m,mw) 30 continue 40 continue c c c -------------- c # down-going: c -------------- c c # determine rotation matrix for interface below cell, using aux2 c [ ax ay ] c [-ay ax ] c do i=ix1,ixm1 c if (imp.eq.1) then i1 = i-1 else i1 = i endif c ax(i) = aux2(i1,inx) ay(i) = aux2(i1,iny) u(i) = (ax(i)*ql(i1,2) + ay(i)*ql(i1,3)) / h(i) v(i) = (-ay(i)*ql(i1,2) + ax(i)*ql(i1,3)) / h(i) enddo c c c c # now split asdq into waves: c c # find a1 thru a3, the coefficients of the 3 eigenvectors: do 80 i = ix1,ixm1 delta(1) = asdq(i,1) delta(2) = ax(i)*asdq(i,2) + ay(i)*asdq(i,3) delta(3) = -ay(i)*asdq(i,2) + ax(i)*asdq(i,3) a1 = ((u(i)+a(i))*delta(1) - delta(2))*(0.50d0/a(i)) a2 = -v(i)*delta(1) + delta(3) a3 = (-(u(i)-a(i))*delta(1) + delta(2))*(0.50d0/a(i)) c c # Compute the waves. c wave(i,1,1) = a1 wave(i,2,1) = ax(i)*a1*(u(i)-a(i)) - ay(i)*a1*v(i) wave(i,3,1) = ay(i)*a1*(u(i)-a(i)) + ax(i)*a1*v(i) s(i,1) = (u(i)-a(i)) * aux2(i1,ilenrat) c wave(i,1,2) = 0.0d0 wave(i,2,2) = -ay(i)*a2 wave(i,3,2) = ax(i)*a2 s(i,2) = u(i) * aux2(i1,ilenrat) c wave(i,1,3) = a3 wave(i,2,3) = ax(i)*a3*(u(i)+a(i)) - ay(i)*a3*v(i) wave(i,3,3) = ay(i)*a3*(u(i)+a(i)) + ax(i)*a3*v(i) s(i,3) = (u(i)+a(i)) * aux2(i1,ilenrat) 80 continue c c c # compute flux difference bmasdq c -------------------------------- c do 100 m=1,3 do 100 i=ix1,ixm1 bmasdq(i,m) = 0.d0 do 90 mw=1,mwaves bmasdq(i,m) = bmasdq(i,m) & + dmin1(s(i,mw), 0.d0)*wave(i,m,mw) 90 continue 100 continue c c return end subroutine setprob implicit double precision (a-h,o-z) c c # Copy this file to your directory and modify to set up problem c # parameters or read other data. c common /sw/ g common/cdisc/ x0,y0,alf,beta,r0,idisc common /comic/ hl,hr,hul,hur g = 1.d0 c # data for flow into cylinder: idisc = 1 x0 = -2.d0 y0 = 0.d0 alf = 1.d0 beta = 0.d0 hl = 4.d0 hr = 1.d0 hur = 0.d0 hul = hur + (hl-hr)*(hur/hr + dsqrt(g*hr + 0.5d0*g*(hl-hr)* & (3.d0 + (hl-hr)/hr))) return end c ============================================ subroutine setaux(maxmx,maxmy,mbc,mx,my,xlower,ylower,dxc,dyc, & maux,aux) c ============================================ c c c # aux(i,j,1) = ax c # aux(i,j,2) = ay where (ax,ay) is unit normal to left face c # aux(i,j,3) = ratio of length of left face to dyc c c # aux(i,j,4) = bx c # aux(i,j,5) = by where (bx,by) is unit normal to bottom face c # aux(i,j,6) = ratio of length of bottom face to dxc c c # aux(i,j,7) = ratio of cell area to dxc*dyc c # (approximately Jacobian of mapping function) c c implicit double precision (a-h,o-z) dimension aux(1-mbc:maxmx+mbc,1-mbc:maxmy+mbc, 7) dimension xccorn(5),yccorn(5),xpcorn(5),ypcorn(5) c dx2 = dxc/2.d0 dy2 = dyc/2.d0 c do 20 j=1-mbc,my+mbc do 20 i=1-mbc,mx+mbc c c # computational points (xc,yc) are mapped to physical c # coordinates (xp,yp) by mapc2p: c c # lower left corner: xccorn(1) = xlower + (i-1)*dxc yccorn(1) = ylower + (j-1)*dyc call mapc2p(xccorn(1),yccorn(1),xpcorn(1),ypcorn(1)) c # upper left corner: xccorn(2) = xccorn(1) yccorn(2) = yccorn(1) + dyc call mapc2p(xccorn(2),yccorn(2),xpcorn(2),ypcorn(2)) c c # upper right corner: xccorn(3) = xccorn(1) + dxc yccorn(3) = yccorn(1) + dyc call mapc2p(xccorn(3),yccorn(3),xpcorn(3),ypcorn(3)) c c # lower right corner: xccorn(4) = xccorn(1) + dxc yccorn(4) = yccorn(1) call mapc2p(xccorn(4),yccorn(4),xpcorn(4),ypcorn(4)) c c # compute normals to left and bottom side: c ax = (ypcorn(2) - ypcorn(1)) ay = -(xpcorn(2) - xpcorn(1)) anorm = dsqrt(ax*ax + ay*ay) aux(i,j,1) = ax/anorm aux(i,j,2) = ay/anorm aux(i,j,3) = anorm/dyc c bx = -(ypcorn(4) - ypcorn(1)) by = (xpcorn(4) - xpcorn(1)) bnorm = dsqrt(bx*bx + by*by) aux(i,j,4) = bx/bnorm aux(i,j,5) = by/bnorm aux(i,j,6) = bnorm/dxc c c # compute area of physical cell from four corners: xpcorn(5) = xpcorn(1) ypcorn(5) = ypcorn(1) area = 0.d0 do ic=1,4 area = area + 0.5d0 * (ypcorn(ic)+ypcorn(ic+1)) * & (xpcorn(ic+1)-xpcorn(ic)) enddo aux(i,j,7) = area / (dxc*dyc) c 20 continue c return end c c c c ================================================= function fdisc(xc,yc) c ================================================= implicit double precision (a-h,o-z) common/cdisc/ x0,y0,alf,beta,r0,idisc c c # for computing cell averages for initial data that has a c # discontinuity along some curve. fdisc should be negative to the c # left of the curve and positive to the right c # idisc specifies the nature of the discontinuity for two c # particular cases (a straight line and circle) but this routine c # can be modified for any other curve. c c # map to quadrilateral grid: call mapc2p(xc,yc,x,y) go to (10,20) idisc c 10 continue c # straight line through (x0,y0) with normal (alf,beta) pointing c # into right state c fdisc = (x-x0)*alf + (y-y0)*beta return c 20 continue c # circle of radius r0: fdisc = (x-x0)**2 + (y-y0)**2 - r0**2 c return end c c ------------------------------------------------------------------ c subroutine bc2amr(val,aux,nrow,ncol,meqn,naux, 1 hx, hy, level, time, 2 xleft, xright, ybot, ytop, 3 xlower, ylower,xupper,yupper, 4 xperiodic, yperiodic) c c c :::::::::: bc2amr ::::::::::::::::::::::::::::::::::::::::::::::; c c Take a grid patch with mesh widths hx,hy, of dimensions nrow by c ncol, and set the values of any piece of c of the patch which extends outside the physical domain c using the boundary conditions. c c ------------------------------------------------ c # Standard boundary condition choices for amr2ez in clawpack c # modified for a general quadrilateral grid in the case mthbc(k)=3 c c # At each boundary k = 1 (left), 2 (right), 3 (top), 4 (bottom): c # mthbc(k) = 0 for user-supplied BC's (must be inserted!) c # = 1 for zero-order extrapolation c # = 2 for periodic boundary coniditions c # = 3 for solid walls, assuming this can be implemented c # by reflecting the data about the boundary and then c # negating the normal component of the velocity. c # On a quadrilateral grid we know the normal to each c # edge and assume this is stored in the aux array: c c # aux(i,j,1) = ax c # aux(i,j,2) = ay where (ax,ay) is unit normal to left face c # aux(i,j,4) = bx c # aux(i,j,5) = by where (bx,by) is unit normal to bottom face c ------------------------------------------------ c c The corners of the grid patch are at c (xleft,ybot) -- lower left corner c (xright,ytop) -- upper right corner c c The physical domain itself is a rectangle bounded by c (xlower,ylower) -- lower left corner c (xupper,yupper) -- upper right corner c c the picture is the following: c c _____________________ (xupper,yupper) c | | c _________ (xright,ytop) | c | | | | c | | | | c | | | | c |___|____| | c (xleft,ybot) | | c | | c |_____________________| c (xlower,ylower) c c c Any cells that lie outside the physical domain are ghost cells whose c values should be set in this routine. This is tested for by comparing c xleft with xlower to see if values need to be set at the left, as in c the figure above, and similarly at the other boundaries. c c Patches are guaranteed to have at least 1 row of cells filled c with interior values so it is possible to extrapolate. c Fix trimbd if you want more than 1 row pre-set. c c Make sure the order the boundaries are specified is correct c so that diagonal corner cells are also properly taken care of. c c Periodic boundaries are set before calling this routine, so if the c domain is periodic in one direction only you c can safely extrapolate in the other direction. c c Don't overwrite ghost cells in periodic directions! c c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::; implicit double precision (a-h,o-z) common /combc2/ mthbc(4) dimension val(nrow,ncol,meqn), aux(nrow,ncol,naux) logical xperiodic, yperiodic hxmarg = hx*.01 hymarg = hy*.01 if (xperiodic .and. yperiodic) go to 499 c c c------------------------------------------------------- c # left boundary: c------------------------------------------------------- if (xleft .ge. xlower-hxmarg) then c # not a physical boundary -- ghost cells lie within another c # grid and values are set elsewhere in amr code. go to 199 endif c c # number of ghost cells lying outside physical domain: nxl = (xlower+hxmarg-xleft)/hx c go to (100,110,120,130) mthbc(1)+1 c 100 continue c # user-specified boundary conditions go here in place of error output do 105 i=1,nxl do 105 j = 1,ncol val(i,j,1) = rhol val(i,j,2) = rhoul val(i,j,3) = 0.d0 val(i,j,4) = el 105 continue go to 199 c 110 continue c # zero-order extrapolation: do 115 m=1,meqn do 115 i=1,nxl do 115 j = 1,ncol val(i,j,m) = val(nxl+1,j,m) 115 continue go to 199 120 continue c # periodic: handled elsewhere in amr go to 199 130 continue c # solid wall (assumes 2'nd component is velocity or momentum in x): do 135 m=1,meqn do 135 i=1,nxl do 135 j = 1,ncol val(i,j,m) = val(2*nxl+1-i,j,m) 135 continue c c # negate the normal velocity: c # (for a general quadrilateral grid) c do 136 i=1,nxl do 136 j = 1,ncol alf = aux(nxl+1,j,1) beta = aux(nxl+1,j,2) unorm = alf*val(2*nxl+1-i,j,2) + beta*val(2*nxl+1-i,j,3) utang = -beta*val(2*nxl+1-i,j,2) + alf*val(2*nxl+1-i,j,3) unorm = -unorm val(i,j,2) = alf*unorm - beta*utang val(i,j,3) = beta*unorm + alf*utang 136 continue go to 199 199 continue c c------------------------------------------------------- c # right boundary: c------------------------------------------------------- if (xright .le. xupper+hxmarg) then c # not a physical boundary -- ghost cells lie within another c # grid and values are set elsewhere in amr code. go to 299 endif c c # number of ghost cells lying outside physical domain: nxr = (xright - xupper + hxmarg)/hx ibeg = max0(nrow-nxr+1, 1) c go to (200,210,220,230) mthbc(2)+1 c 200 continue c # user-specified boundary conditions go here in place of error output write(6,*) & '*** ERROR *** mthbc(2)=0 and no BCs specified in bc2amr' stop go to 299 210 continue c # zero-order extrapolation: do 215 m=1,meqn do 215 i=ibeg,nrow do 215 j = 1,ncol val(i,j,m) = val(ibeg-1,j,m) 215 continue go to 299 220 continue c # periodic: handled elsewhere in amr go to 299 230 continue c # solid wall (assumes 2'nd component is velocity or momentum in x): do 235 m=1,meqn do 235 i=ibeg,nrow do 235 j = 1,ncol val(i,j,m) = val(2*ibeg-1-i,j,m) 235 continue c c # negate the normal velocity: c # (for a general quadrilateral grid) c do 236 i=ibeg,nrow do 236 j = 1,ncol alf = aux(ibeg,j,1) beta = aux(ibeg,j,2) unorm = alf*val(2*ibeg-1-i,j,2) + beta*val(2*ibeg-1-i,j,3) utang = -beta*val(2*ibeg-1-i,j,2) + alf*val(2*ibeg-1-i,j,3) unorm = -unorm val(i,j,2) = alf*unorm - beta*utang val(i,j,3) = beta*unorm + alf*utang 236 continue go to 299 299 continue c c------------------------------------------------------- c # bottom boundary: c------------------------------------------------------- if (ybot .ge. ylower-hymarg) then c # not a physical boundary -- ghost cells lie within another c # grid and values are set elsewhere in amr code. go to 399 endif c c # number of ghost cells lying outside physical domain: nyb = (ylower+hymarg-ybot)/hy c go to (300,310,320,330) mthbc(3)+1 c 300 continue c # user-specified boundary conditions go here in place of error output write(6,*) & '*** ERROR *** mthbc(3)=0 and no BCs specified in bc2amr' stop go to 399 c 310 continue c # zero-order extrapolation: do 315 m=1,meqn do 315 j=1,nyb do 315 i=1,nrow val(i,j,m) = val(i,nyb+1,m) 315 continue go to 399 320 continue c # periodic: handled elsewhere in amr go to 399 330 continue c # solid wall (assumes 3'rd component is velocity or momentum in y): do 335 m=1,meqn do 335 j=1,nyb do 335 i=1,nrow val(i,j,m) = val(i,2*nyb+1-j,m) 335 continue c c # negate the normal velocity: c # (for a general quadrilateral grid) c do 336 j=1,nyb do 336 i=1,nrow alf = aux(i,nyb+1,4) beta = aux(i,nyb+1,5) unorm = alf*val(i,2*nyb+1-j,2) + beta*val(i,2*nyb+1-j,3) utang = -beta*val(i,2*nyb+1-j,2) + alf*val(i,2*nyb+1-j,3) unorm = -unorm val(i,j,2) = alf*unorm - beta*utang val(i,j,3) = beta*unorm + alf*utang 336 continue go to 399 399 continue c c------------------------------------------------------- c # top boundary: c------------------------------------------------------- if (ytop .le. yupper+hymarg) then c # not a physical boundary -- ghost cells lie within another c # grid and values are set elsewhere in amr code. go to 499 endif c c # number of ghost cells lying outside physical domain: nyt = (ytop - yupper + hymarg)/hy jbeg = max0(ncol-nyt+1, 1) c go to (400,410,420,430) mthbc(4)+1 c 400 continue c # user-specified boundary conditions go here in place of error output write(6,*) & '*** ERROR *** mthbc(4)=0 and no BCs specified in bc2amr' stop go to 499 410 continue c # zero-order extrapolation: do 415 m=1,meqn do 415 j=jbeg,ncol do 415 i=1,nrow val(i,j,m) = val(i,jbeg-1,m) 415 continue go to 499 420 continue c # periodic: handled elsewhere in amr go to 499 430 continue c # solid wall (assumes 3'rd component is velocity or momentum in y): do 435 m=1,meqn do 435 j=jbeg,ncol do 435 i=1,nrow val(i,j,m) = val(i,2*jbeg-1-j,m) 435 continue c c # negate the normal velocity: c # (for a general quadrilateral grid) c do 436 j=jbeg,ncol do 436 i=1,nrow alf = aux(i,jbeg,4) beta = aux(i,jbeg,5) unorm = alf*val(i,2*jbeg-1-j,2) + beta*val(i,2*jbeg-1-j,3) utang = -beta*val(i,2*jbeg-1-j,2) + alf*val(i,2*jbeg-1-j,3) unorm = -unorm val(i,j,2) = alf*unorm - beta*utang val(i,j,3) = beta*unorm + alf*utang 436 continue go to 499 499 continue return end c c ===================================================== subroutine mapc2p(xc,yc,xp,yp) c ===================================================== c c # on input, (xc,yc) is a computational grid point c # on output, (xp,yp) is corresponding point in physical space c implicit double precision (a-h,o-z) c c # squared off radial coordinates, xc = fraction of radius, yc = theta c r0 = 1.d0 r1 = 5.d0 pi = 4.d0*datan(1.d0) yc1 = dmin1(yc, dabs(pi/2-yc)) yc1 = dmin1(yc1, dabs(pi-yc)) yc1 = dmin1(yc1, dabs(3*pi/2-yc)) yc1 = dmin1(yc1, dabs(2*pi-yc)) r = r0 + (r1*dsqrt(1.d0 + dtan(yc1)**2) - r0) * xc xp = r * dcos(yc) yp = r * dsin(yc) c return end c c ---------------------------------------------------------------- c program amr2ez c c Use adaptive mesh refinement to solve the hyperbolic 2-d equation: c c u + f(u) + g(u) = 0 c t x y c c or the more general non-conservation law form: c u + A u + B u = 0 c t x y c c using the wave propagation method as in CLAWPACK in combination c with the locally uniform embedded grids of AMR. c Estimate error with Richardson extrap. (in errest.f) c + gradient checking (in errsp.f). Initial conditions set c in (qinit.f), b.c.'s in (physbd.f). c Specify rectangular domain from c (xlower,ylower) to (xupper,yupper). c c No rotated rectangles are used in this version. c Periodic b.c.'s finally implemented. c c ========================================================================= c Copyright 1996, Marsha J. Berger and Randall J. LeVeque c 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 *** 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. c c -------------------------------------- c AMRCLAW Version 0.4, June, 1999 c compatible with CLAWPACK Version 4.0 c Homepage: http://www.amath.washington.edu/~claw/ c -------------------------------------- c c Authors: c c Marsha J. Berger c Courant Institute of Mathematical Sciences c New York University c 251 Mercer St. c New York, NY 10012 c berger@cims.nyu.edu c c Randall J. LeVeque c Applied Mathematics c Box 352420 c University of Washington, c Seattle, WA 98195-2420 c rjl@amath.washington.edu c c ========================================================================= c c c ---------------------------------------------------------------- c implicit double precision (a-h,o-z) include "call.i" common /combc2/ mthbc(4) character * 12 pltfile,infile,outfile,rstfile,dbugfile,matfile character*10 matname2 logical vtime,rest dimension tout(maxout) integer oldmode c c c you may want to turn this on for SUN workstation, or replace c set to signal on overflow, divide by zero, and illegal operation c c oldmode = ieee_handler("set","common",SIGFPE_ABORT) c if (oldmode .ne. 0) then c write(outunit,*)' could not set ieee trapper ' c write(*,*) ' could not set ieee trapper ' c stop c endif c infile = 'amr2ez.data' outfile = 'fort.amr' pltfile = 'fort.ncar' rstfile = 'restart.data' dbugfile = 'fort.debug' matfile = 'fort.nplot' open(inunit, file=infile,status='old',form='formatted') open(outunit, file=outfile,status='unknown',form='formatted') open(dbugunit,file=dbugfile,status='unknown',form='formatted') c c domain variables read(inunit,*) nx read(inunit,*) ny read(inunit,*) mxnest if (mxnest .gt. maxlv) then write(outunit,*) & 'Error *** mxnest > max. allowable levels (maxlv) in common' write(*,*) & 'Error *** mxnest > max. allowable levels (maxlv) in common' stop endif read(inunit,*) (intrat(i),i=1,max(1,mxnest-1)) read(inunit,*) nout if (nout .gt. maxout) then write(outunit,*) 'Error *** nout > maxout in common' write(*,*) 'Error *** nout > maxout in common' stop endif read(inunit,*) outstyle if (outstyle.eq.1) then read(inunit,*) tfinal c # array tout is set below after reading t0 endif if (outstyle.eq.2) then read(inunit,*) (tout(i), i=1,nout) endif if (outstyle.eq.3) then read(inunit,*) iout,nstop nout = 0 endif read(inunit,*) possk(1) read(inunit,*) dtv2 read(inunit,*) cflv1 read(inunit,*) cfl read(inunit,*) nv1 if (outstyle.eq.1 .or. outstyle.eq.2) then nstop = nv1 endif read(inunit,*) method(1) vtime = (method(1) .eq. 1) read(inunit,*) method(2) iorder = method(2) read(inunit,*) method(3) if (method(3) .lt. 0) then write(6,*) '*** ERROR *** method(3) < 0' write(6,*) ' dimensional splitting not supported in amrclaw' stop endif read(inunit,*) method(4) read(inunit,*) method(5) read(inunit,*) mcapa1 read(inunit,*) naux if (naux .gt. maxaux) then write(outunit,*) 'Error *** naux > maxaux in common' write(*,*) 'Error *** naux > maxaux in common' stop endif do iaux = 1, naux read(inunit,*) auxtype(iaux) end do read(inunit,*) nvar read(inunit,*) mwaves if (mwaves .gt. maxwave) then write(outunit,*) 'Error *** mwaves > maxwave in common' write(*,*) 'Error *** mwaves > maxwave in common' stop endif read(inunit,*) (mthlim(mw), mw=1,mwaves) read(inunit,*) t0 read(inunit,*) xlower read(inunit,*) xupper read(inunit,*) ylower read(inunit,*) yupper read(inunit,*) nghost read(inunit,*) mthbc(1) read(inunit,*) mthbc(2) read(inunit,*) mthbc(3) read(inunit,*) mthbc(4) xperdom = (mthbc(1).eq.2 .and. mthbc(2).eq.2) yperdom = (mthbc(3).eq.2 .and. mthbc(4).eq.2) if ((mthbc(1).eq.2 .and. mthbc(2).ne.2) .or. & (mthbc(2).eq.2 .and. mthbc(1).ne.2)) then write(6,*) '*** ERROR *** periodic boundary conditions: ' write(6,*) ' mthbc(1) and mthbc(2) must BOTH be set to 2' stop endif if ((mthbc(3).eq.2 .and. mthbc(4).ne.2) .or. & (mthbc(4).eq.2 .and. mthbc(3).ne.2)) then write(6,*) '*** ERROR *** periodic boundary conditions: ' write(6,*) ' mthbc(3) and mthbc(4) must BOTH be set to 2' stop endif if (outstyle.eq.1) then do i=1,nout tout(i) = t0 + i*(tfinal-t0)/float(nout) enddo endif c restart and checkpointing read(inunit,*) rest read(inunit,*) ichkpt c c refinement variables read(inunit,*) tol read(inunit,*) tolsp read(inunit,*) kcheck read(inunit,*) ibuff read(inunit,*) cut c c style of output c read(inunit,*) printout read(inunit,*) ncarout read(inunit,*) matlabout c c c # read verbose/debugging flags read(inunit,*) dprint read(inunit,*) eprint read(inunit,*) edebug read(inunit,*) gprint read(inunit,*) nprint read(inunit,*) pprint read(inunit,*) rprint read(inunit,*) sprint read(inunit,*) tprint read(inunit,*) uprint c # look for capacity function via auxtypes: mcapa = 0 do 20 iaux = 1, naux if (auxtype(iaux) .eq. "capacity") then if (mcapa .ne. 0) then write(*,*)" only 1 capacity allowed" stop else mcapa = iaux endif endif c # change to Version 4.1 terminology: if (auxtype(iaux) .eq. "leftface") auxtype(iaux) = "xleft" if (auxtype(iaux) .eq. "bottomface") auxtype(iaux) = "yleft" if (.not. (auxtype(iaux) .eq."xleft" .or. . auxtype(iaux) .eq. "yleft".or. . auxtype(iaux) .eq. "capacity".or. . auxtype(iaux) .eq. "center")) then write(*,*)" unknown type for auxiliary variables" write(*,*)" use xleft/yleft/center/capacity" stop endif 20 continue c ::: error checking of input data ::::::::::::::::::::::: if (mcapa .ne. mcapa1) then write(outunit,*) 'Error *** mcapa does not agree with auxtype' write(*,*) 'Error *** mcapa does not agree with auxtype' stop endif if (nvar .gt. maxvar) then write(outunit,*) 'Error *** nvar > maxvar in common' write(*,*) 'Error *** nvar > maxvar in common' stop endif if (2*nghost .gt. min(nx,ny)) then mindim = 2*nghost write(outunit,*) 'Error *** need finer domain >', . mindim, ' cells' write(*,*) 'Error *** need finer domain >', . mindim, ' cells' stop endif if (mcapa .gt. naux) then write(outunit,*) 'Error *** mcapa > naux in input file' write(*,*) 'Error *** mcapa > naux in input file' stop endif if (.not. vtime .and. nout .ne. 0) then write(outunit,*)' cannot specify output times with fixed dt' write(*,*) ' cannot specify output times with fixed dt' stop endif c c if (ncarout) . open(pltunit1,file=pltfile,status='unknown',form='formatted') c c # call user routine to set up problem parameters: call setprob() c matlabu = 0 hxposs(1) = (xupper - xlower) / nx hyposs(1) = (yupper - ylower) / ny cflmax = 0.d0 c c c if (rest) then open(rstunit,file=rstfile,status='old',form='unformatted') rewind rstunit call restrt(nsteps,time,nvar) nstart = nsteps else lentot = 0 lenmax = 0 lendim = 0 rvol = 0.0d0 do 8 i = 1, mxnest 8 rvoll(i) = 0.0d0 evol = 0.0d0 call stst1 call domain (nvar,vtime,nx,ny,naux) call setgrd (nvar,cut,naux) time = 0.0d0 nstart = 0 endif c c print out program parameters for this run c write(outunit,107)tol,tolsp,iorder,kcheck,ibuff,nghost,cut, 1 mxnest,ichkpt,cfl write(outunit,109) xupper,yupper,xlower,ylower,nx,ny, 1 (intrat(i),i=1,mxnest) write(outunit,119) naux write(outunit,129) (iaux,auxtype(iaux),iaux=1,naux) if (mcapa .gt. 0) write(outunit,139) mcapa 107 format(/ * ' amrclaw parameters:',//, * ' error tol ',e12.5,/, * ' spatial error tol ',e12.5,/, * ' order of integrator ',i9,/, * ' error checking interval ',i9,/, * ' buffer zone size ',i9,/, * ' nghost ',i9,/, * ' volume ratio cutoff ',e12.5,/, * ' max. refinement level ',i9,/, * ' user sub. calling times ',i9,/, * ' cfl # (if var. delt) ',e12.5,/) 109 format(' xupper(upper corner) ',e12.5,/, * ' yupper(upper corner) ',e12.5,/, * ' xlower(lower corner) ',e12.5,/, * ' ylower(lower corner) ',e12.5,/, * ' nx = no. cells in x dir.',i9,/, * ' ny = no. cells in y dir.',i9,/, * ' refinement ratios ',6i5,/,/) 119 format(' no. auxiliary vars. ',i9) 129 format(' var ',i5,' of type ', a10) 139 format(' capacity fn. is aux. var',i9) c write(6,*) ' ' write(6,*) 'running amr2ez... ' write(6,*) ' ' call outtre (mstart,printout,nvar,naux) call conck(1,nvar) call valout(1,lfine,time,nvar,naux) c c -------------------------------------------------------- c # tick is the main routine which drives the computation: c -------------------------------------------------------- call tick(nvar,iout,nstart,nstop,cut,vtime,time,ichkpt,naux, & nout,tout) c -------------------------------------------------------- c # Done with computation, cleanup: lentotsave = lentot call cleanup(nvar,naux) if (lentot .ne. 0) then write(outunit,*) lentot," words not accounted for ", & "in memory cleanup" write(*,*) lentot," words not accounted for ", & "in memory cleanup" endif c c report on statistics c open(matunit,file=matfile,status='unknown',form='formatted') write(matunit,*) matlabu-1 write(matunit,*) mxnest write(outunit,*) write(outunit,*) write(outunit,901) lentotsave write(outunit,902) lenmax write(outunit,903) lendim write(outunit,904) rvol do 60 level = 1,mxnest 60 write(outunit,905) level, rvoll(level) write(outunit,906) evol if (evol+rvol .gt. 0.) then ratmet = rvol / (evol+rvol) * 100.0d0 else ratmet = 0.0d0 endif write(outunit,907) ratmet write(outunit,908) cflmax 901 format('current space usage = ',i12) 902 format('maximum space usage = ',i12) 903 format('need space dimension = ',i12,/) 904 format('number of cells advanced for time integration = ',f20.6) 905 format(3x,'# cells advanced on level ',i4,' = ',f20.2) 906 format('number of cells advanced for error estimation = ',f20.6,/) 907 format(' percentage of cells advanced in time = ', f10.2) 908 format(' maximum Courant number seen = ', f10.2) c write(outunit,909) 909 format(//,' ------ end of AMRCLAW integration -------- ') c c stop end c ============================================ subroutine b4step2(maxmx,maxmy,mbc,mx,my,meqn,q, & xlower,ylower,dx,dy,time,dt,maux,aux) c ============================================ c c # called before each call to step c # use to set time-dependent aux arrays or perform other tasks. c c # dummy routine c c implicit double precision (a-h,o-z) dimension q(1-mbc:maxmx+mbc,1-mbc:maxmy+mbc, meqn) c dimension aux(1-mbc:maxmx+mbc,1-mbc:maxmy+mbc, maux) c return end c c ------------------------------------------------------------- c subroutine qad(valbig,mitot,mjtot,nvar, . svdflx,qc1d,lenbc,lratio,hx,hy, . maux,aux,auxc1d,delt,mptr) implicit double precision (a-h, o-z) include "call.i" logical qprint dimension valbig(mitot,mjtot,nvar) dimension qc1d(lenbc,nvar) dimension svdflx(nvar,lenbc) dimension aux(mitot,mjtot,maux) dimension auxc1d(lenbc,maux) c c ::::::::::::::::::::::::::: QAD :::::::::::::::::::::::::::::::::: c solve RP between ghost cell value on fine grid and coarse grid c value that ghost cell overlaps. The resulting fluctuations c are added in to coarse grid value, as a conservation fixup. c Done each fine grid time step. If source terms are present, the c coarse grid value is advanced by source terms each fine time step too. c Side 1 is the left side of the fine grid patch. Then go around clockwise. c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c c # local storage c # note that dimension here are bigger than dimensions used c # in rp2, but shouldn't matter since wave is not used in qad c # and for other arrays it is only the last parameter that is wrong c # ok as long as meqn, mwaves < maxvar parameter (max1dp1 = max1d+1) dimension ql(max1dp1,maxvar), qr(max1dp1,maxvar) dimension wave(max1dp1,maxvar,maxvar), s(max1dp1,maxvar) dimension amdq(max1dp1,maxvar), apdq(max1dp1,maxvar) dimension auxl(max1dp1,maxaux), auxr(max1dp1,maxaux) data qprint/.false./ c c aux is auxiliary array with user parameters needed in Riemann solvers c on fine grid corresponding to valbig c auxc1d is coarse grid stuff from around boundary, same format as qc1d c auxl, auxr are work arrays needed to pass stuff to rpn2 c maux is the number of aux variables, which may be zero. c if (qprint) write(dbugunit,*)" working on grid ",mptr tgrid = rnode(timemult, mptr) nc = mjtot-2*nghost nr = mitot-2*nghost c c-------- c side 1 c-------- c do 10 j = nghost+1, mjtot-nghost if (maux.gt.0) then do 5 ma = 1,maux if (auxtype(ma).eq."xleft") then c # Assuming velocity at left-face, this fix c # preserves conservation in incompressible flow: auxl(j-nghost+1,ma) = aux(nghost+1,j,ma) else c # Normal case -- we set the aux arrays c # from the cell corresponding to q auxl(j-nghost+1,ma) = aux(nghost,j,ma) endif 5 continue endif do 10 ivar = 1, nvar ql(j-nghost+1,ivar) = valbig(nghost,j,ivar) 10 continue lind = 0 index = 0 ncrse = (mjtot-2*nghost)/lratio do 20 jc = 1, ncrse index = index + 1 do 25 l = 1, lratio lind = lind + 1 if (maux.gt.0) then do 24 ma=1,maux auxr(lind,ma) = auxc1d(index,ma) 24 continue endif do 25 ivar = 1, nvar 25 qr(lind,ivar) = qc1d(index,ivar) 20 continue index = ncrse if (qprint) then write(dbugunit,*) 'side 1, ql and qr:' do i=2,nc write(dbugunit,4101) i,qr(i-1,1),ql(i,1) enddo 4101 format(i3,4e16.6) if (maux .gt. 0) then write(dbugunit,*) 'side 1, auxr:' do i=2,nc write(dbugunit,4101) i,(auxr(i-1,ma),ma=1,maux) enddo write(dbugunit,*) 'side 1, auxl:' do i=2,nc write(dbugunit,4101) i,(auxl(i,ma),ma=1,maux) enddo endif endif call rpn2(1,max1dp1-2*nghost,nvar,mwaves,nghost,nc+1-2*nghost, . ql,qr,auxl,auxr,wave,s,amdq,apdq) c c we have the wave. for side 1 add into sdflxm c influx = 0 do 30 j = 1, nc/lratio influx = influx + 1 jfine = (j-1)*lratio do 40 ivar = 1, nvar do 50 l = 1, lratio svdflx(ivar,influx) = svdflx(ivar,influx) . + amdq(jfine+l+1,ivar) * hy * delt . + apdq(jfine+l+1,ivar) * hy * delt 50 continue 40 continue 30 continue c-------- c side 2 c-------- c do 210 i = nghost+1, mitot-nghost if (maux.gt.0) then do 205 ma = 1,maux auxr(i-nghost,ma) = aux(i,mjtot-nghost+1,ma) 205 continue endif do 210 ivar = 1, nvar qr(i-nghost,ivar) = valbig(i,mjtot-nghost+1,ivar) 210 continue lind = 0 ncrse = (mitot-2*nghost)/lratio do 220 ic = 1, ncrse index = index + 1 do 225 l = 1, lratio lind = lind + 1 if (maux.gt.0) then do 224 ma=1,maux if (auxtype(ma).eq."yleft") then c # Assuming velocity at bottom-face, this fix c # preserves conservation in incompressible flow: ifine = (ic-1)*lratio + nghost + l auxl(lind+1,ma) = aux(ifine,mjtot-nghost+1,ma) else auxl(lind+1,ma) = auxc1d(index,ma) endif 224 continue endif do 225 ivar = 1, nvar 225 ql(lind+1,ivar) = qc1d(index,ivar) 220 continue if (qprint) then write(dbugunit,*) 'side 2, ql and qr:' do i=1,nr write(dbugunit,4101) i,ql(i+1,1),qr(i,1) enddo if (maux .gt. 0) then write(dbugunit,*) 'side 2, auxr:' do i = 1, nr write(dbugunit,4101) i, (auxr(i,ma),ma=1,maux) enddo write(dbugunit,*) 'side 2, auxl:' do i = 1, nr write(dbugunit,4101) i, (auxl(i,ma),ma=1,maux) enddo endif endif call rpn2(2,max1dp1-2*nghost,nvar,mwaves,nghost,nr+1-2*nghost, . ql,qr,auxl,auxr,wave,s,amdq,apdq) c c we have the wave. for side 2. add into sdflxp c do 230 i = 1, nr/lratio influx = influx + 1 ifine = (i-1)*lratio do 240 ivar = 1, nvar do 250 l = 1, lratio svdflx(ivar,influx) = svdflx(ivar,influx) . - amdq(ifine+l+1,ivar) * hx * delt . - apdq(ifine+l+1,ivar) * hx * delt 250 continue 240 continue 230 continue c-------- c side 3 c-------- c do 310 j = nghost+1, mjtot-nghost if (maux.gt.0) then do 305 ma = 1,maux auxr(j-nghost,ma) = aux(mitot-nghost+1,j,ma) 305 continue endif do 310 ivar = 1, nvar qr(j-nghost,ivar) = valbig(mitot-nghost+1,j,ivar) 310 continue lind = 0 ncrse = (mjtot-2*nghost)/lratio do 320 jc = 1, ncrse index = index + 1 do 325 l = 1, lratio lind = lind + 1 if (maux.gt.0) then do 324 ma=1,maux if (auxtype(ma).eq."xleft") then c # Assuming velocity at left-face, this fix c # preserves conservation in incompressible flow: jfine = (jc-1)*lratio + nghost + l auxl(lind+1,ma) = aux(mitot-nghost+1,jfine,ma) else auxl(lind+1,ma) = auxc1d(index,ma) endif 324 continue endif do 325 ivar = 1, nvar 325 ql(lind+1,ivar) = qc1d(index,ivar) 320 continue if (qprint) then write(dbugunit,*) 'side 3, ql and qr:' do i=1,nc write(dbugunit,4101) i,ql(i,1),qr(i,1) enddo endif call rpn2(1,max1dp1-2*nghost,nvar,mwaves,nghost,nc+1-2*nghost, . ql,qr,auxl,auxr,wave,s,amdq,apdq) c c we have the wave. for side 3 add into sdflxp C do 330 j = 1, nc/lratio influx = influx + 1 jfine = (j-1)*lratio do 340 ivar = 1, nvar do 350 l = 1, lratio svdflx(ivar,influx) = svdflx(ivar,influx) . - amdq(jfine+l+1,ivar) * hy * delt . - apdq(jfine+l+1,ivar) * hy * delt 350 continue 340 continue 330 continue c-------- c side 4 c-------- c do 410 i = nghost+1, mitot-nghost if (maux.gt.0) then do 405 ma = 1,maux if (auxtype(ma).eq."yleft") then c # Assuming velocity at bottom-face, this fix c # preserves conservation in incompressible flow: auxl(i-nghost+1,ma) = aux(i,nghost+1,ma) else auxl(i-nghost+1,ma) = aux(i,nghost,ma) endif 405 continue endif do 410 ivar = 1, nvar ql(i-nghost+1,ivar) = valbig(i,nghost,ivar) 410 continue lind = 0 ncrse = (mitot-2*nghost)/lratio do 420 ic = 1, ncrse index = index + 1 do 425 l = 1, lratio lind = lind + 1 if (maux.gt.0) then do 424 ma=1,maux auxr(lind,ma) = auxc1d(index,ma) 424 continue endif do 425 ivar = 1, nvar 425 qr(lind,ivar) = qc1d(index,ivar) 420 continue if (qprint) then write(dbugunit,*) 'side 4, ql and qr:' do i=1,nr write(dbugunit,4101) i, ql(i,1),qr(i,1) enddo endif call rpn2(2,max1dp1-2*nghost,nvar,mwaves,nghost,nr+1-2*nghost, . ql,qr,auxl,auxr,wave,s,amdq,apdq) c c we have the wave. for side 4. add into sdflxm c do 430 i = 1, nr/lratio influx = influx + 1 ifine = (i-1)*lratio do 440 ivar = 1, nvar do 450 l = 1, lratio svdflx(ivar,influx) = svdflx(ivar,influx) . + amdq(ifine+l+1,ivar) * hx * delt . + apdq(ifine+l+1,ivar) * hx * delt 450 continue 440 continue 430 continue c # for source terms: if (method(5) .ne. 0) then call src1d(nvar,nghost,lenbc,qc1d,maux,auxc1d,tgrid,delt) endif return end c c c ========================================================= subroutine src2(maxmx,maxmy,meqn,mbc,mx,my,xlower,ylower,dx,dy, & q,maux,aux,t,dt) c ========================================================= implicit double precision (a-h,o-z) dimension q(1-mbc:maxmx+mbc,1-mbc:maxmy+mbc, meqn) dimension aux(1-mbc:maxmx+mbc,1-mbc:maxmy+mbc, maux) c c # dummy source routine... does nothing c return end c c c ========================================================= subroutine src1d(meqn,mbc,mx1d,q1d,maux,aux1d,t,dt) c ========================================================= implicit double precision (a-h,o-z) dimension q1d(mx1d, meqn) dimension aux1d(mx1d, maux) c c # dummy source routine... does nothing c c # This routine should be a simplified version of src2 c # which applies source terms for a 1-d slice of data along the c # edge of a grid. This is called only from qad where the conservative c # fix-up is applied and is used to apply source terms over partial c # time steps to the coarse grid cell values used in solving Riemann c # problems at the interface between coarse and fine grids. c c # If the source terms depend only on q, it should be easy to c # adapt src2 to create this routine, just loop over 1:mx1d. c # If the source terms are more complicated, it may not be easy. c c # The code may work fine without applying source terms in this c # context, so using this dummy routine might be successful even when c # source terms are present. c return end c c -------------------------------------------------------------- c subroutine errf1(rctfine,nvar,rctcrse,mptr,mi2tot,mj2tot, 2 mitot,mjtot,rctflg,sperr) implicit double precision (a-h,o-z) include "call.i" dimension rctfine(mitot,mjtot,nvar) dimension rctcrse(mi2tot,mj2tot,nvar) dimension rctflg(mitot,mjtot,nvar) dimension sperr(mitot,mjtot) logical allowed c c c ::::::::::::::::::::::::::::: ERRF1 :::::::::::::::::::::::::::::::: c c compare error estimates in rctfine, rctcrse. If exceed tol, flag. c We put in a hook which allows us to ignore the error estimates c and not refine if we are outside some prescribed region. c c sperr is the spatial component of the error estimate only c provides simple way for user to specify extra flagging in errsp c c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c c # to allow refinement anywhere: allowed(x,y,level) = .true. c c # This function can be changed to allow refinement only on some portion c # of the grid. c # This is useful if you wish to zoom in on some structure in a c # known location but don't want the same level of refinement elsewhere. c # Points are flagged only if one of the errors is greater than the c # corresponding tolerance. c c # For example, to allow refinement of Level 1 grids everywhere but c # of finer grids only for y >= 0.4: c c allowed(x,y,level) = c & (level.le.1 .or. y.ge.0.4d0) c time = rnode(timemult, mptr) xleft = rnode(cornxlo,mptr) levm = node(nestlevel, mptr) hx = hxposs(levm) ybot = rnode(cornylo,mptr) hy = hyposs(levm) dt = possk(levm) numsp = 0 errmax = 0.0d0 err2 = 0.0d0 c order = dt*float(2**(iorder+1) - 2) order = float(2**(iorder+1) - 2) c if (.not. (edebug)) go to 20 write(outunit,107) mptr 107 format(//,' coarsened grid values for grid ',i4) do 10 jj = nghost+1, mj2tot-nghost j = mj2tot + 1 - jj write(outunit,101) (rctcrse(i,j,1), . i = nghost+1, mi2tot-nghost) 10 continue write(outunit,108) mptr 108 format(//, ' fine grid values for grid ',i4) do 15 jj = nghost+1, mjtot-nghost j = mjtot + 1 - jj write(outunit,101) (rctfine(i,j,1),i=nghost+1,mitot-nghost) 15 continue 101 format(' ',13f6.3) c c zero out the exterior locations so they don't affect err.est. c 20 continue jfine = nghost+1 do 35 j = nghost+1, mj2tot-nghost yofj = ybot + (float(jfine) - .5d0)*hy ifine = nghost+1 c do 30 i = nghost+1, mi2tot-nghost rflag = goodpt xofi = xleft + (float(ifine) - .5d0)*hx term1 = rctfine(ifine,jfine,1) term2 = rctfine(ifine+1,jfine,1) term3 = rctfine(ifine+1,jfine+1,1) term4 = rctfine(ifine,jfine+1,1) c # divide by (aval*order) for relative error aval = (term1+term2+term3+term4)/4.d0 est = dabs((aval-rctcrse(i,j,1))/ order) if (est .gt. errmax) errmax = est err2 = err2 + est*est c write(outunit,102) i,j,est 102 format(' i,j,est ',2i5,e12.5) c rctcrse(i,j,2) = est c if (est .ge. tol .and. allowed(xofi,yofj,levm)) then rflag = badpt endif rctcrse(i,j,1) = rflag ifine = ifine + 2 30 continue jfine = jfine + 2 35 continue c c transfer flagged points on cell centered coarse grid c to cell centered fine grid. count flagged points. c c initialize rctflg to 0.0 (no flags) before flagging c do 40 j = 1, mjtot do 40 i = 1, mitot 40 rctflg(i,j,1) = goodpt c c print out intermediate flagged rctcrse (for debugging) c if (eprint) then err2 = dsqrt(err2/float((mi2tot-2*nghost)*(mj2tot-2*nghost))) write(outunit,103) mptr, levm, errmax, err2 103 format(' grid ',i4,' level ',i4, . ' max. error = ',e15.7,' err2 = ',e15.7) if (edebug) then write(outunit,*) ' flagged points on coarsened grid ', . 'for grid ',mptr do 45 jj = nghost+1, mj2tot-nghost j = mj2tot + 1 - jj write(outunit,106) (nint(rctcrse(i,j,1)), . i=nghost+1,mi2tot-nghost) 106 format(1h ,80i1) 45 continue endif endif c jfine = nghost+1 do 70 j = nghost+1, mj2tot-nghost ifine = nghost+1 do 60 i = nghost+1, mi2tot-nghost if (rctcrse(i,j,1) .eq. goodpt) go to 55 rctflg(ifine,jfine,1) = badpt rctflg(ifine+1,jfine,1) = badpt rctflg(ifine,jfine+1,1) = badpt rctflg(ifine+1,jfine+1,1)= badpt 55 ifine = ifine + 2 60 continue jfine = jfine + 2 70 continue c if (edebug) then write(outunit,*)" spatial error for grid ",mptr do 75 jjfine = nghost+1, mjtot-nghost jfine = mjtot + 1 - jjfine write(outunit,101)(sperr(ifine,jfine), . ifine=nghost+1,mitot-nghost) 75 continue endif do 80 jfine = nghost+1, mjtot-nghost yofj = ybot + (float(jfine) - nghost - .5d0)*hy do 80 ifine = nghost+1, mitot-nghost xofi = xleft + (float(ifine) - nghost - .5d0)*hx if (sperr(ifine,jfine) .gt. tolsp .and. allowed(xofi,yofj,levm)) & then rflag = rctflg(ifine,jfine,1) if (rflag .ne. badpt) then rctflg(ifine,jfine,1) = badpt numsp = numsp + 1 endif endif 80 continue if (eprint) then write(outunit,118) numsp,mptr 118 format( i5,' more pts. flagged for spatial error on grid',i4,/) if (edebug) then do 56 jj = nghost+1, mjtot-nghost j = mjtot + 1 - jj write(outunit,106) & (nint(rctflg(i,j,1)),i=nghost+1,mitot-nghost) 56 continue endif endif return end c c -------------------------------------------------------------- c subroutine advanc (level,nvar,dtlevnew,vtime,naux) c implicit double precision (a-h,o-z) include "call.i" logical vtime c c ::::::::::::::; ADVANC ::::::::::::::::::::::::::::::::::::::::::: c integrate all grids at the input 'level' by one step of its delta(t) c this includes: setting the ghost cells c advancing the solution on the grid c adjusting fluxes for flux conservation step later c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c mptr = lstart(level) hx = hxposs(level) hy = hyposs(level) delt = possk(level) 3 continue nx = node(ndihi,mptr) - node(ndilo,mptr) + 1 ny = node(ndjhi,mptr) - node(ndjlo,mptr) + 1 mitot = nx + 2*nghost mjtot = ny + 2*nghost locnew = node(store1,mptr) locaux = node(storeaux,mptr) time = rnode(timemult,mptr) c call bound(time,nvar,nghost,alloc(locnew),mitot,mjtot,mptr, 1 alloc(locaux),naux) mptr = node(levelptr, mptr) if (mptr .ne. 0) go to 3 c c save coarse level values if there is a finer level for wave fixup if (level+1 .le. mxnest) then if (lstart(level+1) .ne. null) then call saveqc(level+1,nvar,naux) endif endif c dtlevnew = rinfinity mptr = lstart(level) 5 continue locold = node(store2, mptr) locnew = node(store1, mptr) nx = node(ndihi,mptr) - node(ndilo,mptr) + 1 ny = node(ndjhi,mptr) - node(ndjlo,mptr) + 1 time = rnode(timemult,mptr) c mitot = nx + 2*nghost mjtot = ny + 2*nghost c ::: get scratch storage for fluxes and slopes locfp = igetsp(mitot*mjtot*nvar) locfm = igetsp(mitot*mjtot*nvar) locgp = igetsp(mitot*mjtot*nvar) locgm = igetsp(mitot*mjtot*nvar) c c copy old soln. values into next time step's soln. values c since integrator will overwrite it. only for grids not at c the finest level. finest level grids do not maintain copies c of old and new time solution values. c if (level .lt. mxnest) then ntot = mitot * mjtot * nvar cdir$ ivdep do 10 i = 1, ntot 10 alloc(locold + i - 1) = alloc(locnew + i - 1) endif c xlow = rnode(cornxlo,mptr) - nghost*hx ylow = rnode(cornylo,mptr) - nghost*hy rvol = rvol + nx * ny rvoll(level) = rvoll(level) + nx * ny locaux = node(storeaux,mptr) c if (node(ffluxptr,mptr) .ne. 0) then lenbc = 2*(nx/intrat(level-1)+ny/intrat(level-1)) locsvf = node(ffluxptr,mptr) locsvq = locsvf + nvar*lenbc locx1d = locsvq + nvar*lenbc call qad(alloc(locnew),mitot,mjtot,nvar, 1 alloc(locsvf),alloc(locsvq),lenbc, 2 intrat(level-1),hx,hy, 3 naux,alloc(locaux),alloc(locx1d),delt,mptr) endif c call stepgrid(alloc(locnew),alloc(locfm),alloc(locfp), 1 alloc(locgm),alloc(locgp), 2 mitot,mjtot,nghost, 3 delt,dtnew,hx,hy,nvar, 4 xlow,ylow,time,mptr,vtime,naux,alloc(locaux)) if (node(cfluxptr,mptr) .ne. 0) 1 call fluxsv(mptr,alloc(locfm),alloc(locfp), 2 alloc(locgm),alloc(locgp), 3 alloc(node(cfluxptr,mptr)),mitot,mjtot, 4 nvar,listsp(level),delt,hx,hy) if (node(ffluxptr,mptr) .ne. 0) then lenbc = 2*(nx/intrat(level-1)+ny/intrat(level-1)) locsvf = node(ffluxptr,mptr) call fluxad(alloc(locfm),alloc(locfp), 1 alloc(locgm),alloc(locgp), 2 alloc(locsvf),mptr,mitot,mjtot,nvar, 4 lenbc,intrat(level-1),nghost,delt,hx,hy) endif c call reclam(locfp, mitot*mjtot*nvar) call reclam(locfm, mitot*mjtot*nvar) call reclam(locgp, mitot*mjtot*nvar) call reclam(locgm, mitot*mjtot*nvar) c dtlevnew = dmin1(dtlevnew,dtnew) c rnode(timemult,mptr) = rnode(timemult,mptr)+delt mptr = node(levelptr, mptr) if (mptr .ne. 0) go to 5 c return end c c -------------------------------------------------------------- c subroutine bound(time,nvar,ng,valbig,mitot,mjtot,mptr, 1 aux,naux) c implicit double precision (a-h,o-z) include "call.i" dimension valbig(mitot,mjtot,nvar), aux(mitot,mjtot,naux) c c :::::::::::::: BOUND ::::::::::::::::::::::::::::::::::::::::::: c This routine sets the boundary values for a given grid c at level level. c We are setting the values for a strip ng zones wide all c the way around the border, in 4 rectangular strips. c c Outputs from this routine: c The values around the border of the grid are inserted c directly into the enlarged valbig array. c c This routine calls the routine filpatch c which for any block of mesh points on a given level, c intersects that block with all grids on that level and with c the physical boundaries, copies the values into the c appropriate intersecting regions, and interpolates the remaining c cells from coarser grids as required. c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: xleft = rnode(cornxlo, mptr) xright = rnode(cornxhi, mptr) ybot = rnode(cornylo, mptr) ytop = rnode(cornyhi, mptr) ilo = node(ndilo, mptr) ihi = node(ndihi, mptr) jlo = node(ndjlo, mptr) jhi = node(ndjhi, mptr) level = node(nestlevel, mptr) hx = hxposs(level) hy = hyposs(level) c left boundary xl = xleft - ng*hx xr = xleft yb = ybot yt = ytop if ((xl .lt. xlower) .and. xperdom) then call prefilp(level,nvar,valbig,aux,naux,time,mitot,mjtot, 1 1,ng+1, 2 ilo-ng,ilo-1,jlo,jhi) else call filpatch(level,nvar,valbig,aux,naux,time,mitot,mjtot, 1 1,ng+1, 2 ilo-ng,ilo-1,jlo,jhi) endif c right boundary xl = xright xr = xright + ng*hx yb = ybot yt = ytop if ((xr .gt. xupper) .and. xperdom) then call prefilp(level,nvar,valbig,aux,naux,time,mitot,mjtot, 1 mitot-ng+1,ng+1, 2 ihi+1,ihi+ng,jlo,jhi) else call filpatch(level,nvar,valbig,aux,naux,time,mitot,mjtot, 1 mitot-ng+1,ng+1, 2 ihi+1,ihi+ng,jlo,jhi) endif c bottom boundary xl = xleft - ng*hx xr = xright + ng*hx yb = ybot - ng*hy yt = ybot if (((yb .lt. ylower) .and. yperdom) .or. 1 (((xl .lt. xlower) .or. (xr .gt. xupper)) .and. xperdom) ) then call prefilp(level,nvar,valbig,aux,naux,time,mitot,mjtot, 1 1,1, 2 ilo-ng,ihi+ng,jlo-ng,jlo-1) else call filpatch(level,nvar,valbig,aux,naux,time,mitot,mjtot, 1 1,1, 2 ilo-ng,ihi+ng,jlo-ng,jlo-1) endif c top boundary xl = xleft - ng*hx xr = xright + ng*hx yb = ytop yt = ytop + ng*hy if (((yt .gt. yupper) .and. yperdom) .or. 1 (((xl .lt. xlower) .or. (xr .gt. xupper)) .and. xperdom) ) then call prefilp(level,nvar,valbig,aux,naux,time,mitot,mjtot, 1 1,mjtot-ng+1, 2 ilo-ng,ihi+ng,jhi+1,jhi+ng) else call filpatch(level,nvar,valbig,aux,naux,time,mitot,mjtot, 1 1,mjtot-ng+1, 2 ilo-ng,ihi+ng,jhi+1,jhi+ng) endif c c external boundary conditions c if (.not. (xperdom .and. yperdom)) then xl = xleft - ng*hx yb = ybot - ng*hy xr = xright + ng*hx yt = ytop + ng*hy call bc2amr( valbig,aux,mitot,mjtot,nvar,naux, 1 hx,hy,level,time,xl,xr,yb,yt, 3 xlower,ylower,xupper,yupper, 4 xperdom,yperdom) endif c return end c c ------------------------------------------------------------- c subroutine stepgrid(q,fm,fp,gm,gp,mitot,mjtot,mbc,dt,dtnew,dx,dy, & nvar,xlow,ylow,time,mptr,vtime,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 ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: implicit double precision (a-h,o-z) external rpn2,rpt2 include "call.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, vtime, 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 do i = 1, mitot do j = 1, mjtot write(outunit,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 call b4step2(mx,my,mbc,mx,my,nvar,q, & xlowmbc,ylowmbc,dx,dy,time,dt,maux,aux) 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 write(outunit,*) ' Courant # of grid ',mptr, ' is ',cflgrid c cflmax = dmax1(cflmax,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 if (method(5).eq.1) then c # with source term: use Godunov splitting call src2(mx,my,nvar,mbc,mx,my,xlower,ylower,dx,dy, & q,maux,aux,time,dt) endif c c 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 (vtime) then 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 else dtnew = dt endif c # give a warning if Courant number too large... c if (cflgrid .gt. cflv1) then write(*,810) cflgrid write(outunit,810) cflgrid 810 format('*** WARNING *** Courant number =', d12.4, & ' is larger than cflv(1) ') endif c return end c c ---------------------------------------------------------------- c subroutine auxcoarsen(auxdub,midub,mjdub,auxbgc, 1 mi2tot,mj2tot,naux,auxtype) implicit double precision (a-h, o-z) dimension auxdub(midub, mjdub, naux) dimension auxbgc(mi2tot,mj2tot,naux) character*10 auxtype(naux) c :::::::::::::::::::::::: COARSEN :::::::::::::::::::::::::::::::: c coarsen = coarsen the fine grid auxiliary data (with double the usual c number of ghost cells to prepare coarsened data c for error estimation. c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: do 50 iaux = 1, naux if (auxtype(iaux) .eq. "center" .or. . auxtype(iaux) .eq. "capacity") then do 20 j = 1, mj2tot jfine = 2*(j-1) + 1 do 20 i = 1, mi2tot ifine = 2*(i-1) + 1 auxbgc(i,j,iaux) = (auxdub(ifine,jfine,iaux) + & auxdub(ifine+1,jfine,iaux)+ & auxdub(ifine,jfine+1,iaux) + & auxdub(ifine+1,jfine+1,iaux))/4.d0 20 continue elseif (auxtype(iaux) .eq. "xleft") then do 10 j = 1, mj2tot jfine = 2*(j-1) + 1 do 10 i = 1, mi2tot ifine = 2*(i-1) + 1 auxbgc(i,j,iaux) = (auxdub(ifine,jfine,iaux) + & auxdub(ifine,jfine+1,iaux)) /2.d0 10 continue elseif (auxtype(iaux) .eq. "yleft") then do 15 j = 1, mj2tot jfine = 2*(j-1) + 1 do 15 i = 1, mi2tot ifine = 2*(i-1) + 1 auxbgc(i,j,iaux) = (auxdub(ifine,jfine,iaux) + & auxdub(ifine+1,jfine,iaux))/2.d0 15 continue endif 50 continue return end c subroutine fixcapaq(val,aux,mitot,mjtot,valc,auxc,mic,mjc, & nvar,naux,lratio) implicit double precision (a-h,o-z) include "call.i" c c ::::::::::::::::::::::: FIXCAPAQ :::::::::::::::::::::::::::::: c new fine grid solution q was linearly interpolated. but want c to conserve kappa*q, not q. calculate the discrepancy c in kappa*q using this q, and modify q to account for it. c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: dimension val(mitot,mjtot,nvar), valc(mic,mjc,nvar) dimension aux(mitot,mjtot,naux), auxc(mic,mjc,naux) dcapamax = 0.d0 do 10 ic = 2, mic-1 do 10 jc = 2, mjc-1 do 15 ivar = 1, nvar capaqfine = 0.d0 do 20 ico = 1, lratio ifine = (ic-2)*lratio + nghost + ico do 20 jco = 1, lratio jfine = (jc-2)*lratio + nghost + jco capaqfine = capaqfine + aux(ifine,jfine,mcapa)* & val(ifine,jfine,ivar) 20 continue dcapaq = auxc(ic,jc,mcapa)*valc(ic,jc,ivar)- & capaqfine/(lratio*lratio) dcapamax = dmax1(dcapamax,dabs(dcapaq)) do 30 ico = 1, lratio ifine = (ic-2)*lratio + nghost + ico do 30 jco = 1, lratio jfine = (jc-2)*lratio + nghost + jco val(ifine,jfine,ivar) = val(ifine,jfine,ivar) + & dcapaq/aux(ifine,jfine,mcapa) 30 continue 15 continue 10 continue c write(6,*)" max discrepancy ", dcapamax return end c c----------------------------------------------------------------------- c subroutine estdt(val,mitot,mjtot,nvar,dx,dy,dt,nghost,aux,naux) c c :::::::::::::::::::::::: ESTDT :::::::::::::::::::::::::::; c estimate the initial time step for the given values c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::; implicit double precision (a-h, o-z) dimension val(mitot,mjtot,nvar) dimension aux(mitot,mjtot,naux) c c return end c c ---------------------------------------------------------- c function igetsp (nwords) c implicit double precision (a-h,o-z) include "call.i" c c ::::::::::::::::::::::::::: IGETSP :::::::::::::::::::::::::::: c c allocate contiguous space of length nword in main storage array c alloc. user code (or pointer to the owner of this storage) c is mptr. lenf = current length of lfree list. c c ::::::::::::::::::::::::::: IGETSP :::::::::::::::::::::::::::: c c find first fit from free space list c itake = 0 do 20 i = 1, lenf if (lfree(i,2) .lt. nwords) go to 20 itake = i go to 25 20 continue go to 900 c c anything left? c 25 left = lfree(itake,2) - nwords igetsp = lfree(itake,1) iendtake = lfree(itake,1) + nwords if (lendim .lt. iendtake) lendim = iendtake c c the following code which is ignored for now adds the new if (left .le. 0) go to 30 lfree(itake,2) = left lfree(itake,1) = iendtake go to 99 c c item is totally removed. move next items in list up one. c 30 lenf = lenf - 1 do 40 i = itake, lenf lfree(i,1) = lfree(i+1,1) 40 lfree(i,2) = lfree(i+1,2) go to 99 c 900 write(outunit,901) nwords write(*,901) nwords 901 format(' require ',i10,' words - either none left or not big', 1 ' enough space') write(outunit,902) ((lfree(i,j),j=1,2),i=1,lenf) write(*,902) ((lfree(i,j),j=1,2),i=1,lenf) 902 format(' free list: ',//,2x,50(i10,4x,i10,/,2x)) stop c 99 lentot = lentot + nwords if (lenmax .lt. lentot) lenmax = lentot if (sprint) write(outunit,100) nwords, igetsp, lentot, lenmax 100 format(' allocating ',i8,' words in location ',i8, 1 ' lentot, lenmax ', 2i10) return end c c ----------------------------------------------------------- c subroutine reclam (index, nwords) c c ::::::::::::::::::::::::: RECLAM ::::::::::::::::::::::::::: c c return of space. add to free list. c iplace points to next item on free list with larger index than c the item reclaiming, unless said item is greater then c everything on the list. c c ::::::::::::::::::::::::::::::::::;::::::::::::::::::::::::: c implicit double precision (a-h,o-z) include "call.i" c do 20 i = 1, lenf iplace = i if (lfree(i,1) .gt. index) go to 30 20 continue write(outunit,902) write(*,902) 902 format(' no insertion pointer into freelist. error stop') stop c c check previous segment for merging c 30 iprev = iplace - 1 if (lfree(iprev,1)+lfree(iprev,2) .lt. index) go to 40 lfree(iprev,2) = lfree(iprev,2) + nwords go to 50 c c check after segment - no previous merge case c 40 nexti = index + nwords if (lfree(iplace,1).ne. nexti) go to 70 lfree(iplace,1) = index lfree(iplace,2) = lfree(iplace,2) + nwords go to 99 c c check following segment - yes previous merge case c 50 nexti = index + nwords if (lfree(iplace,1) .ne. nexti) go to 99 c c forward merge as well, bump all down 1 c lfree(iprev,2) = lfree(iprev,2)+lfree(iplace,2) ipp1 = iplace + 1 do 60 i = ipp1, lenf lfree(i-1,1) = lfree(i,1) 60 lfree(i-1,2) = lfree(i,2) lenf = lenf - 1 go to 99 c c no merges case - insert and bump future segments up to make room c 70 if (lenf .eq. lfdim) go to 900 do 80 ii = iplace, lenf i = lenf + 1 - ii + iplace lfree(i,1) = lfree(i-1,1) 80 lfree(i,2) = lfree(i-1,2) lenf = lenf + 1 lfree(iplace,1) = index lfree(iplace,2) = nwords go to 99 c 900 write(outunit,901) lfdim write(*,901) lfdim 901 format(' free list full with ',i5,' items') stop c 99 lentot = lentot - nwords if (sprint) write(outunit,100) nwords, index, lentot 100 format(' reclaiming ',i8,' words at loc. ',i8,' lentot ',i10) return end c c -------------------------------------------------- c subroutine birect(mptr1) c implicit double precision (a-h,o-z) include "call.i" c c ::::::::::::: BIRECT ::::::::::::::::::::::::::::::::::::::: c check each grid, starting with mptr1 (either newstl or lstart) c to see that it has no more than max1d points in either dimensions. c needed so that scratch array space in stepgrid not exceeded. c c also check for too small grids - but has never happened. c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c mptr = mptr1 level = node(nestlevel,mptr) hx = hxposs(level) hy = hyposs(level) c 10 continue cxlo = rnode(cornxlo,mptr) cxhi = rnode(cornxhi,mptr) cylo = rnode(cornylo,mptr) cyhi = rnode(cornyhi,mptr) nx = node(ndihi,mptr) - node(ndilo,mptr) + 1 ny = node(ndjhi,mptr) - node(ndjlo,mptr) + 1 minsize = 2*nghost c c check number of rows first - if too many, bisect grid with vertical c line down the middle. make sure new grid corners are anchored c on coarse grid point. make sure if bisecting coarse grid that c new grids have even number of points c if (nx + 2*nghost .gt. max1d) then nxl = nx/2 if (level .gt. 1) then lratio = intrat(level-1) else lratio = 2 endif nxl = (nxl/lratio)*lratio nxr = nx - nxl cxmid = cxlo + nxl*hx mptrnx = nodget(dummy) node(levelptr,mptrnx) = node(levelptr,mptr) node(levelptr,mptr) = mptrnx rnode(cornxhi,mptr) = cxmid node(ndihi,mptrnx) = node(ndihi,mptr) node(ndihi,mptr) = node(ndilo,mptr) + nxl - 1 node(ndilo,mptrnx) = node(ndihi,mptr) + 1 node(ndjhi,mptrnx) = node(ndjhi,mptr) node(ndjlo,mptrnx) = node(ndjlo,mptr) rnode(cornxlo,mptrnx) = cxmid rnode(cornylo,mptrnx) = cylo rnode(cornyhi,mptrnx) = cyhi rnode(cornxhi,mptrnx) = cxhi rnode(timemult,mptrnx) = rnode(timemult,mptr) node(nestlevel,mptrnx) = node(nestlevel,mptr) go to 10 c c check number of columns next - if too many, bisect grid with horizontal c line down the middle c else if (ny + 2*nghost .gt. max1d) then nyl = ny/2 if (level .gt. 1) then lratio = intrat(level-1) else lratio = 2 endif nyl = (nyl/lratio)*lratio nyr = ny - nyl cymid = cylo + nyl*hy mptrnx = nodget(dummy) node(levelptr,mptrnx) = node(levelptr,mptr) node(levelptr,mptr) = mptrnx rnode(cornyhi,mptr) = cymid node(ndjhi,mptrnx) = node(ndjhi,mptr) node(ndjhi,mptr) = node(ndjlo,mptr) + nyl - 1 node(ndjlo,mptrnx) = node(ndjhi,mptr) + 1 node(ndihi,mptrnx) = node(ndihi,mptr) node(ndilo,mptrnx) = node(ndilo,mptr) rnode(cornxlo,mptrnx) = cxlo rnode(cornylo,mptrnx) = cymid rnode(cornyhi,mptrnx) = cyhi rnode(cornxhi,mptrnx) = cxhi node(nestlevel,mptrnx) = node(nestlevel,mptr) rnode(timemult,mptrnx) = rnode(timemult,mptr) go to 10 c c grid ok - check the next c else mptr = node(levelptr,mptr) if (mptr.ne.0) go to 10 endif return end c c --------------------------------------------------------- c subroutine check(nsteps,time,nvar,naux) c c :::::::::::::::::::::: CHECK ::::::::::::::::::::::::::::::::; c check point routine - can only call at end of coarse grid cycle c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::; implicit double precision (a-h,o-z) character chkname*12 include "call.i" c c ### make the file name showing the time step c chkname = 'fort.chkxxxx' nstp = nsteps do 20 ipos = 12, 9, -1 idigit = mod(nstp,10) chkname(ipos:ipos) = char(ichar('0') + idigit) nstp = nstp / 10 20 continue open(unit=chkunit,file=chkname,status='unknown', . form='unformatted') c c ### dump the data c write(chkunit) lenmax,lendim,memsize,(alloc(i),i=1,lendim) write(chkunit) hxposs,hyposs,possk,icheck write(chkunit) lfree,lenf write(chkunit) rnode,node,lstart,newstl,listsp,tol, 1 ibuff,mstart,ndfree,lfine,iorder,mxnest, 2 intrat,iregsz,jregsz,kcheck,nsteps,time,matlabu write(chkunit) evol, rvol, rvoll, lentot c close(chkunit) c return end c c --------------------------------------------------------- c subroutine cleanup(nvar,naux) c c :::::::::::::::::::::: CLEANUP ::::::::::::::::::::::::::::::::; c this is just a check to make sure all storage was accounted for. c routine is called after all the data has been checkpointed. c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::; implicit double precision (a-h,o-z) include "call.i" c c ## clean up storage to double check that everything taken care of c ## done after the checkpoint so pointers sitll work on restart do 120 level = 1, lfine call putsp(1,level,nvar,naux) mptr = lstart(level) 110 nx = node(ndihi,mptr) - node(ndilo,mptr) + 1 ny = node(ndjhi,mptr) - node(ndjlo,mptr) + 1 mitot = nx + 2*nghost mjtot = ny + 2*nghost nwords = mitot*mjtot*nvar call reclam(node(store1, mptr), nwords) if (level .lt. mxnest) . call reclam(node(store2, mptr), nwords) if (naux .gt. 0) . call reclam(node(storeaux, mptr), mitot*mjtot*naux) mptr = node(levelptr, mptr) if (mptr .ne. 0) go to 110 120 continue return end c c ----------------------------------------------------------- c subroutine colate (badpts, len, lcheck, 1 iflags,domflags,isize,jsize,npts) c implicit double precision (a-h,o-z) dimension badpts(2,len) integer iflags (0:isize+1,0:jsize+1) integer domflags(0:isize+1,0:jsize+1) include "call.i" c c c ************************************************************* c c colate = takes the error plane with flagged pts at level lcheck c and puts their (i,j) cell centered c indices into the badpts array. c To insure proper nesting, get rid of flagged point c that don't fit into properly nested domain (in iflags2) c c ************************************************************* c c # if pt. flagged but domain not flagged, turn it off c # note that this results in flags of 1, not 2 of 3. if (dprint) then write(outunit,*)" from colate: iflags" do 48 jj = 1, jsize j = jsize + 1 - jj write(outunit,101)(iflags(i,j),i=1,isize) 48 continue write(outunit,*)" from colate: domflags" do 49 jj = 1, jsize j = jsize + 1 - jj write(outunit,101)(domflags(i,j),i=1,isize) 101 format(80i1) 49 continue endif do 10 j = 1, jsize do 10 i = 1, isize iflags(i,j) = min(iflags(i,j),domflags(i,j)) 10 continue index = 0 c c give points the indices from integer region space. do 20 j = 1, jsize do 20 i = 1, isize if (iflags(i,j) .ne. goodpt) then index = index + 1 badpts(1,index) = float(i)-.5 badpts(2,index) = float(j)-.5 endif 20 continue c 99 npts = index if (gprint) then write(outunit,100) npts, lcheck 100 format( i5,' flagged points colated on level ',i4) endif return end c c ------------------------------------------------------------- c subroutine errest (nvar,naux,numbad,lcheck,iflags,isize,jsize) c implicit double precision (a-h,o-z) include "call.i" dimension iflags (0:isize+1,0:jsize+1) logical vtime data vtime/.false./ c c # no sense computing new time step if just for error estimation, c # so vtime set to false here. c :::::::::::::::::::::::::: ERREST ::::::::::::::::::::::::::::::::::: c for all grids at level lcheck: c estimate the error by taking a large (2h,2k) step based on the c values in the old storage loc., then take one regular (and for c now wasted) step based on the new info. compare using an c error relation for a pth order accurate integration formula. c flag error plane as either bad (needs refinement), or good. c c call the regular integrator on a grid twice as coarse. c initialize such a grid directly, instead of trick dimensioning. c this is to make other l1 type estimates easier to program. c ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c numbad = 0 do 4 i = 1, isize do 4 j = 1, jsize 4 iflags(i,j) = 0 c hx = hxposs(lcheck) hy = hyposs(lcheck) hx2 = 2.d0*hx hy2 = 2.d0*hy mptr = lstart(lcheck) 5 continue nx = node(ndihi,mptr) - node(ndilo,mptr) + 1 ny = node(ndjhi,mptr) - node(ndjlo,mptr) + 1 mitot = nx + 2*nghost mjtot = ny + 2*nghost locnew = node(store1,mptr) locold = node(store2,mptr) locaux = node(storeaux,mptr) mi2tot = nx/2 + 2*nghost mj2tot = ny/2 + 2*nghost time = rnode(timemult,mptr) dt = possk(lcheck) tpre = time - dt c c prepare double the stencil size worth of boundary values, c then coarsen them for the giant step integration. c midub = nx+4*nghost mjdub = ny+4*nghost locdub = igetsp(midub*mjdub*(nvar+naux)) locbgc = igetsp(mi2tot*mj2tot*(nvar+naux)) node(errptr,mptr) = locbgc ng2 = 2*nghost c # transfer soln. into grid with twice the ghost cells call copysol(alloc(locdub),alloc(locold),nvar,mitot,mjtot, 1 nghost,midub,mjdub,ng2) c if (naux .gt. 0) then locaxb = locdub + midub*mjdub*nvar locaxc = locbgc + mi2tot*mj2tot*nvar xl = rnode(cornxlo, mptr) yb = rnode(cornylo, mptr) maxmx = midub - 4*nghost mx = maxmx maxmy = mjdub - 4*nghost my = maxmy call setaux(maxmx,maxmy,2*nghost,mx,my,xl,yb,hx,hy, & naux,alloc(locaxb)) call auxcoarsen(alloc(locaxb),midub,mjdub, 1 alloc(locaxc),mi2tot,mj2tot,naux,auxtype) else locaxb = 1 endif c # fill it - use enlarged (before coarsening) aux arrays call bound(tpre,nvar,ng2,alloc(locdub),midub,mjdub,mptr, 1 alloc(locaxb),naux) c coarsen by 2 in every direction call coarsen(alloc(locdub),midub,mjdub, 1 alloc(locbgc),mi2tot,mj2tot,nvar) call reclam(locdub,midub*mjdub*(nvar+naux)) c c c We now fill bndry values at time t = time, in preparation c for calculating the solution on the grid mptr for error estimation. c locbig = igetsp(mitot*mjtot*nvar) node(tempptr,mptr) = locbig c # straight copy into scratch array so don't mess up latest soln. do 10 i = 1, mitot*mjtot*nvar 10 alloc(locbig+i-1) = alloc(locnew+i-1) call bound(time,nvar,nghost,alloc(locbig),mitot,mjtot,mptr, 1 alloc(locaux),naux) c mptr = node(levelptr,mptr) if (mptr .ne. 0) go to 5 c hx = hxposs(lcheck) hy = hyposs(lcheck) hx2 = 2.d0*hx hy2 = 2.d0*hy dt = possk(lcheck) dt2 = 2. * dt mptr = lstart(lcheck) 25 continue nx = node(ndihi,mptr) - node(ndilo,mptr) + 1 ny = node(ndjhi,mptr) - node(ndjlo,mptr) + 1 mitot = nx+ 2*nghost mjtot = ny+ 2*nghost mi2tot = nx/2 + 2*nghost mj2tot = ny/2 + 2*nghost c c # this scratch storage will be used both for regular and half c # sized grid. different dimensions in stepgrid - do not reuse. locfp = igetsp(mitot*mjtot*nvar) locfm = igetsp(mitot*mjtot*nvar) locgp = igetsp(mitot*mjtot*nvar) locgm = igetsp(mitot*mjtot*nvar) locaux = node(storeaux,mptr) c locbgc = node(errptr,mptr) c locaxc = locbgc+nvar*mi2tot*mj2tot c should we set to 1 if naux=0? c evol = evol + (nx/2)*(ny/2) xlow = rnode(cornxlo,mptr) - nghost*hx2 ylow = rnode(cornylo,mptr) - nghost*hy2 call stepgrid(alloc(locbgc),alloc(locfm),alloc(locfp), 1 alloc(locgm),alloc(locgp), 2 mi2tot,mj2tot,nghost, 3 dt2,dtnew2,hx2,hy2,nvar, 4 xlow,ylow,tpre,mptr,vtime,naux,alloc(locaxc)) c c the one giant step based on old values is done. now take c one regular step based on new values. c evol = evol + nx * ny xlow = rnode(cornxlo,mptr) - nghost*hx ylow = rnode(cornylo,mptr) - nghost*hy locbig = node(tempptr,mptr) loctmp = node(store2,mptr) c c estimate spatial component of error - use old values before c integration to get accurate boundary gradients c locerrsp = igetsp(mitot*mjtot) call errsp(alloc(locbig), alloc(locerrsp), mitot,mjtot, & nvar, mptr,nghost, eprint, outunit) call stepgrid(alloc(locbig),alloc(locfm),alloc(locfp), 1 alloc(locgm),alloc(locgp), 2 mitot,mjtot,nghost, 3 dt,dtnew,hx,hy,nvar, 4 xlow,ylow,time,mptr,vtime,naux,alloc(locaux)) c call reclam(locfp, mitot*mjtot*nvar) call reclam(locfm, mitot*mjtot*nvar) call reclam(locgp, mitot*mjtot*nvar) call reclam(locgm, mitot*mjtot*nvar) c call errf1(alloc(locbig),nvar, 1 alloc(locbgc),mptr,mi2tot,mj2tot, 2 mitot,mjtot,alloc(loctmp),alloc(locerrsp)) call reclam(locbgc,mi2tot*mj2tot*(nvar+naux)) call reclam(locerrsp, mitot*mjtot) call reclam(locbig,mitot*mjtot*nvar) c mptr = node(levelptr, mptr) if (mptr .ne. 0) go to 25 c c copy flagged arrays in individual grids (now stored in loctmp) c into 1 big iflag array (iflags) c mptr = lstart(lcheck) 41 continue nx = node(ndihi,mptr) - node(ndilo,mptr) + 1 ny = node(ndjhi,mptr) - node(ndjlo,mptr) + 1 mitot = nx + 2*nghost mjtot = ny + 2*nghost loctmp = node(store2, mptr) call setflags(iflags,isize,jsize, . alloc(loctmp),nvar,mitot,mjtot,mptr) mptr = node(levelptr,mptr) if (mptr .ne. 0) go to 41 if (eprint) then write(outunit,*)" flagged points before buffering on level", . lcheck do 47 jj = 1, jsize j = jsize + 1 - jj write(outunit,100)(iflags(i,j),i=1,isize) 100 format(80i1) 47 continue endif c c project finer grids to insure level nesting numpro = 0 if (lcheck+2 .le. maxlv) then lrat2 = intrat(lcheck)*intrat(lcheck+1) call projec(lcheck,numpro,iflags,isize,jsize,lrat2) endif if (eprint) then write(outunit,*)" flagged points after projecting to level", . lcheck write(outunit,*) " with ",numpro," additional points projected" do 49 jj = 1, jsize j = jsize + 1 - jj write(outunit,100)(iflags(i,j),i=1,isize) 49 continue endif c c diffuse flagged points in all 4 directions to make buffer zones c note that this code flags with a same value as true flagged c points, not a different number. c c # first get scratch work space (not that other scratch c # arrays above have been reclaimed. c ldom3 = igetsp((isize+2)*(jsize+2)) c do 55 inum = 1, ibuff call shiftset(iflags,alloc(ldom3),+1,0,isize,jsize) call shiftset(iflags,alloc(ldom3),-1,0,isize,jsize) call shiftset(iflags,alloc(ldom3),0,+1,isize,jsize) call shiftset(iflags,alloc(ldom3),0,-1,isize,jsize) 55 continue if (eprint) then write(outunit,*)" flagged points after buffering on level", . lcheck do 48 jj = 1, jsize j = jsize + 1 - jj write(outunit,100)(iflags(i,j),i=1,isize) 48 continue endif c c count up c numbad = 0 do 82 i = 1, isize do 82 j = 1, jsize if (iflags(i,j) .ne. goodpt) numbad = numbad + 1 82 continue write(outunit,116) numbad, lcheck 116 format(i5,' points flagged on level ',i4) call reclam(ldom3,(isize+2)*(jsize+2)) return end c c ---------------------------------------------------------- c subroutine errsp(rect,sperr,mitot,mjtot,nvar,mptr,ng, & eprint, outunit) c c ::::::::::::::::::::: ERRSP :::::::::::::::::::::::::::::::::: c estimate spatial only component of the error c rect = grid values including ghost cells c sperr = user-supplied function indicating where to flag cells c (if sperr(i,j) > tolsp cell will be flagged) c In the default version, the maximum magnitude of the jump c between grid cells over all components is computed. c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c implicit double precision (a-h, o-z) c include "call.i" dimension rect(mitot,mjtot,nvar) dimension sperr(mitot,mjtot) integer outunit logical eprint do 20 j = ng+1, mjtot-ng do 10 i = ng+1, mitot-ng dq = 0.d0 do 5 m = 1,nvar dqi = dabs(rect(i+1,j,m) - rect(i-1,j,m)) dqj = dabs(rect(i,j+1,m) - rect(i,j-1,m)) dq = dmax1(dq,dqi, dqj) 5 continue sperr(i,j) = dq 10 continue 20 continue return end c c ----------------------------------------------------------- c subroutine gfixup(lbase, lfnew, nvar, naux) c implicit double precision (a-h,o-z) include "call.i" c c ::::::::::::::::::::::::: GFIXUP ::::::::::::::::::::::::::::::::; c interpolate initial values for the newly created grids. c the start of each level is located in newstl array. c since only levels greater than lbase were examined, start c looking there. c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::; c c reclaim old storage (position 8) and list space 15 and 16 c before allocating new storage. remember, finest level grids c (if level = mxnest so that error never estimated) don't have c 2 copies of solution values at old and new times. c c call putsp(lbase,lbase,nvar,naux) level = lbase + 1 1 if (level .gt. lfine) go to 4 call putsp(lbase,level,nvar,naux) mptr = lstart(level) 2 if (mptr .eq. 0) go to 3 nx = node(ndihi,mptr) - node(ndilo,mptr) + 1 ny = node(ndjhi,mptr) - node(ndjlo,mptr) + 1 mitot = nx + 2*nghost mjtot = ny + 2*nghost nwords = mitot*mjtot*nvar if (level .lt. mxnest) . call reclam(node(store2, mptr), nwords) node(store2, mptr) = 0 mptr = node(levelptr, mptr) go to 2 3 level = level + 1 go to 1 c 4 lcheck = lbase + 1 5 if (lcheck .gt. maxlv) go to 99 hx = hxposs(lcheck) hy = hyposs(lcheck) c c interpolate level lcheck c mptr = newstl(lcheck) 10 if (mptr .eq. 0) go to 80 nx = node(ndihi,mptr) - node(ndilo,mptr) + 1 ny = node(ndjhi,mptr) - node(ndjlo,mptr) + 1 mitot = nx + 2*nghost mjtot = ny + 2*nghost corn1 = rnode(cornxlo,mptr) corn2 = rnode(cornylo,mptr) loc = igetsp(mitot * mjtot * nvar) node(store1, mptr) = loc if (naux .gt. 0) then locaux = igetsp(mitot * mjtot * naux) maxmx = mitot - 2*nghost mx = maxmx maxmy = mjtot - 2*nghost my = maxmy call setaux(maxmx,maxmy,nghost,mx,my,corn1,corn2,hx,hy, & naux,alloc(locaux)) else locaux = 1 endif node(storeaux, mptr) = locaux time = rnode(timemult, mptr) c c We now fill in the values for grid mptr using filvar. It uses c piecewise linear interpolation to obtain values from the c (lcheck - 1) grid, then overwrites those with whatever (lcheck) c grids are available. We take advantage of the fact that the c (lcheck - 1) grids have already been set, and that the distance c between any point in mptr and a (lcheck - 1) and (lcheck - 2) c interface is at least one (lcheck - 1) cell wide. c c # make a coarsened patch with ghost cells so can use c # grid interpolation routines, but only set "interior". c # extra 2 cells so that can use linear interp. on c # "interior" of coarser patch to fill fine grid. mic = nx/intrat(lcheck-1) + 2 mjc = ny/intrat(lcheck-1) + 2 ivalc = igetsp(mic*mjc*(nvar+naux)) ivalaux = ivalc + nvar*mic*mjc xl = rnode(cornxlo,mptr) xr = rnode(cornxhi,mptr) yb = rnode(cornylo,mptr) yt = rnode(cornyhi,mptr) hx = hxposs(lcheck) hy = hyposs(lcheck) ilo = node(ndilo, mptr) ihi = node(ndihi, mptr) jlo = node(ndjlo, mptr) jhi = node(ndjhi, mptr) call filval(alloc(loc),mitot,mjtot,hx,hy,lcheck,time, 1 alloc(ivalc),alloc(ivalaux),mic,mjc, 2 xl,xr,yb,yt,nvar, 3 mptr,ilo,ihi,jlo,jhi, 4 alloc(locaux),naux) call reclam(ivalc,mic*mjc*(nvar+naux)) mptr = node(levelptr, mptr) go to 10 c c done filling new grids at level. move them into lstart from newstl c (so can use as source grids for filling next level). can also c get rid of loc. 7 storage for old level. c 80 mptr = lstart(lcheck) 85 if (mptr .eq. 0) go to 90 nx = node(ndihi,mptr) - node(ndilo,mptr) + 1 ny = node(ndjhi,mptr) - node(ndjlo,mptr) + 1 mitot = nx + 2*nghost mjtot = ny + 2*nghost call reclam(node(store1,mptr),mitot*mjtot*nvar) if (naux .gt. 0) then call reclam(node(storeaux,mptr),mitot*mjtot*naux) endif mold = mptr mptr = node(levelptr,mptr) call putnod(mold) go to 85 90 lstart(lcheck) = newstl(lcheck) lcheck = lcheck + 1 go to 5 c 99 lfine = lfnew c c initialize 2nd (old time) storage block for new grids not at top level c levend = min(lfine,mxnest-1) do 110 level = lbase+1, levend mptr = lstart(level) 105 if (mptr .eq. 0) go to 110 nx = node(ndihi,mptr) - node(ndilo,mptr) + 1 ny = node(ndjhi,mptr) - node(ndjlo,mptr) + 1 mitot = nx + 2*nghost mjtot = ny + 2*nghost nwords = mitot*mjtot*nvar node(store2,mptr) = igetsp(nwords) mptr = node(levelptr,mptr) go to 105 110 continue c c grid structure now complete again. safe to print, etc. assuming c things initialized to zero in nodget. c return end c c ------------------------------------------------------------------ c subroutine filval(val,mitot,mjtot,hx,hy,lev,time, 1 valc,auxc,mic,mjc, 2 xleft,xright,ybot,ytop,nvar, 3 mptr,ilo,ihi,jlo,jhi,aux,naux) implicit double precision (a-h,o-z) include "call.i" dimension val(mitot,mjtot,nvar), valc(mic,mjc,nvar) dimension aux(mitot,mjtot,nvar), auxc(mic,mjc,nvar) dimension dudx(max1d), dudy(max1d) c c :::::::::::::::::::::::::::::: FILVAL :::::::::::::::::::::::::: c c create and fill coarser (lev-1) patch with one extra coarse cell all c around, plus the ghost cells . will interpolate from this patch to grid mptr c without needing special boundary code. c c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c levc = lev - 1 lratio = intrat(levc) hxcrse = hx*lratio hycrse = hy*lratio xl = xleft - hxcrse xr = xright + hxcrse yb = ybot - hycrse yt = ytop + hycrse c c set integer indices for coarser patch enlarged by 1 cell c (can stick out of domain). proper nesting will insure this one c call is sufficient. iclo = ilo/lratio - 1 jclo = jlo/lratio - 1 ichi = (ihi+1)/lratio - 1 + 1 jchi = (jhi+1)/lratio - 1 + 1 ng = 0 c ::: mcapa is the capacity function index if (mcapa .eq. 0) then if (xperdom .or. yperdom) then call preintcopy(valc,mic,mjc,nvar,iclo,ichi,jclo,jchi,levc) else call intcopy(valc,mic,mjc,nvar,iclo,ichi,jclo,jchi,levc,1,1) endif else if (xperdom .or. yperdom) then call preicall(valc,auxc,mic,mjc,nvar,naux,iclo,ichi,jclo, & jchi,levc) else call icall(valc,auxc,mic,mjc,nvar,naux,iclo,ichi,jclo,jchi, & levc,1,1) endif endif call bc2amr(valc,auxc,mic,mjc,nvar,naux, 1 hxcrse,hycrse,levc,time, 2 xl,xr,yb,yt, 3 xlower,ylower,xupper,yupper, 4 xperdom,yperdom) c c prepare slopes - use min-mod limiters c do 30 j = 2, mjc-1 do 30 ivar = 1, nvar do 10 i = 2, mic-1 slp = valc(i+1,j,ivar) - valc(i,j,ivar) slm = valc(i,j,ivar) - valc(i-1,j,ivar) slopex = dmin1(dabs(slp),dabs(slm))* . dsign(1.0d0,valc(i+1,j,ivar) - valc(i-1,j,ivar)) c # if there's a sign change, set slope to 0. if ( slm*slp .gt. 0.d0) then dudx(i) = slopex else dudx(i) = 0.d0 endif slp = valc(i,j+1,ivar) - valc(i,j,ivar) slm = valc(i,j,ivar) - valc(i,j-1,ivar) slopey = dmin1(dabs(slp),dabs(slm))* . dsign(1.0d0,valc(i,j+1,ivar) - valc(i,j-1,ivar)) if ( slm*slp .gt. 0.d0) then dudy(i) = slopey else dudy(i) = 0.d0 endif 10 continue c c interp. from coarse cells to fine grid c do 20 ico = 1,lratio xoff = (float(ico) - .5)/lratio - .5 do 20 jco = 1,lratio jfine = (j-2)*lratio + nghost + jco yoff = (float(jco) - .5)/lratio - .5 do 20 i = 2, mic-1 ifine = (i-2)*lratio + nghost + ico val(ifine,jfine,ivar) = valc(i,j,ivar) 1 + xoff*dudx(i) + yoff*dudy(i) 20 continue c 30 continue if (mcapa .ne. 0) then call fixcapaq(val,aux,mitot,mjtot,valc,auxc,mic,mjc, & nvar,naux,lratio) endif c c overwrite interpolated values with fine grid values, if available. c call intcopy(val,mitot,mjtot,nvar,ilo-nghost,ihi+nghost, & jlo-nghost,jhi+nghost,lev,1,1) return end c c --------------------------------------------------------------- c subroutine filpatch(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 :::::::::::::::::::::::::::::::::::::::;:::::::::::::::::::::::; implicit double precision (a-h,o-z) include "call.i" logical set, sticksout dimension valbig(mitot,mjtot,nvar), aux(mitot,mjtot,naux) iadflag(i,j) = locuse + i-1+(j-1)*nrowp ivalc(i,j,ivar) = loccrse + (i - 1) + nrowc*(j - 1) & + nrowc*ncolc*(ivar-1) sticksout(iplo,iphi,jplo,jphi) = & (iplo .lt. 0 .or. jplo .lt. 0 .or. & iphi .ge. iregsz(levc) .or. jphi .ge. jregsz(levc)) c c We begin by filling values for grids at level level. If all values can be c filled in this way, we return; hxf = hxposs(level) hyf = hyposs(level) xlp = xlower + ilo*hxf xrp = xlower + (ihi+1)*hxf ybp = ylower + jlo*hyf ytp = ylower + (jhi+1)*hyf nrowp = ihi - ilo + 1 ncolp = jhi - jlo + 1 locuse = igetsp(nrowp*ncolp) c 2/27/02 : added naux to intfil call. call intfil & (valbig,mitot,mjtot,time,locuse,nrowst,ncolst, & ilo,ihi,jlo,jhi,level,nvar,naux) 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 (from bound), after c all 4 boundary pieces filled. call trimbd(alloc(locuse),nrowp,ncolp,set,il,ir,jb,jt) if (set) then call reclam(locuse,nrowp*ncolp) return else if (level .eq. 1) then write(outunit,*)" error in filpatch - level 1 not set" write(outunit,900) nrowst,ncolst write(*,*)" error in filpatch - level 1 not set" write(*,900) nrowst,ncolst 900 format("starting 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 recursive formulation for interpolating. c IS THIS TRUE ? - the fine grid patch remaining unfilled is always c anchored to a coarse cell. 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 lratio = intrat(levc) c iplo = (isl-lratio/2+nghost*lratio)/lratio - nghost c jplo = (jsb-lratio/2+nghost*lratio)/lratio - nghost c iphi = (isr+lratio/2)/lratio c jphi = (jst+lratio/2)/lratio iplo = (isl-lratio +nghost*lratio)/lratio - nghost jplo = (jsb-lratio +nghost*lratio)/lratio - nghost iphi = (isr+lratio )/lratio jphi = (jst+lratio )/lratio 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) loccrse = igetsp(ntot) 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,alloc(locauxc)) endif if ((xperdom .or. yperdom) .and. & sticksout(iplo,iphi,jplo,jphi)) then call prefil2(levc,nvar,alloc(loccrse),alloc(locauxc),naux, 1 time,nrowc,ncolc,1,1, 2 iplo,iphi,jplo,jphi) else call filpatch2(levc,nvar,alloc(loccrse),alloc(locauxc),naux, 1 time,nrowc,ncolc,1,1, 2 iplo,iphi,jplo,jphi) endif c interpolate back up 20 continue do 100 iff = 1,nrowp ic = 2 + (iff - (isl - ilo) - 1)/lratio eta1 = (-0.5d0+dble(mod(iff-1,lratio)))/dble(lratio) do 100 jf = 1,ncolp jc = 2 + (jf -(jsb-jlo)-1)/lratio eta2 = (-0.5d0+dble(mod(jf -1,lratio)))/dble(lratio) flag = alloc(iadflag(iff,jf)) if (flag .eq. 0.0) then c xif = xlp + (.5 + float(iff-1))*hxf c yjf = ybp + (.5 + float(jf -1))*hyf c ic=idint((xif-xlc+.5*hxc)/hxc) c jc=idint((yjf-ybc+.5*hyc)/hyc) c xc = xlc + (float(ic) - .5)*hxc c yc = ybc + (float(jc) - .5)*hyc c eta1 = (xif - xc)/hxc c eta2 = (yjf - yc)/hyc do 101 ivar = 1,nvar c valc00 = alloc(ivalc(ic,jc,ivar)) c valc10 = alloc(ivalc(ic+1,jc,ivar)) c valc01 = alloc(ivalc(ic,jc+1,ivar)) c valc11 = alloc(ivalc(ic+1,jc+1,ivar)) valp10 = alloc(ivalc(ic+1,jc,ivar)) valm10 = alloc(ivalc(ic-1,jc,ivar)) valc = alloc(ivalc(ic ,jc,ivar)) valp01 = alloc(ivalc(ic ,jc+1,ivar)) valm01 = alloc(ivalc(ic ,jc-1,ivar)) dupc = valp10 - valc dumc = valc - valm10 ducc = valp10 - valm10 du = dmin1(dabs(dupc),dabs(dumc)) du = dmin1(2.d0*du,.5d0*dabs(ducc)) fu = dmax1(0.d0,dsign(1.d0,dupc*dumc)) dvpc = valp01 - valc dvmc = valc - valm01 dvcc = valp01 - valm01 dv = dmin1(dabs(dvpc),dabs(dvmc)) dv = dmin1(2.d0*dv,.5d0*dabs(dvcc)) fv = dmax1(0.d0,dsign(1.d0,dvpc*dvmc)) valint = valc + eta1*du*dsign(1.d0,ducc)*fu . + eta2*dv*dsign(1.d0,dvcc)*fv c valint = (1. - eta2)* c & ((1. - eta1)*valc00 + eta1*valc10) c & + eta2*((1. - eta1)*valc01 + eta1*valc11) valbig(iff+nrowst-1,jf+ncolst-1,ivar) = valint 101 continue endif 100 continue call reclam(loccrse,ntot) call reclam(locuse,nrowp*ncolp) return end c c --------------------------------------------------------------- c subroutine filpatch2(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 :::::::::::::::::::::::::::::::::::::::;:::::::::::::::::::::::; implicit double precision (a-h,o-z) include "call.i" logical set, sticksout dimension valbig(mitot,mjtot,nvar), aux(mitot,mjtot,naux) iadflag(i,j) = locuse + i-1+(j-1)*nrowp ivalc(i,j,ivar) = loccrse + (i - 1) + nrowc*(j - 1) & + nrowc*ncolc*(ivar-1) sticksout(iplo,iphi,jplo,jphi) = & (iplo .lt. 0 .or. jplo .lt. 0 .or. & iphi .ge. iregsz(levc) .or. jphi .ge. jregsz(levc)) 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 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 hmargin = dmin1(hxf,hyf)/10.d0 call intfil & (valbig,mitot,mjtot,time,locuse,nrowst,ncolst, & ilo,ihi,jlo,jhi,level,nvar,naux) c c only call bc2amr for coarse patches (where bc2amr sets the whole c patch. otherwise, orig. caller of filpatch2 will take care of this 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 c into account the points filled by c the boundary conditions. bc2amr will be called later, after all 4 c boundary pieces filled. call trimbd(alloc(locuse),nrowp,ncolp,set,il,ir,jb,jt) if (set) then 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) call reclam(locuse,nrowp*ncolp) return else if (level .eq. 1) then write(outunit,*)" error in filpatch2 - level 1 not set" write(outunit,900) nrowst,ncolst write(*,*)" error in filpatch2 - level 1 not set" 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 lratio = intrat(levc) iplo = (isl-lratio+nghost*lratio)/lratio - nghost jplo = (jsb-lratio+nghost*lratio)/lratio - nghost iphi = (isr+lratio)/lratio jphi = (jst+lratio)/lratio 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) loccrse = igetsp(ntot) 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,alloc(locauxc)) endif if ((xperdom .or. yperdom) .and. & sticksout(iplo,iphi,jplo,jphi)) then call prefil3(levc,nvar,alloc(loccrse),alloc(locauxc),naux, 1 time,nrowc,ncolc,1,1, 2 iplo,iphi,jplo,jphi) else call filpatch3(levc,nvar,alloc(loccrse),alloc(locauxc),naux, 1 time,nrowc,ncolc,1,1, 2 iplo,iphi,jplo,jphi) endif c interpolate back up 20 continue do 100 iff = 1,nrowp ic = 2 + (iff-(isl-ilo)-1)/lratio eta1 = (-0.5d0+dble(mod(iff-1,lratio)))/dble(lratio) do 100 jf = 1,ncolp jc = 2 + (jf -(jsb-jlo)-1)/lratio eta2 = (-0.5d0+dble(mod(jf -1,lratio)))/dble(lratio) flag = alloc(iadflag(iff,jf)) if (flag .eq. 0.0) then c xif = xlp + (.5 + float(iff-1))*hxf c yjf = ybp + (.5 + float(jf -1))*hyf c ic=idint((xif-xlc+.5*hxc)/hxc) c jc=idint((yjf-ybc+.5*hyc)/hyc) c xc = xlc + (float(ic) - .5)*hxc c yc = ybc + (float(jc) - .5)*hyc c eta1 = (xif - xc)/hxc c eta2 = (yjf - yc)/hyc do 101 ivar = 1,nvar valp10 = alloc(ivalc(ic+1,jc,ivar)) valm10 = alloc(ivalc(ic-1,jc,ivar)) valc = alloc(ivalc(ic ,jc,ivar)) valp01 = alloc(ivalc(ic ,jc+1,ivar)) valm01 = alloc(ivalc(ic ,jc-1,ivar)) dupc = valp10 - valc dumc = valc - valm10 ducc = valp10 - valm10 du = dmin1(dabs(dupc),dabs(dumc)) du = dmin1(2.d0*du,.5d0*dabs(ducc)) fu = dmax1(0.d0,dsign(1.d0,dupc*dumc)) dvpc = valp01 - valc dvmc = valc - valm01 dvcc = valp01 - valm01 dv = dmin1(dabs(dvpc),dabs(dvmc)) dv = dmin1(2.d0*dv,.5d0*dabs(dvcc)) fv = dmax1(0.d0,dsign(1.d0,dvpc*dvmc)) valint = valc + eta1*du*dsign(1.d0,ducc)*fu . + eta2*dv*dsign(1.d0,dvcc)*fv c valc00 = alloc(ivalc(ic,jc,ivar)) c valc10 = alloc(ivalc(ic+1,jc,ivar)) c valc01 = alloc(ivalc(ic,jc+1,ivar)) c valc11 = alloc(ivalc(ic+1,jc+1,ivar)) c valint = (1. - eta2)* c & ((1. - eta1)*valc00 + eta1*valc10) c & + eta2*((1. - eta1)*valc01 + eta1*valc11) valbig(iff+nrowst-1,jf+ncolst-1,ivar) = valint 101 continue endif 100 continue 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) call reclam(loccrse,ntot) call reclam(locuse,nrowp*ncolp) return end c c --------------------------------------------------------------- c subroutine filpatch3(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 :::::::::::::::::::::::::::::::::::::::;:::::::::::::::::::::::; implicit double precision (a-h,o-z) include "call.i" logical set dimension valbig(mitot,mjtot,nvar), aux(mitot,mjtot,naux) iadflag(i,j) = locuse + i-1+(j-1)*nrowp 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 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,locuse,nrowst,ncolst, & ilo,ihi,jlo,jhi,level,nvar,naux) c c only call bc2amr for coarse patches (where bc2amr sets the whole c patch. otherwise, orig. caller of filpatch3 will take care of this c if (nrowp .eq. mitot .and. ncolp .eq. mjtot) then call bc2amr(valbig,aux,mitot,mjtot,nvar,naux, 1 hxf,hyf,level,tim