|
movetopo_geo.f.html |
 |
|
Source file: movetopo_geo.f
|
|
Directory: /Users/rjl/git/rjleveque/clawpack-4.6.3/geoclaw/2d/lib
|
|
Converted: Mon Jan 21 2013 at 20:15:57
using clawcode2html
|
|
This documentation file will
not reflect any later changes in the source file.
|
c ============================================
subroutine movetopo(maxmx,maxmy,mbc,mx,my,
& xlow,ylow,dx,dy,t,dt,maux,aux,dtopo,
& xlowdtopo,ylowdtopo,xhidtopo,yhidtopo,t0dtopo,tfdtopo,
& dxdtopo,dydtopo,dtdtopo,mxdtopo,mydtopo,mtdtopo,mdtopo,
& minleveldtopo,maxleveldtopo,topoaltered)
c ============================================
c
c This routine changes the topography dynamically.
c It also changes the actual topo array after the motion has ended.
c Topography is stored in aux(i,j,1).
c
c
use geoclaw_module
use topo_module
implicit double precision (a-h,o-z)
dimension aux(1-mbc:maxmx+mbc,1-mbc:maxmy+mbc, maux)
dimension dtopo(mxdtopo,mydtopo,mtdtopo)
logical topoaltered
include "call.i"
dimension auxorig(-1:mx+mbc,-1:my+mbc, maux)
c dimension auxorig(-1:max1d+mbc,-1:max1d+mbc, maxaux)
t0=t !# start of coming timestep
tf=t+dt !# end of coming timestep
t2=t0+0.5d0*dt
tt=t2 !# this is the time in the fault file that will be used for this timestep
c write(26,*) 'MOVETOPO: tt, dt, t0dtopo: ',tt,dt,t0dtopo
if (topoaltered.and.t0.le.tfdtopo) then
write(*,*) 'MOVETOPO: WARNING!'
write(*,*) 'FAULT MOTION HAS COMPLETED PREMATURELY'
endif
if (topoaltered) return
if (tf.lt.t0dtopo) return
if (t0.gt.tfdtopo) go to 28
c # change topography
write(*,*) 'MOVETOPO: setting dtopo at time = ',t
c write(26,*) 'MOVETOPO: setting dtopo at time = ',t
if (tt.ge.tfdtopo) then
kkm = mtdtopo
kkp = mtdtopo
tau = 1.d0
go to 14
endif
if (tt.le.t0dtopo) then
kkm = 1
kkp = 1
tau= 1.d0
go to 14
endif
dkk = (tt-t0dtopo)/dtdtopo
kkm = dint(dkk)+1
kkp = dint(dkk)+2
kkm = max(kkm,1)
kkm = min(kkm,mtdtopo)
kkp = max(kkp,1)
kkp = min(kkp,mtdtopo)
tdtopom = t0dtopo+dtdtopo*(kkm-1)
tdtopop = tdtopom+dtdtopo
tau = 1.d0-max(0.d0,((tt-tdtopom)/dtdtopo))
tau = max(tau,0.d0)
14 continue
c # recreate original topography:
call setaux(mx,my,mbc,mx,my,xlow,ylow,dx,dy,
& maux,auxorig)
c call setaux(max1d,max1d,mbc,mx,my,xlow,ylow,dx,dy,
c & maxaux,auxorig)
c=======loop through the computational grid row by row====================
do j=1-mbc,my+mbc
ycell = ylow + (j-0.5d0)*dy
yjm = ylow + (j-1.d0)*dy
yjp = ylow + dble(j)*dy
c if (yjm.ge.ylowdtopo.and.yjp.le.yhidtopo) then
if (yjp.gt.ylowdtopo.and.yjm.lt.yhidtopo) then
yjpc=min(yjp,yhidtopo)
yjmc=max(yjm,ylowdtopo)
dyc=yjpc-yjmc
ycellc=0.5d0*(yjpc+yjmc)
c=========sweep through the columns of the computational row================
do i=1-mbc,mx+mbc
xcell = xlow + (i- 0.5d0)*dx
xim = xlow + (i- 1.0d0)*dx
xip = xlow + dble(i)*dx
c if (xim.ge.xlowdtopo.and.xip.le.xhidtopo) then
if (xip.gt.xlowdtopo.and.xim.lt.xhidtopo) then
xipc=min(xip,xhidtopo)
ximc=max(xim,xlowdtopo)
dxc=xipc-ximc
xcellc=0.5d0*(xipc+ximc)
c==========alter aux(i,j,1) if it is in the dtopo region=====================
aux(i,j,1) = auxorig(i,j,1)
c==================dynamically alter the topography to the interpolated value=
dauxijm=0.d0
dauxijm = topointegral(ximc,xcellc,xipc,yjmc,
& ycellc,yjpc,xlowdtopo,ylowdtopo,dxdtopo,dydtopo,
& mxdtopo,mydtopo,dtopo(1,1,kkm),1)
dauxijm = dauxijm/(dxc*dyc*aux(i,j,2))
dauxijp=0.d0
dauxijp = topointegral(ximc,xcellc,xipc,yjmc,
& ycellc,yjpc,xlowdtopo,ylowdtopo,dxdtopo,dydtopo,
& mxdtopo,mydtopo,dtopo(1,1,kkp),1)
dauxijp = dauxijp/(dxc*dyc*aux(i,j,2))
aux(i,j,1)=aux(i,j,1)+tau*dauxijm+(1.d0-tau)*dauxijp
endif
enddo
endif
enddo
c============= end of altering topography dynamically========================
c # if the aux cell is outside the physical domain,
c # boundary conditions must be reset
if (xlowdtopo.lt.xlower.and.
& xlow+(0.5d0-mbc)*dx.lt.xlower) then
do j=1-mbc,my+mbc
ycell = ylow +(j-0.5d0)*dy
yjm = ylow + (j-1.0d0)*dy
yjp = ylow + dble(j)*dy
xcell = xlow + (1- 0.5d0)*dx
xim = xlow
xip = xlow + dx
if ((xim.ge.xlowdtopo.and.xip.le.xhidtopo).and.
& (yjm.ge.ylowdtopo.and.yjp.le.yhidtopo)) then
do i=1-mbc,0
aux(i,j,1) = aux(1,j,1)
enddo
endif
enddo
endif
if (xhidtopo.gt.xupper.and.
& xlow+(mx+mbc-0.5d0)*dx.gt.xupper) then
do j=1-mbc,my+mbc
ycell = ylow+(j-0.5d0)*dy
yjm = ylow+(j-1.0d0)*dy
yjp = ylow + dble(j)*dy
xcell = xlow + (mx-0.5d0)*dx
xim = xlow + (mx-1.d0)*dx
xip = xlow + dble(mx)*dx
if ((xim.ge.xlowdtopo.and.xip.le.xhidtopo).and.
& (yjm.ge.ylowdtopo.and.yjp.le.yhidtopo)) then
do i=mx+1,mx+mbc
aux(i,j,1) = aux(mx,j,1)
enddo
endif
enddo
endif
if (ylowdtopo.lt.ylower.and.
& ylow+(.5d0-mbc)*dy.lt.ylower) then
do i=1-mbc,mx+mbc
ycell = ylow+(1-0.5d0)*dy
yjm = ylow
yjp = ylow + dy
xcell = xlow + (i- 0.5d0)*dx
xim = xlow + (i- 1.d0)*dx
xip = xlow + dble(i)*dx
if ((xim.ge.xlowdtopo.and.xip.le.xhidtopo).and.
& (yjm.ge.ylowdtopo.and.yjp.le.yhidtopo)) then
do j=1-mbc,0
aux(i,j,1) = aux(i,1,1)
enddo
endif
enddo
endif
if (yhidtopo.gt.yupper.and.
& ylow+(my+mbc-0.5d0)*dy.gt.yupper) then
do i=1-mbc,mx+mbc
ycell = ylow+(my-0.5d0)*dy
yjm = ylow+(my-1.d0)*dy
yjp = ylow+ dble(my)*dy
xcell = xlow + (i- 0.5d0)*dx
xim = xlow + (i- 1.d0)*dx
xip = xlow + dble(i)*dx
if ((xim.ge.xlowdtopo.and.xip.le.xhidtopo).and.
& (yjm.ge.ylowdtopo.and.yjp.le.yhidtopo)) then
do j=my+1,my+mbc
aux(i,j,1) = aux(i,my,1)
enddo
endif
enddo
endif
c===========================B.C.'s reset========================================
c============= statically change topowork array for next call to setaux====
28 continue
if (tt-dt.gt.tfdtopo) then
c if (.false.) then
topoaltered=.true. !# so this procedure only done once
write(*,*) 'MOVETOPO: Altering topography arrays at t=',t
write(*,*) 'MOVETOPO: Topography Motion Complete'
do m=1,mtopofiles
if (xlowtopo(m).lt.xhidtopo.and.xhitopo(m).gt.xlowdtopo
& .and.ylowtopo(m).lt.yhidtopo.and.
& yhitopo(m).gt.ylowdtopo) then
do ib=1,mxtopo(m)
xcell=xlowtopo(m) + (ib-.5d0)*dxtopo(m)
xim=xlowtopo(m) + (ib-1.d0)*dxtopo(m)
xip=xlowtopo(m) + ib*dxtopo(m)
if (xim.lt.xhidtopo.and.xip.gt.xlowdtopo) then
ximc=max(xim,xlowdtopo)
xipc=min(xip,xhidtopo)
dxc=xipc-ximc
xcellc=0.5d0*(ximc+xipc)
do jb=1,mytopo(m)
ycell=ylowtopo(m) + (jb-.5d0)*dytopo(m)
yjm=ylowtopo(m) + (jb-1.d0)*dytopo(m)
yjp = ylowtopo(m) + jb*dytopo(m)
if (yjm.lt.yhidtopo.and.yjp.gt.ylowdtopo) then
yjmc=max(yjm,ylowdtopo)
yjpc=min(yjp,yhidtopo)
dyc = yjpc-yjmc
ycellc= 0.5d0*(yjmc+yjpc)
ztopoij = 0.d0
ztopoij = topointegral(ximc,xcellc,xipc,yjmc,ycellc,
& yjpc,xlowdtopo,ylowdtopo,dxdtopo,dydtopo,
& mxdtopo,mydtopo,dtopo(1,1,mtdtopo),1)
c # topointegral actually integrates,
c #so need to divide by area
c #physical area of cell depends on coordinate system
ztopoij=ztopoij/(dxc*dyc)
if (icoordsys.eq.2) then
deg2rad = pi/180.d0
capac_area = deg2rad*Rearth**2*
& (sin(yjp*deg2rad)-sin(yjm*deg2rad))/dyc
ztopoij=ztopoij/capac_area
endif
jbr = mytopo(m)-jb+1
ij = i0topo(m) + (jbr-1)*mxtopo(m) +ib -1
topowork(ij) = topowork(ij)+ztopoij
endif
enddo
endif
enddo
endif
enddo
endif
c=================================================================================
return
end