
      subroutine interpolate(gn, q, p, mask, nx, dd, i)
      integer          nd, nx, i
      parameter        (nd=8)
      double precision gn(-(nd+1):(nd+1),-(nd+1):(nd+1),0:nx,0:1)
      double precision q(-(nd+1):(nd+1))  
      double precision p(-(nd+1):(nd+1))
      integer          mask(-(nd+1):(nd+1),-(nd+1):(nd+1))
      double precision dd

      double precision z
      integer          b1,b2,b3,b4
      parameter       (z=1.)
      parameter       (b1=-(nd+1),b2=(nd+1),b3=-(nd+1)/2,b4=(nd+1)/2)  
      integer          D, U, k, l
      integer          jg(-(nd+1):(nd+1),-(nd+1):(nd+1))
      integer          ig(-(nd+1):(nd+1),-(nd+1):(nd+1))
      integer          ig2(-(nd+1):(nd+1),-(nd+1):(nd+1))
      integer          jg2(-(nd+1):(nd+1),-(nd+1):(nd+1))
      integer          k1(-(nd+1):(nd+1),-(nd+1):(nd+1))
      integer          k2(-(nd+1):(nd+1),-(nd+1):(nd+1))
      integer          l1(-(nd+1):(nd+1),-(nd+1):(nd+1))
      integer          l2(-(nd+1):(nd+1),-(nd+1):(nd+1))
      double precision qf(-(nd+1):(nd+1),-(nd+1):(nd+1))
      double precision pf(-(nd+1):(nd+1),-(nd+1):(nd+1))
      double precision cp(-(nd+1):(nd+1),-(nd+1):(nd+1))
      double precision cq(-(nd+1):(nd+1),-(nd+1):(nd+1))
      double precision  cp1(-(nd+1):(nd+1),-(nd+1):(nd+1))
      double precision cp2(-(nd+1):(nd+1),-(nd+1):(nd+1))
      double precision  cq1(-(nd+1):(nd+1),-(nd+1):(nd+1))
      double precision cq2(-(nd+1):(nd+1),-(nd+1):(nd+1))
      double precision q1(-(nd+1):(nd+1),-(nd+1):(nd+1))
      double precision q2(-(nd+1):(nd+1),-(nd+1):(nd+1))
      double precision  p1(-(nd+1):(nd+1),-(nd+1):(nd+1))
      double precision p2(-(nd+1):(nd+1),-(nd+1):(nd+1))
      double precision  v1(-(nd+1):(nd+1),-(nd+1):(nd+1))
      double precision  v2(-(nd+1):(nd+1),-(nd+1):(nd+1)) 
      double precision q1p(-(nd+1):(nd+1),-(nd+1):(nd+1))
      double precision q2p(-(nd+1):(nd+1),-(nd+1):(nd+1))
      double precision p1p(-(nd+1):(nd+1),-(nd+1):(nd+1))
      double precision p2p(-(nd+1):(nd+1),-(nd+1):(nd+1))   
      double precision r1(-(nd+1):(nd+1),-(nd+1):(nd+1))
      double precision r2(-(nd+1):(nd+1),-(nd+1):(nd+1)) 
      double precision r3(-(nd+1):(nd+1),-(nd+1):(nd+1))
      double precision  r4(-(nd+1):(nd+1),-(nd+1):(nd+1)) 
      double precision rx(-(nd+1):(nd+1),-(nd+1):(nd+1))
      double precision  qd1(-(nd+1):(nd+1),-(nd+1):(nd+1)) 
      double precision qd2(-(nd+1):(nd+1),-(nd+1):(nd+1))
      double precision pd1(-(nd+1):(nd+1),-(nd+1):(nd+1))
      double precision pd2(-(nd+1):(nd+1),-(nd+1):(nd+1))      
c     double precision :: q_s2n, p_s2n, cubic
c      double precision :: q_s2n, p_s2n, cubic 

      INTERFACE
       PURE DOUBLE PRECISION FUNCTION Q_S2N(QS,PS)
        DOUBLE PRECISION, INTENT(IN) ::  QS,PS
       END FUNCTION Q_S2N 

       PURE DOUBLE PRECISION FUNCTION P_S2N(QS,PS)
        DOUBLE PRECISION, INTENT(IN) :: QS,PS
       END FUNCTION P_S2N 

       PURE DOUBLE PRECISION FUNCTION CUBIC(yp,yp1,yp2,yp3,yp4,f1
     &,f2,f3,f4)
        DOUBLE PRECISION, INTENT(IN) ::yp,yp1,yp2,yp3,yp4,f1,f2
     &  ,f3,f4 
       END FUNCTION CUBIC 
      END INTERFACE


c     Fourth order interpolation at the patch boundary
c     Upper is the grid we are looking at (rectangular)
c     Lower is the other patch, with circular arc coordinate lines

c     do both grids
  
      do D = 0,1
     
      if (D.eq.0) U = 1
      if (D.eq.1) U = 0
      
c    Find the coordinate of the fictitious point in the lower grid

      forall (l=b1:b2) qf(b1,l)=q_s2n(q(b1),p(l))
      forall (l=b1:b2) pf(b1,l)=p_s2n(q(b1),p(l))
      forall (l=b1:b2) qf(b2,l)=q_s2n(q(b2),p(l))
      forall (l=b1:b2) pf(b2,l)=p_s2n(q(b2),p(l)) 
      forall (k=-nd:nd)qf(k,b1)=q_s2n(q(k),p(b1))
      forall (k=-nd:nd)pf(k,b1)=p_s2n(q(k),p(b1))
      forall (k=-nd:nd)qf(k,b2)=q_s2n(q(k),p(b2))
      forall (k=-nd:nd)pf(k,b2)=p_s2n(q(k),p(b2))

c    All the operations for mask(k,l)=11
  
c    Guard point index and outer neighbor index
      forall(l=b3:b4) ig(b1,l)=int(qf(b1,l)/dd)
      forall(l=b3:b4) jg(b1,l)=int(pf(b1,l)/dd) 
     &                         + int((sign(z,pf(b1,l))-1)/2)  
      forall(l=b3:b4) ig(b2,l)=int(qf(b2,l)/dd)
      forall(l=b3:b4) jg(b2,l)=int(pf(b2,l)/dd) 
     &                          + int((sign(z,pf(b2,l))-1)/2)       

      forall(l=b3:b4) ig2(b1,l)=ig(b1,l) 
     &                          - int(sign(z,q(ig(b1,l))))
      forall(l=b3:b4) jg2(b1,l)=jg(b1,l) 
     &                          - int(sign(z,p(jg(b1,l))))
      forall(l=b3:b4) ig2(b2,l)=ig(b2,l) 
     &                          - int(sign(z,q(ig(b2,l))))
      forall(l=b3:b4) jg2(b2,l)=jg(b2,l) 
     &                          - int(sign(z,p(jg(b2,l))))


