
      subroutine setmask(mask,q,p,nd)


      integer mask(-(nd+1) : (nd+1), -(nd+1) : (nd+1))
      double precision q(-(nd+1):(nd+1))
      double precision p(-(nd+1):(nd+1))
      integer nd

      integer i,j

c     Set mask=0 for interor points

      forall (i=-nd:nd, j=-nd:nd) mask(i,j) = 0

c     Set mask=11 for x edge and interpolation along constant x'
     
      forall (i=-(nd+1):-(nd+1),j=-(nd+1)/2:(nd+1)/2) mask(i,j)=11
      forall (i=(nd+1):(nd+1), j=-(nd+1)/2:(nd+1)/2) mask(i,j)=11

c     Set mask=12 for x edge  and interpolation along constant y'

      forall (i=-(nd+1):-(nd+1), j=((nd+1)/2)+1:nd) mask(i,j)=12
      forall (i=-(nd+1):-(nd+1), j=-nd:-(((nd+1)/2)+1)) mask(i,j)=12  
      forall (i=(nd+1):(nd+1), j=((nd+1)/2)+1:nd) mask(i,j)=12
      forall (i=(nd+1):(nd+1), j=-nd:-(((nd+1)/2)+1)) mask(i,j)=12   
    
c     Set mask=21 for y edge and interpolation along constant x'

      forall (j=-(nd+1):-(nd+1), i=((nd+1)/2)+1:nd) mask(i,j)=21
      forall (j=-(nd+1):-(nd+1), i=-nd:-(((nd+1)/2)+1)) mask(i,j)=21  
      forall (j=(nd+1):(nd+1), i=((nd+1)/2)+1:nd) mask(i,j)=21
      forall (j=(nd+1):(nd+1), i=-nd:-(((nd+1)/2)+1)) mask(i,j)=21       


c     Set mask=22 for y edge and interpolation along constant y'
      
      forall (j=-(nd+1):-(nd+1),i=-(nd+1)/2:(nd+1)/2) mask(i,j)=22
      forall (j=(nd+1):(nd+1), i=-(nd+1)/2:(nd+1)/2) mask(i,j)=22

c     Set mask=3 for vertices 
      
      mask(-(nd+1),-(nd+1))=3
      mask(-(nd+1),(nd+1))=3
      mask((nd+1),(nd+1))=3
      mask((nd+1),-(nd+1))=3

      do i=-(nd+1), (nd+1)
        do j=-(nd+1), (nd+1)
              print *,mask(i,j)
      enddo
      enddo    
              
      end
  
