subroutine gmiddle(go,gn,x,q,p,nd,nx,dd,dx,du,i)
      double precision go(-(nd+1):(nd+1),-(nd+1):(nd+1),0:nx,0:1)
      double precision gn(-(nd+1):(nd+1),-(nd+1):(nd+1),0:nx,0:1)  
      double precision x(0:nx)
      double precision q(-(nd+1):(nd+1))
      double precision p(-(nd+1):(nd+1))
      integer nd,nx
      double precision dd, dx, du
      integer i

      integer D, k,l
      double precision dold, dnew, d1, d2
      double precision xp, xq, xr, xs, xc
      double precision rp, rq, rr, rs 
      double precision gp(-nd:nd,-nd:nd), gq(-nd:nd,-nd:nd),
     &                 gr(-nd:nd,-nd:nd), gs(-nd:nd,-nd:nd) 
      double precision l2(-nd:nd,-nd:nd), l2im1(-nd:nd,-nd:nd),
     &                  l2i(-nd:nd,-nd:nd), factor(-nd:nd,-nd:nd),
     &                  integral

c     Calculate the characteristic deviation
 
      dold = 0.25 * du * (1. - x(i-1)) ** 2
      dnew = 0.25 * du * (1 - x(i)) ** 2

      d1 = dold/dx
      d2 = dnew/dx
 
      xp = x(i-1) - dold 
      xq = x(i) - dnew
      xr = x(i-1) + dold
      xs = x(i)  + dnew

      rp = xp / (1. - xp)
      rr = xr / (1. - xr)
      rq = xq / (1. - xq)
      rs = xs / (1. - xs)

      xc = 0.5 * (xp + xs)

c     integral = \int_\sigma du dr 1/r^2

      integral = 2. * log((rq * rr)/ (rp * rs))

c     Calculate gn for both patches independently

      do D = 0,1

c     Calculate the laplacian

      forall (l=-nd:nd,k=-nd:nd) factor(l,k) = -0.25 * (1 + q(k) ** 2
     &                                        + p(l) ** 2) ** 2
       
      forall (l=-nd:nd,k=-nd:nd) l2im1(l,k) = factor(l,k) * 
     &                                        ( 1. / dd ** 2) *
     &  (gn(k+1,l,i-1,D)
     &                                         + gn(k-1,l,i-1,D) 
     &                                         + gn(k,l+1,i-1,D)
     &                                         + gn(k,l-1,i-1,D)
     &                                         - 4. * gn(k,l,i-1,D))

      forall (l=-nd:nd,k=-nd:nd)  l2i(l,k) =  factor(l,k) * (1. / dd
     &                                         ** 2) * (go(k+1,l,i,D)
     &                                         + go(k-1,l,i,D) 
     &                                         + go(k,l+1,i,D)
     &                                         + go(k,l-1,i,D)
     &                                         - 4. * go(k,l,i,D))

     
      forall (l=-nd:nd,k=-nd:nd) l2(l,k) = 0.5 * (l2im1(l,k) / x(i-1) **
     &                                     2 + l2i(l,k) / x(i) ** 2)


c    interpolate the g = r*f value


      forall (l=-nd:nd,k=-nd:nd) gp(l,k) = (1. - d1) * gn(k,l,i-1,D) +
     &                                     d1 * gn(k,l,i-2,D)
   
      forall (l=-nd:nd,k=-nd:nd) gr(l,k) = (1. - d1) * go(k,l,i-1,D) +
     &                                     d1 * go(k,l,i,D)
    
      forall (l=-nd:nd,k=-nd:nd) gs(l,k) = (1. -d2) * go(k,l,i,D) +
     &                                     d2 * go(k,l,i+1,D)

      
c    calculate the gq box

      forall (l=-nd:nd,k=-nd:nd) gq(l,k) = gp(l,k) + gs(l,k) - gr(l,k)
     &                                     - 0.5 * integral * l2(l,k) * 
     &                                      xc ** 2
    

c    extrapolate to obtain the new gamma grid value

      forall (l=-nd:nd,k=-nd:nd) gn(k,l,i,D) = 
     &                           (gq(l,k) - gn(k,l,i-1,D) * d2) /
     &                           (1. - d2)

      end do 
      return 
   
       end