fgridinterp_geo.f.html | ![]() |
Source file: fgridinterp_geo.f | |
Directory: /Users/rjl/git/rjleveque/clawpack-4.6.3/geoclaw/2d/lib | |
Converted: Mon Jan 21 2013 at 20:15:56 using clawcode2html | |
This documentation file will not reflect any later changes in the source file. |
c====================================================================== subroutine fgridinterp(fgrid,xlowfg,ylowfg, & xhifg,yhifg,dxfg,dyfg,mxfg,myfg,t,mvarsfg,q,meqn, & mxc,myc,mbc,dxc,dyc,nvar,xlowc,ylowc,maux,aux, & ioutarrivaltimes,ioutsurfacemax,maxcheck) c====================================================================== use geoclaw_module implicit double precision (a-h,o-z) dimension q(1-mbc:mxc+mbc,1-mbc:myc+mbc, meqn) dimension aux(1-mbc:mxc+mbc,1-mbc:myc+mbc, maux) dimension fgrid(1:mxfg,1:myfg,mvarsfg) c=====================FGRIDINTERP======================================= c # This routine interpolates q and aux on a computational grid c # to a fgrid not necessarily aligned with the computational grid c # using bilinear interpolation defined on computation grid c======================================================================= xhic=xlowc + dxc*mxc yhic=ylowc + dyc*myc tol=drytolerance arrivaltol=1.d-2 indb=meqn+1 indeta=meqn+2 if (maxcheck.gt.0) then indarrive=0 indetamin=0 indetamax=0 indetanow=0 if (ioutarrivaltimes.gt.0) then indarrive = 1 endif if (ioutsurfacemax.gt.0) then indetanow = indarrive+1 indetamin = indarrive+2 indetamax = indarrive+3 endif endif do ifg=1,mxfg xfg=xlowfg + (ifg-1)*dxfg do jfg=1,myfg yfg=ylowfg + (jfg-1)*dyfg if (.not.((xfg.lt.xlowc.or.xfg.gt.xhic).or. & (yfg.lt.ylowc.or.yfg.gt.yhic))) then c # find where xfg,yfg is in the computational grid ic1 = int((xfg-(xlowc+0.5d0*dxc))/(dxc))+1 jc1 = int((yfg-(ylowc+0.5d0*dyc))/(dyc))+1 if (ic1.eq.mxc) ic1=mxc-1 if (jc1.eq.myc) jc1=myc-1 ic2= ic1+1 jc2= jc1+1 xc1= xlowc + dxc*(ic1-0.5d0) yc1= ylowc + dyc*(jc1-0.5d0) xc2= xlowc + dxc*(ic2-0.5d0) yc2= ylowc + dyc*(jc2-0.5d0) c # interpolate bilinear used to interpolate to xfg,yfg c # define constant parts of bilinear xterm=(xfg-xc1)/dxc yterm=(yfg-yc1)/dyc xyterm= xterm*yterm if (maxcheck.eq.0) then do m=1,meqn z11=q(ic1,jc1,m) z21=q(ic2,jc1,m) z12=q(ic1,jc2,m) z22=q(ic2,jc2,m) a=z21-z11 b=z12-z11 d=z11 c=z22-(a+b+d) fgrid(ifg,jfg,m) = a*xterm + b*yterm + c*xyterm + d enddo z11=aux(ic1,jc1,1) z21=aux(ic2,jc1,1) z12=aux(ic1,jc2,1) z22=aux(ic2,jc2,1) a=z21-z11 b=z12-z11 d=z11 c=z22-(a+b+d) fgrid(ifg,jfg,indb) = a*xterm + b*yterm + c*xyterm + d endif c # This next output variable is the surface using bilinear interpolation, c # using a surface that only uses the wet eta points near the shoreline z11=aux(ic1,jc1,1)+q(ic1,jc1,1) z21=aux(ic2,jc1,1)+q(ic2,jc1,1) z12=aux(ic1,jc2,1)+q(ic1,jc2,1) z22=aux(ic2,jc2,1)+q(ic2,jc2,1) h11=q(ic1,jc1,1) h21=q(ic2,jc1,1) h12=q(ic1,jc2,1) h22=q(ic2,jc2,1) depthindicator= min(h11,h12,h21,h22) totaldepth= h11+h22+h21+h12 if (depthindicator.lt.tol.and.totaldepth.gt.4.d0*tol) then !near shoreline if (h11.lt.tol) then z11w= (h12*z12 + h21*z21 + h22*z22)/totaldepth z11=z11w endif if (h12.lt.tol) then z12w= (h11*z11 + h21*z21 + h22*z22)/totaldepth z12=z12w endif if (h21.lt.tol) then z21w= (h11*z11 + h12*z12 + h22*z22)/totaldepth z21=z21w endif if (h22.lt.tol) then z22w= (h12*z12 + h21*z21 + h11*z11)/totaldepth z22=z22w endif endif if (totaldepth.le.4.d0*tol) then z22=d_nan() endif a=z21-z11 b=z12-z11 d=z11 c=z22-(a+b+d) c #If eta max/min are saved on this grid initialized if necessary if (ioutsurfacemax.gt.0.and.maxcheck.eq.2) then if (.not.(fgrid(ifg,jfg,indetamin).eq. & fgrid(ifg,jfg,indetamin))) fgrid(ifg,jfg,indetamin)=0.d0 if (.not.(fgrid(ifg,jfg,indetamax).eq. & fgrid(ifg,jfg,indetamax))) fgrid(ifg,jfg,indetamax)=0.d0 endif c # check wich task to perform if (maxcheck.eq.0) then fgrid(ifg,jfg,indeta) = a*xterm + b*yterm + c*xyterm + d fgrid(ifg,jfg,mvarsfg) = t elseif (maxcheck.eq.1.and.ioutsurfacemax.gt.0) then fgrid(ifg,jfg,indetanow) = a*xterm + b*yterm + c*xyterm + d elseif (maxcheck.eq.2.and.ioutsurfacemax.gt.0) then fgrid(ifg,jfg,indetamin) = & min(fgrid(ifg,jfg,indetamin),fgrid(ifg,jfg,indetanow)) fgrid(ifg,jfg,indetamax) = & max(fgrid(ifg,jfg,indetamax),fgrid(ifg,jfg,indetanow)) endif c #If arrival times are saved on this grid if (maxcheck.eq.1.and.ioutarrivaltimes.gt.0) then check=fgrid(ifg,jfg,indarrive) if (.not.(check.eq.check)) then !# check=NaN: Waves haven't arrived previously if (dabs(fgrid(ifg,jfg,indeta)).gt.arrivaltol) then fgrid(ifg,jfg,indarrive)= t endif endif endif endif enddo enddo return end