c    find the upper grid points used in the cubic interpolation
    
      forall(l=b3:b4) k1(b1,l) = b1 - 1 * int(sign(z,q(b1)))
      forall(l=b3:b4) k2(b1,l) = b1 - 2 * int(sign(z,q(b1)))
      forall(l=b3:b4) l1(b1,l) = b1 - 1 * int(sign(z,p(l)))
      forall(l=b3:b4) l2(b1,l) = b1 - 2 * int(sign(z,p(l)))
      forall(l=b3:b4) k1(b2,l) = b2 - 1 * int(sign(z,q(b2)))
      forall(l=b3:b4) k2(b2,l) = b2 - 2 * int(sign(z,q(b2)))
      forall(l=b3:b4) l1(b2,l) = b2 - 1 * int(sign(z,p(l)))
      forall(l=b3:b4) l2(b2,l) = b2 - 2 * int(sign(z,p(l)))
   
c    calculate the intersection points depending on which combination of
c    crossing lines we use

      forall(l=b3:b4) q1(b1,l) = (1. + sqrt(1. - 4. * p(l)**2 
     &                            * q(ig(b1,l))**2))/
     &                           (2.* q(ig(b1,l)))
      forall(l=b3:b4) p1p(b1,l) = -2.* p(l) * q(ig(b1,l))**2
     &                              /(1.+ sqrt(1.- 4. * p(l)**2
     &                               * q(ig(b1,l))**2))
      forall(l=b3:b4) q2(b1,l) = (1. + sqrt(1. - 4. * p(l)**2 
     &                            * q(ig2(b1,l))**2))/
     &                           (2.* q(ig2(b1,l)))
      forall(l=b3:b4) p2p(b1,l) = -2.* p(l) * q(ig2(b1,l))**2
     &                            /(1. + sqrt(1. - 4.* p(l)**2
     &                             * q(ig2(b1,l))**2))                  
      forall(l=b3:b4) q1(b2,l) = (1. + sqrt(1. - 4. * p(l)**2 
     &                           * q(ig(b2,l))**2))/
     &                           (2.* q(ig(b2,l)))
      forall(l=b3:b4) p1p(b2,l) = -2.* p(l) * q(ig(b2,l))**2
     &                             /(1.+ sqrt(1.- 4. * p(l)**2
     &                              * q(ig(b2,l))**2))
      forall(l=b3:b4) q2(b2,l) = (1. + sqrt(1. - 4. * p(l)**2 
     &                           * q(ig2(b2,l))**2))/
     &                           (2.* q(ig2(b2,l)))
      forall(l=b3:b4) p2p(b2,l) = -2.* p(l) * q(ig2(b2,l))**2
     &                             /(1. + sqrt(1. - 4.* p(l)**2
     &                               * q(ig2(b2,l))**2))

c   Do the lower interpolation depending on which direction is used


      forall(l=b3:b4) v1(b1,l) = cubic(p1p(b1,l),p(jg(b1,l)-1),
     &                           p(jg(b1,l)),p(jg(b1,l)+1),
     &                 p(jg(b1,l)+2),gn(ig(b1,l),jg(b1,l)-1,i,U),
     &                           gn(ig(b1,l),jg(b1,l),i,U),
     &                           gn(ig(b1,l),jg(b1,l)+1,i,U)
     &                   ,gn(ig(b1,l),jg(b1,l)+2,i,U))
      forall(l=b3:b4) v2(b1,l) = cubic(p2p(b1,l),p(jg(b1,l)-1),
     &                            p(jg(b1,l)),p(jg(b1,l)+1)
     &                           ,p(jg(b1,l)+2),gn(ig2(b1,l),
     &                            jg(b1,l)-1,i,U),
     &                           gn(ig2(b1,l),jg(b1,l),i,U),
     &                           gn(ig2(b1,l),jg(b1,l)+1,i,U)
     &                           ,gn(ig2(b1,l),jg(b1,l)+2,i,U))

     
      forall(l=b3:b4) v1(b2,l) = cubic(p1p(b2,l),p(jg(b2,l)-1)
     &                           ,p(jg(b2,l)),p(jg(b2,l)+1),
     &                p(jg(b2,l)+2),gn(ig(b2,l),jg(b2,l)-1,i,U),
     &                           gn(ig(b2,l),jg(b2,l),i,U),
     &                           gn(ig(b2,l),jg(b2,l)+1,i,U)
     &                           ,gn(ig(b2,l),jg(b2,l)+2,i,U))

      forall(l=b3:b4) v2(b2,l) = cubic(p2p(b2,l),p(jg(b2,l)-1)
     &                           ,p(jg(b2,l)),p(jg(b2,l)+1)
     &                   ,p(jg(b2,l)+2),gn(ig2(b2,l),jg(b2,l)-1,i,U),
     &                           gn(ig2(b2,l),jg(b2,l),i,U),
     &                           gn(ig2(b2,l),jg(b2,l)+1,i,U)
     &                           ,gn(ig2(b2,l),jg(b2,l)+2,i,U))

c   Do the upper interpolation depending on which direction is used
 
      forall(l=b3:b4) gn(b1,l,i,D) = cubic(q(b1),q(k2(b1,l)),
     &                              q(k1(b1,l)),q1(b1,l),q2(b1,l)
     &                               ,gn(k2(b1,l),l,i,D),
     &                           gn(k1(b1,l),l,i,D),v1(b1,l),v2(b1,l))
      forall(l=b3:b4) gn(b2,l,i,D) = cubic(q(b2),q(k2(b2,l)),
     &                                q(k1(b2,l)),
     &                                q1(b2,l),q2(b2,l)
     &                          ,gn(k2(b2,l),l,i,D),gn(k1(b2,l),l,i,D),
     &                                v1(b2,l),v2(b2,l))      


c   All the operations for mask(k,l)=12


