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