subroutine srcpol(xx,yy,alpxy)

c===========================================================================*
c   SUN AND LATE-TYPE STAR BABCOCK-LEIGHTON DYNAMOS:                  
c
c   computes spatial dependency of nonlinear poloidal source term
c   at r/R=yy, cos(theta)=xx;
c
c   The source term is concentrated immediately beneath surface
c   See sec 2.3 and eq. (7) in DC99. It is computed only in 0.92<r/R<1
c   to avoid underflow in the error function calculation
c===========================================================================
c   Refs: Dikpati & Charbonneau 1999 , ApJ, 518, 508-ff
c
c   Last revised 10/07/2001
c===========================================================================
      implicit  none
      include  "arrays.global.f"
c Input
      real    xx,yy
c Output
      real    alpxy
c Local
      real pi,rs2,rs3,ds2,ds3
      parameter(pi=3.1415926536)
c -- these are for the Babcock-Leighton source term
      parameter(rs2=0.95,rs3=1.0,ds2=0.01,ds3=0.01)
      real tt2,tt3,g2r,g3r,thta,sinth,ftheta
      real erf
c---------------------------------------------------------------------------

      if(yy.gt.0.92.and.yy.le.1.)then
c --> below surface: Babcock-Leighton source (and avoiding IEEE underflows)

         tt2=(yy-rs2)/(ds2)
         g2r=0.5*(1.+erf(tt2))
         tt3=(yy-rs3)/(ds3)
c        g3r=0.5*(1.-erf(tt3))
c --> following line to be consistent with DC99
         g3r=    (1.-erf(tt3))

         thta=acos(xx)
         sinth=sin(thta)
         ftheta=xx*sinth

         alpxy=g2r*g3r *ftheta
      else
c --> exterior and deep interior: no source
	 alpxy=0.

      endif

      return
      end