c   The guard point index and the outer neighbor index

      forall(l=(b4+1):nd) ig(b1,l) = int(qf(b1,l)/dd) 
     &                             + int((sign(z,qf(b1,l))-1)/2) 
      forall(l=(b4+1):nd) jg(b1,l) = int(pf(b1,l)/dd)
      forall(l=(b4+1):nd) ig(b2,l) = int(qf(b2,l)/dd) 
     &                               + int((sign(z,qf(b2,l))-1)/2) 
      forall(l=(b4+1):nd) jg(b2,l) = int(pf(b2,l)/dd) 
 
      forall(l=-nd:(b3-1)) ig(b1,l) = int(qf(b1,l)/dd) 
     &                                + int((sign(z,qf(b1,l))-1)/2)
      forall(l=-nd:(b3-1)) jg(b1,l) = int(pf(b1,l)/dd)
      forall(l=-nd:(b3-1)) ig(b2,l) = int(qf(b2,l)/dd) 
     &                                + int((sign(z,qf(b2,l))-1)/2) 
      forall(l=-nd:(b3-1)) jg(b2,l) = int(pf(b2,l)/dd) 

      forall(l=(b4+1):nd) ig2(b1,l) = ig(b1,l) 
     &                                - int(sign(z,q(ig(b1,l))))
      forall(l=(b4+1):nd) jg2(b1,l) = jg(b1,l) 
     &                                 - int(sign(z,p(jg(b1,l))))      
      forall(l=(b4+1):nd) ig2(b2,l) = ig(b2,l) 
     &                                - int(sign(z,q(ig(b2,l))))
      forall(l=(b4+1):nd) jg2(b2,l) = jg(b2,l) 
     &                                - int(sign(z,p(jg(b2,l))))      

      forall(l=-nd:(b3-1)) ig2(b1,l) = ig(b1,l) 
     &                                - int(sign(z,q(ig(b1,l))))
      forall(l=-nd:(b3-1)) jg2(b1,l) = jg(b1,l) 
     &                                - int(sign(z,p(jg(b1,l))))      
      forall(l=-nd:(b3-1)) ig2(b2,l) = ig(b2,l) 
     &                                - int(sign(z,q(ig(b2,l))))
      forall(l=-nd:(b3-1)) jg2(b2,l) = jg(b2,l) 
     &                                - int(sign(z,p(jg(b2,l))))      


c   Find the upper grid grid points used in the cubic interpolation

      forall(l=(b4+1):nd) k1(b1,l) = b1 - 1 * int(sign(z,q(b1)))
      forall(l=(b4+1):nd) k2(b1,l) = b1 - 2 * int(sign(z,q(b1)))
      forall(l=(b4+1):nd) k1(b2,l) = b2 - 1 * int(sign(z,q(b2)))
      forall(l=(b4+1):nd) k2(b2,l) = b2 - 2 * int(sign(z,q(b2))) 

      forall(l=-nd:(b3-1)) k1(b1,l) = b1 - 1 * int(sign(z,q(b1)))
      forall(l=-nd:(b3-1)) k2(b1,l) = b1 - 2 * int(sign(z,q(b1)))
      forall(l=-nd:(b3-1)) k1(b2,l) = b2 - 1 * int(sign(z,q(b2)))
      forall(l=-nd:(b3-1)) k2(b2,l) = b2 - 2 * int(sign(z,q(b2))) 

c   Calculate the intersection points depending on which combination of
c   crossing lines we use 

      forall(l=(b4+1):nd) q1(b1,l) = sign(z,q(b1)) 
     &                    * sqrt(-p(l)/p(jg(b1,l)) - p(l)**2)
      forall(l=(b4+1):nd) q1p(b1,l) = - (p(jg(b1,l))
     &                      /p(l)) * q1(b1,l)
      forall(l=(b4+1):nd) q2(b1,l) = sign(z,q(b1)) 
     &                    * sqrt(-p(l)/p(jg2(b1,l)) - q(b1)**2)
      forall(l=(b4+1):nd) q2p(b1,l) = - (p(jg2(b1,l))
     &                               /p(l)) * q2(b1,l)
      forall(l=(b4+1):nd) q1(b2,l) = sign(z,q(b2)) 
     &                      * sqrt(-p(l)/p(jg(b2,l)) - p(l)**2)
      forall(l=(b4+1):nd) q1p(b2,l) = - (p(jg(b2,l))
     &                           /p(l)) * q1(b2,l)
      forall(l=(b4+1):nd) q2(b2,l) = sign(z,q(b2)) 
     &                    * sqrt(-p(l)/p(jg2(b2,l)) - q(b2)**2)
      forall(l=(b4+1):nd) q2p(b2,l) = - (p(jg2(b2,l))
     &                              /p(l)) * q2(b2,l)

      forall(l=-nd:(b3-1)) q1(b1,l) = sign(z,q(b1)) 
     &                  * sqrt(-p(l)/p(jg(b1,l)) - p(l)**2)
      forall(l=-nd:(b3-1)) q1p(b1,l) = - (p(jg(b1,l))
     &                     /p(l)) * q1(b1,l)
      forall(l=-nd:(b3-1)) q2(b1,l) = sign(z,q(b1)) 
     &                     * sqrt(-p(l)/p(jg2(b1,l)) - q(b1)**2)
      forall(l=-nd:(b3-1)) q2p(b1,l) = - (p(jg2(b1,l))
     &                       /p(l)) * q2(b1,l)
      forall(l=-nd:(b3-1)) q1(b2,l) = sign(z,q(b2)) 
     &                      * sqrt(-p(l)/p(jg(b2,l)) - p(l)**2)
      forall(l=-nd:(b3-1)) q1p(b2,l) = - (p(jg(b2,l))
     &                        /p(l)) * q1(b2,l)
      forall(l=-nd:(b3-1)) q2(b2,l) = sign(z,q(b2)) * sqrt(-p(l)
     &                        /p(jg2(b2,l)) - q(b2)**2)
      forall(l=-nd:(b3-1)) q2p(b2,l) = - (p(jg2(b2,l))
     &                        /p(l)) * q2(b2,l)

