       subroutine separ(xd)
c
c  seperates a field into hurricane component and remainder
c

       parameter( nmx=24,nmx1=nmx+1,nmx2=nmx*2,nmx6=nmx*6)
       parameter (imx=150, jmx=150)
       dimension xr(nmx),xd(imx,jmx)
cc
c      common  /coor/ xv,yv,xold,yold,xcorn,ycorn,factr,ix,iy
c      common  /ifact/nnn,rovect(nmx),rb,ienv
        common /posit/  xold,yold,xcorn,ycorn,rox,xv,yv,rb,ihour
        common /total/ uh(imx,jmx),vh(imx,jmx),ds(imx,jmx),del(imx,jmx),
     *              tha(imx,jmx),f(imx,jmx),ddel,dtha,dzbdx,dzbdy
       common /xxx/  xf(imx,jmx),xc,yc,dx,dy
c
c new arrays
        dimension b(nmx),w(nmx),ab(nmx,nmx1),ipvt(nmx)
     1      ,wrk(nmx6),iwrk(nmx2)
        common /matrix/ a(nmx,nmx),capd2
        common /vect/rovect(nmx),xvect(nmx),yvect(nmx)
c
        data xr/24*0./
c
c  xc,yc are hurricane coordinates
c  ro  is radius at which hurricane component of field goes to zero
c  xr array contains the field values of 12 equally spaced points
c     on circle of radius ro centered at xc,yc
c
c  set ro to be max value of rovect
c
c
        ro=0.
        do 22 i=1,nmx
        ro=amax1(ro,rovect(i))
22       continue
          pi = 4.*atan(1.0)
       pi180 = 4.*atan(1.0)/180.
       fact =  cos(yold*pi180)
cc
cc   xc is the i position of the center of the old vortex
cc   yc is the j position of the center of the old vortex
cc   ddel is the long. in radians of the outer nest
cc   dtha is the lat.  in radians of the outer nest
cc
c no fact here
c      dx=fact*ddel/pi180
c
       dx=ddel/pi180
       dy=dtha/pi180
cc
cc
       xc = (xold-xcorn)
       yc = (yold-ycorn)
       is=int((xc-ro/fact)/dx) +1.
       ie=int((xc+ro/fact)/dx + 1.)
       js=int((yc-ro)/dy) +1.
       je=int((yc+ro)/dy + 1.)
c
        do 1 j = 1 , jmx
        do 1 i = 1 , imx
          xf(i,j)  = xd(i,j)
1       continue
c
c  subroutine bound computes field values of array xr using
c         bilinear interpolation
c
c 
        call bound(nmx,xr,rovect)
c
c  xrop(nmx) are the interpolated values of the disturbance
c   field at the rovect pts
c
c romax is the maximum value in rovect(nmx). within the loop a local
c ro is computed for use in the separation. at the start of the loop
c ro is again set to romax to define the domain.
c
c
c
        w=0.
        romax=ro
c
        do 10 ix=is,ie
        do 11 jy=js,je
             ro=romax
             x=del(ix,jy)/pi180 -xcorn
             y=tha(ix,jy)/pi180 -ycorn
c
              delx=(x-xc)*fact
              dely=(y-yc)
             dr=sqrt((delx)**2 +(dely)**2)
             if(dr.gt.ro)goto11
               if(delx.ne.0.) theta=atan((dely)/(delx))
               if(delx.eq.0..and.dely.lt.0.)theta=270.*pi180
               if(delx.eq.0..and.dely.gt.0.)theta=90. *pi180
             if(delx.lt.0.)theta=theta+pi
             if(theta.lt.0.)theta=2.*pi+theta
             n1=int(theta*nmx/(2.*pi))
             if(n1.gt.nmx)print*,n1,theta*57.296
             if(n1.lt.0)print*,n1,theta*57.296
             n2=n1+2
             if(n2.gt.nmx)n2=n2-nmx
             delth=theta- 2.*pi*float(n1)/float(nmx)
c
             ro=delth*float(nmx)/(2.*pi)*(rovect(n2)-rovect(n1+1))
     1             +rovect(n1+1)
             if(dr.gt.ro)goto11
             xro=delth*float(nmx)/(2.*pi)*(xr(n2)-xr(n1+1)) +xr(n1+1)
cc
c now add new code to compute distance from each gridpt. to rovect pts
c
             do 12 ip=1,nmx
             dpij= (fact*(x-xvect(ip)))**2 +(y-yvect(ip))**2
             b(ip)=exp(-dpij/capd2)
12           continue
c
c
             do 44 ip=1,nmx
               do 43 jp=1,nmx
43               ab(ip,jp)=a(ip,jp)
               ab(ip,nmx1)=b(ip)
44           continue
c
c solve system using constrained least squares method
c
             call wnnls(ab,nmx,0,nmx,nmx,0,1.,w,rnm,md,iwrk,wrk)
c
             temp=0.
             do 20 ip=1,nmx
             temp=temp +w(ip)*xr(ip)
20           continue
c            xh(ix,jy)=xf(ix,jy)-temp
             xd(ix,jy)=temp
11     continue
10     continue
       return
       end
