copyright (C) 2001 MSC-RPN COMM %%%MC2%%% ***s/r nu_del2_n *subroutine nu_del2_n (rfd,sfd,lminx,lmaxx,lminy,lmaxy,lnk,nu, 6 $ is,js,m,n) implicit none * integer lminx,lmaxx,lminy,lmaxy,lnk,is,js,m,n real rfd (lminx:lmaxx,lminy:lmaxy,lnk), $ sfd (lminx:lmaxx,lminy:lmaxy,lnk) real*8 nu * * *AUTHORs C. Girard & M. Desgagne * *OBJECT * * ******************************************************************* * * * * * OPERATOR nu_DEL2_n * * * * * * n successive calls to this subroutine produces * * * * * * ---------------- * * * | n | * * * the equivalent of a diffusion: | -(-nu_DEL2) | * * * | | * * * ---------------- * * * * * * * * * * * * Each call applies a 9-point filter * * * * * * [ Shuman, M.W.R. #57, p.357-361, eq #5. ] * * * * * * to the difference: rfd - sfd. * * * * * * * * ******************************************************************* * *EXTERNALS #include "lcldim.cdk"
* * ** integer i,j,k,id,jd,iff,jf real wk(minx:maxx,miny:maxy) real*8 c1,c2,c3,one,two,four parameter(one=1.d0,two=2.d0,four=4.d0) *---------------------------------------------------------------------- * id = 1+is*west jd = 1+js*south iff= ldni-east jf = ldnj-north * c1 = nu*(one-two*nu) c2 = nu**2 c3 = nu*four*(nu-one) * !$omp do do k=1,lnk if (m.eq.1) then do j=jd-1,jf+1 do i=id-1,iff+1 sfd(i,j,k) = rfd(i,j,k) end do end do else if (m.eq.2) then do j=jd-1,jf+1 do i=id-1,iff+1 sfd(i,j,k) = rfd(i,j,k) - sfd(i,j,k) end do end do else do j=jd-1+south,jf+1-north do i=id-1+west,iff+1-east sfd(i,j,k) = rfd(i,j,k) - sfd(i,j,k) end do end do endif if (m.eq.n) then do j=jd,jf do i=id,iff rfd(i,j,k)= rfd(i,j,k) + $ c1*(sfd(i ,j+1,k)+sfd(i+1,j ,k) + $ sfd(i ,j-1,k)+sfd(i-1,j ,k))+ $ c2*(sfd(i-1,j+1,k)+sfd(i+1,j+1,k) + $ sfd(i-1,j-1,k)+sfd(i+1,j-1,k))+ $ c3* sfd(i ,j ,k) end do end do else do j=jd,jf do i=id,iff wk(i,j) = c1*(sfd(i ,j+1,k)+sfd(i+1,j ,k) + $ sfd(i ,j-1,k)+sfd(i-1,j ,k))+ $ c2*(sfd(i-1,j+1,k)+sfd(i+1,j+1,k) + $ sfd(i-1,j-1,k)+sfd(i+1,j-1,k))+ $ c3* sfd(i ,j ,k) end do end do do j=jd,jf do i=id,iff sfd(i,j,k)= rfd(i,j,k) + wk(i,j) end do end do endif end do !$omp enddo * *---------------------------------------------------------------------- return end