c   Do the lower interpolation depending on which direction is used

      forall(l=(b4+1):nd) v1(b1,l) = cubic (q1p(b1,l),q(ig(b1,l)-1),
     &                                 q(ig(b1,l)),q(ig(b1,l)+1)
     &                               ,q(ig(b1,l)+2),
     &                                gn(ig(b1,l)-1,jg(b1,l),i,U)
     &                               ,gn(ig(b1,l),jg(b1,l),i,U),
     &                                gn(ig(b1,l)+1,jg(b1,l),i,U)
     &                               ,gn(ig(b1,l)+2,jg(b1,l),i,U))

      forall(l=(b4+1):nd) v2(b1,l) = cubic (q2p(b1,l),q(ig(b1,l)-1),
     &                               q(ig(b1,l)),q(ig(b1,l)+1)
     &                     ,q(ig(b1,l)+2),gn(ig(b1,l)-1,jg2(b1,l),i,U)
     &                               ,gn(ig(b1,l),jg2(b1,l),i,U),
     &                                gn(ig(b1,l)+1,jg2(b1,l),i,U),
     &                                gn(ig(b1,l)+2,jg2(b1,l),i,U))

   
      forall(l=-nd:(b3-1)) v1(b1,l) = cubic (q1p(b1,l),q(ig(b1,l)-1),
     &                               q(ig(b1,l)),q(ig(b1,l)+1)
     &                      ,q(ig(b1,l)+2),gn(ig(b1,l)-1,jg(b1,l),i,U)
     &                               ,gn(ig(b1,l),jg(b1,l),i,U),
     &                                gn(ig(b1,l)+1,jg(b1,l),i,U)
     &                               ,gn(ig(b1,l)+2,jg(b1,l),i,U))

      forall(l=-nd:(b3-1)) v2(b1,l) = cubic (q2p(b1,l),q(ig(b1,l)-1),
     &                                q(ig(b1,l)),q(ig(b1,l)+1)
     &                   ,q(ig(b1,l)+2),gn(ig(b1,l)-1,jg2(b1,l),i,U)
     &                               ,gn(ig(b1,l),jg2(b1,l),i,U),
     &                                gn(ig(b1,l)+1,jg2(b1,l),i,U),
     &                                gn(ig(b1,l)+2,jg2(b1,l),i,U))

      forall(l=(b4+1):nd) v1(b2,l) = cubic (q1p(b2,l),q(ig(b2,l)-1),
     &                                q(ig(b2,l)),q(ig(b2,l)+1)
     &                               ,q(ig(b2,l)+2),
     &                                gn(ig(b2,l)-1,jg(b2,l),i,U)
     &                               ,gn(ig(b2,l),jg(b2,l),i,U),
     &                                gn(ig(b2,l)+1,jg(b2,l),i,U)
     &                               ,gn(ig(b2,l)+2,jg(b2,l),i,U))

      forall(l=(b4+1):nd) v2(b2,l) = cubic (q2p(b2,l),q(ig(b2,l)-1),
     &                                q(ig(b2,l)),q(ig(b2,l)+1)
     &                               ,q(ig(b2,l)+2),
     &                                gn(ig(b2,l)-1,jg2(b2,l),i,U)
     &                               ,gn(ig(b2,l),jg2(b2,l),i,U),
     &                                gn(ig(b2,l)+1,jg2(b2,l),i,U),
     &                                gn(ig(b2,l)+2,jg2(b2,l),i,U))

   
      forall(l=-nd:(b3-1)) v1(b2,l) = cubic (q1p(b2,l),q(ig(b2,l)-1),
     &                                q(ig(b2,l)),q(ig(b2,l)+1)
     &                               ,q(ig(b2,l)+2),
     &                                gn(ig(b2,l)-1,jg(b2,l),i,U)
     &                               ,gn(ig(b2,l),jg(b2,l),i,U),
     &                                gn(ig(b2,l)+1,jg(b2,l),i,U)
     &                               ,gn(ig(b2,l)+2,jg(b2,l),i,U))

      forall(l=-nd:(b3-1)) v2(b2,l) = cubic (q2p(b2,l),q(ig(b2,l)-1),
     &                                q(ig(b2,l)),q(ig(b2,l)+1)
     &                      ,q(ig(b2,l)+2),gn(ig(b2,l)-1,jg2(b2,l),i,U)
     &                               ,gn(ig(b2,l),jg2(b2,l),i,U),
     &                                gn(ig(b2,l)+1,jg2(b2,l),i,U),
     &                                gn(ig(b2,l)+2,jg2(b2,l),i,U))


c   Do the upper interpolation depending on which direction is used
    
      forall(l=(b4+1):nd) gn(b1,l,i,D) = cubic(q(b1),q(k2(b1,l)),
     &                                q(k1(b1,l)),q1(b1,l),q2(b1,l)
     &                               ,gn(k2(b1,l),l,i,D),
     &                       gn(k1(b1,l),l,i,D),v1(b1,l),v2(b1,l))
      forall(l=(b4+1):nd) gn(b2,l,i,D) = cubic(q(b2),q(k2(b2,l)),
     &                                q(k1(b2,l)),q1(b2,l),q2(b2,l)
     &                               ,gn(k2(b2,l),l,i,D),
     &                         gn(k1(b2,l),l,i,D),v1(b2,l),v2(b2,l))
      forall(l=-nd:(b3-1)) gn(b1,l,i,D) = cubic(q(b1),q(k2(b1,l)),
     &                                q(k1(b1,l)),q1(b1,l),q2(b1,l)
     &                               ,gn(k2(b1,l),l,i,D),
     &                          gn(k1(b1,l),l,i,D),v1(b1,l),v2(b1,l))
      forall(l=-nd:(b3-1)) gn(b2,l,i,D) = cubic(q(b2),q(k2(b2,l)),
     &                                q(k1(b2,l)),q1(b2,l),q2(b2,l)
     &                               ,gn(k2(b2,l),l,i,D),
     &                           gn(k1(b2,l),l,i,D),v1(b2,l),v2(b2,l))


c   All the operations for mask(k,l)=21

c   The guard point index and the outer neighbor index
     
      forall(k=(b4+1):nd) ig(k,b1) = int(qf(k,b1)/dd)
      forall(k=(b4+1):nd) jg(k,b1) = int(pf(k,b1)/dd) 
     &                               + int((sign(z,pf(k,b1))-1)/2) 
      forall(k=(b4+1):nd) ig(k,b2) = int(qf(k,b2)/dd)
      forall(k=(b4+1):nd) jg(k,b2) = int(pf(k,b2)/dd) 
     &                               + int((sign(z,pf(k,b2))-1)/2) 
      forall(k=-nd:(b3-1)) ig(k,b1) = int(qf(k,b1)/dd)
      forall(k=-nd:(b3-1)) jg(k,b1) = int(pf(k,b1)/dd) 
     &                               + int((sign(z,pf(k,b1))-1)/2) 
      forall(k=-nd:(b3-1)) ig(k,b2) = int(qf(k,b2)/dd)
      forall(k=-nd:(b3-1)) jg(k,b2) = int(pf(k,b2)/dd) 
     &                               + int((sign(z,pf(k,b2))-1)/2) 

      forall(k=(b4+1):nd) ig2(k,b1) = ig(k,b1) 
     &                                - int(sign(z,q(ig(k,b1))))
      forall(k=(b4+1):nd) jg2(k,b1) = jg(k,b1) 
     &                                - int(sign(z,p(jg(k,b1))))     
      forall(k=(b4+1):nd) ig2(k,b2) = ig(k,b2) 
     &                                - int(sign(z,q(ig(k,b2))))
      forall(k=(b4+1):nd) jg2(k,b2) = jg(k,b2) 
     &                                - int(sign(z,p(jg(k,b2))))
      forall(k=-nd:(b3-1)) ig2(k,b1) = ig(k,b1) 
     &                                - int(sign(z,q(ig(k,b1)))) 
      forall(k=-nd:(b3-1)) jg2(k,b1) = jg(k,b1) 
     &                                - int(sign(z,p(jg(k,b1))))     
      forall(k=-nd:(b3-1)) ig2(k,b2) = ig(k,b2) 
     &                                - int(sign(z,q(ig(k,b2))))
      forall(k=-nd:(b3-1)) jg2(k,b2) = jg(k,b2) 
     &                                - int(sign(z,p(jg(k,b2))))


