c $Header: /home/proj/vtfng/vtf/amroc/rim/src/slope.f,v 1.2 2004/12/01 20:51:35 ralf Exp $
subroutine slope(w,dwx,dwy,dx,dy,nx,ny,nbc)
implicit double precision (a-h,o-z)
c
c IMPORTANT: the dimentions of w MUST match the dimentions
c off all arrays in 'rim2d.i'
c
parameter (MBC=2, NXM=512, NYM=512)
dimension w(1-MBC:NXM+MBC,1-MBC:NYM+MBC),
& dwx(1-MBC:NXM+MBC,1-MBC:NYM+MBC),
& dwy(1-MBC:NXM+MBC,1-MBC:NYM+MBC)
c
vcon=1d-14
c
do j=0,ny+1
do i=0,nx+1
c----------------------------------------right-left -------------
arx = (w(i+1,j)-w(i,j))/dx
alx = (w(i,j)-w(i-1,j))/dx
corxlr = (arx-alx)*(arx-alx)/(alx*alx+arx*arx+vcon)
dwlr = 0.5d0*(alx+arx)*(1d0-corxlr)
if(alx.lt.0d0.and.arx.gt.0d0) dwlr=0d0
c
c----------------------------------------top-bottom -------------
aty = (w(i,j+1)-w(i,j))/dy
aby = (w(i,j)-w(i,j-1))/dy
corybt = (aty-aby)*(aty-aby)/(aty*aty+aby*aby+vcon)
dwbt = 0.5d0*(aty+aby)*(1d0-corybt)
if(aby.lt.0d0.and.aty.gt.0d0) dwbt=0d0
c
c---------------integration along the cell boundary---------------
c
wl = w(i,j) - 0.5d0*dx*dwlr
wr = w(i,j) + 0.5d0*dx*dwlr
wb = w(i,j) - 0.5d0*dy*dwbt
wt = w(i,j) + 0.5d0*dy*dwbt
c
dwx(i,j) =(wr-wl)/dx
dwy(i,j) =(wt-wb)/dy
enddo
enddo
c
return
end