subroutine srcpol(xx,yy,alpxy)

c===========================================================================*
c   SUN AND LATE-TYPE STAR ALPHA-OMEGA MEAN-FIELD 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 base-CZ alpha effect source term
      parameter(rs2=0.95,rs3=0.8,ds2=0.01,ds3=0.05)
      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: Classical mean-field alpha-effect source
c     (and avoiding IEEE underflows)

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

c -- following for cos(theta) alpha-effect
         ftheta=xx
c -- following for equatorial alpha-effect
c        thta=acos(xx)
c        sinth=sin(4.*thta)
c        ftheta=-sin(4.*thta)
c        if(thta.le.pi/4.) ftheta=0. 
c -- following for sin^2(theta)cos(theta) alpha-effect
c        thta=acos(xx)
c        sinth=sin(thta)
c        ftheta=2.598*xx*sinth**2

c        alpxy=1.46175*g1r*g3r *ftheta
         alpxy=1.178*g1r*g3r *ftheta
c        alpxy=1.*g1r*g3r *ftheta
      else
c --> exterior and deep interior: no source
	 alpxy=0.

      endif

      return
      end