c   Find the upper grid points used in the cubic interpolation
   
      forall(k=(b4+1):nd) l1(k,b1) = b1 - 1 * int(sign(z,p(b1)))
      forall(k=(b4+1):nd) l2(k,b1) = b1 - 2 * int(sign(z,p(b1)))
      forall(k=(b4+1):nd) l1(k,b2) = b2 - 1 * int(sign(z,p(b2)))
      forall(k=(b4+1):nd) l2(k,b2) = b2 - 2 * int(sign(z,p(b2)))
      forall(k=-nd:(b3-1)) l1(k,b1) = b1 - 1 * int(sign(z,p(b1)))
      forall(k=-nd:(b3-1)) l2(k,b1) = b1 - 2 * int(sign(z,p(b1)))
      forall(k=-nd:(b3-1)) l1(k,b2) = b2 - 1 * int(sign(z,p(b2)))
      forall(k=-nd:(b3-1)) l2(k,b2) = b2 - 2 * int(sign(z,p(b2)))

c   Calculate the intersection points depending on which combination of crossing 
c   lines we use

           
      forall(k=(b4+1):nd) p1(k,b1) = sign(z,p(b1)) * sqrt(q(k)
     &                               /q(ig(k,b1)) - q(k)**2)
      forall(k=(b4+1):nd) p1p(k,b1) = -(q(ig(k,b1))
     &                               /q(k)) * p1(k,b1)
      forall(k=(b4+1):nd) p2(k,b1) = sign(z,p(b1)) * sqrt(q(k)
     &                               /q(ig2(k,b1)) - q(k)**2)
      forall(k=(b4+1):nd) p2p(k,b1) = -(q(ig2(k,b1))/q(k)) * p2(k,b1)
      forall(k=(b4+1):nd) p1(k,b2) = sign(z,p(b2)) * sqrt(q(k)
     &                              /q(ig(k,b2)) - q(k)**2)
      forall(k=(b4+1):nd) p1p(k,b2) = -(q(ig(k,b2))/q(k)) 
     &                                * p1(k,b2)
      forall(k=(b4+1):nd) p2(k,b2) = sign(z,p(b2)) 
     &                            * sqrt(q(k)/q(ig2(k,b2)) - q(k)**2)
      forall(k=(b4+1):nd) p2p(k,b2) = -(q(ig2(k,b2))/q(k)) * p2(k,b2)   

      forall(k=-nd:(b3-1)) p1(k,b1) = sign(z,p(b1)) 
     &                             * sqrt(q(k)/q(ig(k,b1)) - q(k)**2)
      forall(k=-nd:(b3-1)) p1p(k,b1) = -(q(ig(k,b1))/q(k)) * p1(k,b1)
      forall(k=-nd:(b3-1)) p2(k,b1) = sign(z,p(b1)) 
     &                            * sqrt(q(k)/q(ig2(k,b1)) - q(k)**2)
      forall(k=-nd:(b3-1)) p2p(k,b1) = -(q(ig2(k,b1))/q(k)) * p2(k,b1)
      forall(k=-nd:(b3-1)) p1(k,b2) = sign(z,p(b2)) 
     &                            * sqrt(q(k)/q(ig(k,b2)) - q(k)**2)
      forall(k=-nd:(b3-1)) p1p(k,b2) = -(q(ig(k,b2))/q(k)) * p1(k,b2)
      forall(k=-nd:(b3-1)) p2(k,b2) = sign(z,p(b2)) 
     &                            * sqrt(q(k)/q(ig2(k,b2)) - q(k)**2)
      forall(k=-nd:(b3-1)) p2p(k,b2) = -(q(ig2(k,b2))/q(k)) * p2(k,b2)   

