subroutine srcpol(xx,yy,alpxy)

c===========================================================================*
c   SUN AND LATE-TYPE STAR ALPHA-OMEGA MEAN-FIELD INTERFACE 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 distributed within convection zone.
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
c    and for full CZ alpha-effect
c     parameter(rs2=0.95,rs3=0.98,ds2=0.01,ds3=0.01)
c -- these are for base-CZ alpha effect source term
      parameter(rs2=0.95,rs3=0.725,ds2=0.01,ds3=0.01)
      real tt1,tt2,tt3,g1r,g2r,g3r,thta,sinth,ftheta
      real erf
c---------------------------------------------------------------------------

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

         tt1=(yy-0.7)/(0.01)
         g1r=0.5*(1.+erf(tt1))
         tt3=(yy-rs3)/(ds3)
         g3r=0.5*(1.-erf(tt3))

c -- following for cos(theta) alpha-effect
c        ftheta=xx
c -- following for cos(theta) alpha-effect
c        thta=acos(xx)
c        sinth=sin(thta)
c        ftheta=2.598*xx*sinth**2
c -- following for equatorial alpha-effect
         thta=acos(xx)
         sinth=sin(4.*thta)
         if(thta.le.pi/4.)then
            ftheta=0.
         else
            ftheta=-sinth
         endif
         alpxy=1.08292*g1r*g3r *ftheta
      else
c --> exterior and deep interior: no source
	 alpxy=0.

      endif

      return
      end