subroutine angvel(xx,yy,omg,domgdr,domgdm)

c===========================================================================*
c   SUN AND LATE-TYPE STAR BABCOCK-LEIGHTON DYNAMOS:                  
c
c   computes angular velocity and derivatives at r/R=yy, cos(theta)=xx
c
c   Differential rotation is solar-like, with a "tachocline" centered
c   on the core-envelope interface (r/R=rzc), with half-width w/R=ybl
c   See 2.1 in DC99
c
c   omg is angular velocity [Omega(r,mu) ]
c   domgdr=d Omega/d r       [radial derivative of Omega ]
c   domgdm=d Omega/d mu      [latitudinal derivative of Omega ]
c===========================================================================
c   Refs: Dikpati & Charbonneau 1999 , ApJ, 518, 508-ff
c         Charbonneau et al. 1999, ApJ, 
c
c   Last revised 10/07/2001
c===========================================================================
      implicit  none
      include  "arrays.global.f"
c Input
      real    xx,yy
c Output
      real    omg,domgdr,domgdm
c Local
      real pi,a2,a4,omgc,omge,tt,del,frac,dfrac
      parameter(pi=3.1415926536)
c -- these are for the differential rotation pattern
      parameter(a2=0.1264,a4=0.1591)
      real     mu02,mu04,mu2,mu4
      parameter(mu02=(0.5735764)**2,mu04=mu02*mu02)
      parameter(omgc=1.-a2*mu02-a4*mu04)
      real     erf

c---------- now angular velocity; solar-like parameterization
      if(yy.ge.1.)then
c --> exterior
         omg=0.
	 domgdr=0.
	 domgdm=0.

      else
c --> interior
         mu2=xx*xx
         mu4=mu2*mu2
         omge=1.-a2*mu2-a4*mu4
         del=omge-omgc
         tt=(yy-rzc)/ybl
	 frac=0.5*( 1.+erf(tt) )
	 dfrac=1./(ybl*sqrt(pi)) *exp(-tt**2)

	 omg = omgc+frac*del
         domgdr = dfrac*del                    
         domgdm = -frac*(2.*a2*xx+4.*a4*mu2*xx) 
      endif

      return
      end