c     Do the lower interpolation depending on which direction is used
   
      forall(k=(b4+1):nd) v1(k,b1) = cubic(p1p(k,b1),p(jg(k,b1)-1),
     &                               p(jg(k,b1)),p(jg(k,b1)+1)
     &                               ,p(jg(k,b1)+2),
     &                                gn(ig(k,b1),jg(k,b1)-1,i,U),
     &                               gn(ig(k,b1),jg(k,b1),i,U),
     &                               gn(ig(k,b1),jg(k,b1)+1,i,U),
     &                               gn(ig(k,b1),jg(k,b1)+2,i,U))         
      forall(k=(b4+1):nd) v2(k,b1) = cubic(p2p(k,b1),p(jg(k,b1)-1),
     &                                p(jg(k,b1)),p(jg(k,b1)+1)
     &                               ,p(jg(k,b1)+2),
     &                                gn(ig2(k,b1),jg(k,b1)-1,i,U),
     &                               gn(ig2(k,b1),jg(k,b1),i,U),
     &                               gn(ig2(k,b1),jg(k,b1)+1,i,U),
     &                               gn(ig2(k,b1),jg(k,b1)+2,i,U))

      forall(k=(b4+1):nd) v1(k,b2) = cubic(p1p(k,b2),p(jg(k,b2)-1),
     &                               p(jg(k,b2)),p(jg(k,b2)+1)
     &                               ,p(jg(k,b2)+2),
     &                                gn(ig(k,b2),jg(k,b2)-1,i,U),
     &                               gn(ig(k,b2),jg(k,b2),i,U),
     &                               gn(ig(k,b2),jg(k,b2)+1,i,U),
     &                               gn(ig(k,b2),jg(k,b2)+2,i,U))         
      forall(k=(b4+1):nd) v2(k,b2) = cubic(p2p(k,b2),p(jg(k,b2)-1),
     &                               p(jg(k,b2)),p(jg(k,b2)+1)
     &                               ,p(jg(k,b2)+2),
     &                                gn(ig2(k,b2),jg(k,b2)-1,i,U),
     &                               gn(ig2(k,b2),jg(k,b2),i,U),
     &                               gn(ig2(k,b2),jg(k,b2)+1,i,U),
     &                               gn(ig2(k,b2),jg(k,b2)+2,i,U))

      forall(k=-nd:(b3-1)) v1(k,b1) = cubic(p1p(k,b1),p(jg(k,b1)-1),
     &                               p(jg(k,b1)),p(jg(k,b1)+1)
     &                               ,p(jg(k,b1)+2),
     &                               gn(ig(k,b1),jg(k,b1)-1,i,U),
     &                               gn(ig(k,b1),jg(k,b1),i,U),
     &                               gn(ig(k,b1),jg(k,b1)+1,i,U),
     &                               gn(ig(k,b1),jg(k,b1)+2,i,U))         
      forall(k=-nd:(b3-1)) v2(k,b1) = cubic(p2p(k,b1),p(jg(k,b1)-1),
     &                                p(jg(k,b1)),p(jg(k,b1)+1)
     &                               ,p(jg(k,b1)+2),
     &                               gn(ig2(k,b1),jg(k,b1)-1,i,U),
     &                               gn(ig2(k,b1),jg(k,b1),i,U),
     &                               gn(ig2(k,b1),jg(k,b1)+1,i,U),
     &                               gn(ig2(k,b1),jg(k,b1)+2,i,U))

      forall(k=(b4+1):nd) v1(k,b2) = cubic(p1p(k,b2),p(jg(k,b2)-1),
     &                                p(jg(k,b2)),p(jg(k,b2)+1)
     &                               ,p(jg(k,b2)+2),
     &                                gn(ig(k,b2),jg(k,b2)-1,i,U),
     &                               gn(ig(k,b2),jg(k,b2),i,U),
     &                               gn(ig(k,b2),jg(k,b2)+1,i,U),
     &                               gn(ig(k,b2),jg(k,b2)+2,i,U))         
      forall(k=(b4+1):nd) v2(k,b2) = cubic(p2p(k,b2),p(jg(k,b2)-1),
     &                               p(jg(k,b2)),p(jg(k,b2)+1)
     &                               ,p(jg(k,b2)+2),
     &                                gn(ig2(k,b2),jg(k,b2)-1,i,U),
     &                               gn(ig2(k,b2),jg(k,b2),i,U),
     &                               gn(ig2(k,b2),jg(k,b2)+1,i,U),
     &                               gn(ig2(k,b2),jg(k,b2)+2,i,U))


c     Do the upper interpolation depending on which direction is used

      forall(k=(b4+1):nd) gn(k,b1,i,D) = cubic(p(b1),p(l2(k,b1)),
     &                                    p(l1(k,b1)),p1(k,b1)
     &                                   ,p2(k,b1),gn(k,l2(k,b1),i,D)
     &                                   ,gn(k,l1(k,b1),i,D)
     &                                   ,v1(k,b1),v2(k,b1)) 
      forall(k=(b4+1):nd) gn(k,b2,i,D) = cubic(p(b2),p(l2(k,b2))
     &                                   ,p(l1(k,b2)),p1(k,b2)
     &                                   ,p2(k,b2),gn(k,l2(k,b2),i,D)
     &                                   ,gn(k,l1(k,b2),i,D)
     &                                   ,v1(k,b2),v2(k,b2))

      forall(k=-nd:(b3-1)) gn(k,b1,i,D) = cubic(p(b1),p(l2(k,b1))
     &                                    ,p(l1(k,b1)),p1(k,b1)
     &                                   ,p2(k,b1),gn(k,l2(k,b1),i,D)
     &                                    ,gn(k,l1(k,b1),i,D)
     &                                   ,v1(k,b1),v2(k,b1)) 
      forall(k=-nd:(b3-1)) gn(k,b2,i,D) = cubic(p(b2),p(l2(k,b2))
     &                                    ,p(l1(k,b2)),p1(k,b2)
     &                                   ,p2(k,b2),gn(k,l2(k,b2),i,D)
     &                                   ,gn(k,l1(k,b2),i,D)
     &                                   ,v1(k,b2),v2(k,b2))


c     All the operations for mask(k,l)=22
    
c     The guard point index and the outer neighbor index

      forall(k=b3:b4) ig(k,b1) = int(qf(k,b1)/dd) 
     &                           + int((sign(z,qf(k,b1)) - 1)/2)
      forall(k=b3:b4) jg(k,b1) = int(pf(k,b1)/dd)
      forall(k=b3:b4) ig(k,b2) = int(qf(k,b2)/dd) 
     &                           + int((sign(z,qf(k,b2)) - 1)/2)
      forall(k=b3:b4) jg(k,b2) = int(pf(k,b2)/dd) 

      forall(k=b3:b4) ig2(k,b1) = ig(k,b1) - int(sign(z,q(ig(k,b1))))
      forall(k=b3:b4) jg2(k,b1) = jg(k,b1) - int(sign(z,p(jg(k,b1))))
      forall(k=b3:b4) ig2(k,b2) = ig(k,b2) - int(sign(z,q(ig(k,b2))))
      forall(k=b3:b4) jg2(k,b2) = jg(k,b2) - int(sign(z,p(jg(k,b2))))  

c     Find the upper grid points used in the cubic interpolation

      forall(k=b3:b4) l1(k,b1) = b1 - 1 * int(sign(z,p(b1)))       
      forall(k=b3:b4) l2(k,b1) = b1 - 2 * int(sign(z,p(b1)))
      forall(k=b3:b4) l1(k,b2) = b2 - 1 * int(sign(z,p(b2)))       
      forall(k=b3:b4) l2(k,b2) = b2 - 2 * int(sign(z,p(b2))) 

c     Calculate the intersection points depending on which combination of
c      crossing lines we use

      forall(k=b3:b4) p1(k,b1) = (-1. - sqrt(1. - 4. * q(k)**2 
     &                            * (p(jg(k,b1)))**2))/
     &                           (2.* p(jg(k,b1)))
      forall(k=b3:b4) q1p(k,b1) = 2.*q(k)*(p(jg(k,b1)))**2
     &                            /(1. + sqrt(1. - 4. * q(k)
     &                            **2 * (p(jg(k,b1)))**2))
      forall(k=b3:b4) p2(k,b1) = (-1. - sqrt(1. - 4. * q(k)**2 
     &                            * p(jg2(k,b1))**2))/
     &                            (2.*p(jg2(k,b1)))
      forall(k=b3:b4) q2p(k,b1) = 2.*q(k)*(p(jg2(k,b1)))**2
     &                            /(1. +sqrt(1. - 4. * 
     &                            q(k)**2 * (p(jg2(k,b1)))**2))
    
      forall(k=b3:b4) p1(k,b2) = (-1. - sqrt(1. - 4. * q(k)**2 
     &                           * p(jg(k,b2))**2))/
     &                           (2.* p(jg(k,b2)))
      forall(k=b3:b4) q1p(k,b2) = 2.*q(k)*(p(jg(k,b2)))**2
     &                            /(1. + sqrt(1. - 4. * q(k)
     &                            **2 * p(jg(k,b2))**2))
      forall(k=b3:b4) p2(k,b2) = (-1. - sqrt(1. - 4. * q(k)**2 
     &                            * p(jg2(k,b2))**2))/
     &                            (2.*p(jg2(k,b2)))
      forall(k=b3:b4) q2p(k,b2) = 2.*q(k)*p(jg2(k,b2))**2
     &                            /(1. +sqrt(1. - 4. * 
     &                            q(k)**2 * p(jg2(k,b2))**2))

c    Do the lower interpolation depending on which direction we use


      forall(k=b3:b4) v1(k,b1) = cubic(q1p(k,b1), q(ig(k,b1)-1), 
     &                            q(ig(k,b1))
     &                           , q(ig(k,b1)+1),q(ig(k,b1)+2)
     &                           ,gn(ig(k,b1)-1,jg(k,b1),i,U), 
     &                             gn(ig(k,b1),jg(k,b1),i,U)
     &                           ,gn(ig(k,b1)+1,jg(k,b1),i,U),
     &                             gn(ig(k,b1)+2,jg(k,b1),i,U))     
      forall(k=b3:b4) v1(k,b2) = cubic(q1p(k,b2), q(ig(k,b2)-1), 
     &                           q(ig(k,b2))
     &                           , q(ig(k,b2)+1),q(ig(k,b2)+2)
     &                           ,gn(ig(k,b2)-1,jg(k,b2),i,U), 
     &                            gn(ig(k,b2),jg(k,b2),i,U)
     &                           ,gn(ig(k,b2)+1,jg(k,b2),i,U)
     &                           ,gn(ig(k,b2)+2,jg(k,b2),i,U))     
      forall(k=b3:b4) v2(k,b1) = cubic(q2p(k,b1), q(ig(k,b1)-1), 
     &                            q(ig(k,b1))
     &                           , q(ig(k,b1)+1),q(ig(k,b1)+2)
     &                           ,gn(ig(k,b1)-1,jg2(k,b1),i,U), 
     &                            gn(ig(k,b1),jg2(k,b1),i,U)
     &                           ,gn(ig(k,b1)+1,jg2(k,b1),i,U),
     &                            gn(ig(k,b1)+2,jg2(k,b1),i,U))     
      forall(k=b3:b4) v2(k,b2) = cubic(q2p(k,b2), q(ig(k,b2)-1), 
     &                            q(ig(k,b2))
     &                           , q(ig(k,b2)+1),q(ig(k,b2)+2)
     &                           ,gn(ig(k,b2)-1,jg2(k,b2),i,U), 
     &                             gn(ig(k,b2),jg2(k,b2),i,U)
     &                           ,gn(ig(k,b2)+1,jg2(k,b2),i,U),
     &                            gn(ig(k,b2)+2,jg2(k,b2),i,U))     



c     Do the upper interpolation depending on which direction is used

      forall(k=b3:b4) gn(k,b1,i,D) = cubic(p(b1), p(l2(k,b1))
     &                               , p(l1(k,b1)),p1(k,b1)
     &                               ,p2(k,b1),gn(k,l2(k,b1),i,D), 
     &                                gn(k,l1(k,b1),i,D)
     &                                ,v1(k,b1),v2(k,b1))       

      forall(k=b3:b4) gn(k,b2,i,D) = cubic(p(b2), p(l2(k,b2)), 
     &                               p(l1(k,b2)),p1(k,b2)
     &                               ,p2(k,b2),gn(k,l2(k,b2),i,D), 
     &                                gn(k,l1(k,b2),i,D)
     &                               ,v1(k,b2),v2(k,b2))

c     All the operations for mask(k,l)=3

c     guard point index and the outer neighbor index

      ig(b1,b1)=int(qf(b1,b1)/dd)
      ig(b1,b2)=int(qf(b1,b2)/dd)
      ig(b2,b1)=int(qf(b2,b1)/dd)
      ig(b2,b2)=int(qf(b2,b2)/dd)
     
      jg(b1,b1)=int(qf(b1,b1)/dd)
      jg(b1,b2)=int(qf(b1,b2)/dd)
      jg(b2,b1)=int(qf(b2,b1)/dd)
      jg(b2,b2)=int(qf(b2,b2)/dd)

      ig2(b1,b1)= ig(b1,b1) - int(sign (z,q(ig(b1,b1))))
      ig2(b1,b2)= ig(b1,b2) - int(sign (z,q(ig(b1,b2))))         
      ig2(b2,b1)= ig(b2,b1) - int(sign (z,q(ig(b2,b1))))
      ig2(b2,b2)= ig(b2,b2) - int(sign (z,q(ig(b2,b2))))

      jg2(b1,b2)= jg(b1,b1) - int(sign (z,p(ig(b1,b1))))
      jg2(b1,b2)= jg(b1,b2) - int(sign (z,p(ig(b1,b2))))
      jg2(b2,b1)= jg(b2,b1) - int(sign (z,p(ig(b2,b1))))
      jg2(b2,b2)= jg(b2,b2) - int(sign (z,p(ig(b2,b2))))


c     Find the upper grid points used in the cubic interpolation

      k1(b1,b1) = b1 - 1 * int(sign(z,q(b1)))
      k1(b1,b2) = b1 - 1 * int(sign(z,q(b1)))       
      k1(b2,b1) = b2 - 1 * int(sign(z,q(b2)))
      k1(b2,b2) = b2 - 1 * int(sign(z,q(b2)))
      
      k2(b1,b1) = b1 - 2 * int(sign(z,q(b1)))
      k2(b1,b2) = b1 - 2 * int(sign(z,q(b1)))       
      k2(b2,b1) = b2 - 2 * int(sign(z,q(b2)))
      k2(b2,b2) = b2 - 2 * int(sign(z,q(b2)))      

      l1(b1,b1) = b1 - 1 * int(sign(z,p(b1)))
      l1(b1,b2) = b2 - 1 * int(sign(z,p(b2)))     
      l1(b2,b1) = b1 - 1 * int(sign(z,p(b1)))
      l1(b2,b2) = b2 - 1 * int(sign(z,p(b2)))

      l2(b1,b1) = b1 - 2 * int(sign(z,p(b1)))
      l2(b1,b2) = b2 - 2 * int(sign(z,p(b2)))     
      l2(b2,b1) = b1 - 2 * int(sign(z,p(b1)))
      l2(b2,b2) = b2 - 2 * int(sign(z,p(b2)))


c     Do the lower interpolation depending on which direction is used

      v1(b1,b1) = gn(ig(b1,b1),jg(b1,b1),i,U)
      v1(b1,b2) = gn(ig(b1,b2),jg(b1,b2),i,U)
      v1(b2,b1) = gn(ig(b2,b1),jg(b2,b1),i,U)
      v1(b2,b2) = gn(ig(b2,b2),jg(b2,b2),i,U)

      v1(b1,b1) = gn(ig2(b1,b1),jg2(b1,b1),i,U)
      v1(b1,b2) = gn(ig2(b1,b2),jg2(b1,b2),i,U)
      v1(b2,b1) = gn(ig2(b2,b1),jg2(b2,b1),i,U)
      v1(b2,b2) = gn(ig2(b2,b2),jg2(b2,b2),i,U)

c     Do the upper interpolation depending on which direction is used

      qd1(b1,b1) = q_s2n(q(ig(b1,b1)),p(jg(b1,b1)))
      qd1(b1,b2) = q_s2n(q(ig(b1,b2)),p(jg(b1,b2)))
      qd1(b2,b1) = q_s2n(q(ig(b2,b1)),p(jg(b2,b1)))
      qd1(b2,b2) = q_s2n(q(ig(b2,b2)),p(jg(b2,b2)))

      pd1(b1,b1) = p_s2n(q(ig(b1,b1)),p(jg(b1,b1)))
      pd1(b1,b2) = p_s2n(q(ig(b1,b2)),p(jg(b1,b2)))
      pd1(b2,b1) = p_s2n(q(ig(b2,b1)),p(jg(b2,b1)))
      pd1(b2,b2) = p_s2n(q(ig(b2,b2)),p(jg(b2,b2)))      

      qd2(b1,b1) = q_s2n(q(ig2(b1,b1)),p(jg2(b1,b1)))
      qd2(b1,b2) = q_s2n(q(ig2(b1,b2)),p(jg2(b1,b2)))
      qd2(b2,b1) = q_s2n(q(ig2(b2,b1)),p(jg2(b2,b1)))
      qd2(b2,b2) = q_s2n(q(ig2(b2,b2)),p(jg2(b2,b2)))

      pd2(b1,b1) = p_s2n(q(ig2(b1,b1)),p(jg2(b1,b1)))
      pd2(b1,b2) = p_s2n(q(ig2(b1,b2)),p(jg2(b1,b2)))
      pd2(b2,b1) = p_s2n(q(ig2(b2,b1)),p(jg2(b2,b1)))
      pd2(b2,b2) = p_s2n(q(ig2(b2,b2)),p(jg2(b2,b2)))      

      r1(b1,b1) = sqrt(q(k1(b1,b1))**2+p(l1(b1,b1))**2)
      r1(b1,b2) = sqrt(q(k1(b1,b2))**2+p(l1(b1,b2))**2)     
      r1(b2,b1) = sqrt(q(k1(b2,b1))**2+p(l1(b2,b1))**2)
      r1(b2,b2) = sqrt(q(k1(b2,b2))**2+p(l1(b2,b2))**2)

      r2(b1,b1) = sqrt(q(k2(b1,b1))**2+p(l2(b1,b1))**2)
      r2(b1,b2) = sqrt(q(k2(b1,b2))**2+p(l2(b1,b2))**2)     
      r2(b2,b1) = sqrt(q(k2(b2,b1))**2+p(l2(b2,b1))**2)
      r2(b2,b2) = sqrt(q(k2(b2,b2))**2+p(l2(b2,b2))**2)

      rx(b1,b1) = sqrt(q(b1)**2+p(b1)**2) 
      rx(b1,b2) = sqrt(q(b1)**2+p(b2)**2) 
      rx(b2,b1) = sqrt(q(b2)**2+p(b1)**2) 
      rx(b2,b2) = sqrt(q(b2)**2+p(b2)**2) 

      r3(b1,b1) = sqrt(qd1(b1,b1)**2 + pd1(b1,b1)**2) 
      r3(b1,b2) = sqrt(qd1(b1,b2)**2 + pd1(b1,b2)**2)
      r3(b2,b1) = sqrt(qd1(b2,b1)**2 + pd1(b2,b1)**2)
      r3(b2,b2) = sqrt(qd1(b2,b2)**2 + pd1(b2,b2)**2)

      r4(b1,b1) = sqrt(qd2(b1,b1)**2 + pd2(b1,b1)**2) 
      r4(b1,b2) = sqrt(qd2(b1,b2)**2 + pd2(b1,b2)**2)
      r4(b2,b1) = sqrt(qd2(b2,b1)**2 + pd2(b2,b1)**2)
      r4(b2,b2) = sqrt(qd2(b2,b2)**2 + pd2(b2,b2)**2)

      gn(b1,b1,i,D) = cubic(rx(b1,b1),r1(b1,b1),
     &                      r2(b1,b1),r3(b1,b1),r4(b1,b1)
     &                      ,gn(k1(b1,b1),l1(b1,b1),i,D),
     &                       gn(k2(b1,b1),l2(b1,b1),i,D)
     &                      ,v1(b1,b1),v2(b1,b1))                         

      gn(b1,b2,i,D) = cubic(rx(b1,b2),r1(b1,b2),
     &                      r2(b1,b2),r3(b1,b2),r4(b1,b2)
     &                      ,gn(k1(b1,b2),l1(b1,b2),i,D),
     &                       gn(k2(b1,b2),l2(b1,b2),i,D)
     &                      ,v1(b1,b2),v2(b1,b2))

      gn(b2,b1,i,D) = cubic(rx(b2,b1),r1(b2,b1),
     &                      r2(b2,b1),r3(b2,b1),r4(b2,b1)
     &                      ,gn(k1(b2,b1),l1(b2,b1),i,D),
     &                       gn(k2(b2,b1),l2(b2,b1),i,D)
     &                      ,v1(b2,b1),v2(b2,b1))

      gn(b2,b2,i,D) = cubic(rx(b2,b2),r1(b2,b2),
     &                      r2(b2,b2),r3(b2,b2),r4(b2,b2)
     &                      ,gn(k1(b2,b2),l1(b2,b2),i,D),
     &                       gn(k2(b2,b2),l2(b2,b2),i,D)
     &                      ,v1(b2,b2),v2(b2,b2))


      end do
      return 
      end
 