subroutine circmd(xx,yy,ur,uth,durdr,duthdt)

c===========================================================================*
c   SUN AND LATE-TYPE STAR BABCOCK-LEIGHTON DYNAMOS:                  
c
c   computes components of meridional circulation speed and derivatives
c   at r/R=yy, cos(theta)=xx
c
c   The meridional circulation has one cell per quadrant, poleward
c   at surface r/R=1, penetrates slightly below the interface.
c   Mathematical expression given in Dikpati & Charbonneau 1999,
c   ApJ, 518, 508, eqs. (5a)--(5f), specialized to case qq=0.
c   radial flow derivative is computed by second-order finite
c   differences, theta-derivative analytically
c   The flow has a polytropic density profile built into it,
c   so no explicit density computation is needed.
c   Original flow derivation in van Ballegooijen & Choudhuri 1988
c
c   NOTE: Eqs. (5a) and (5b) in DC99 are missing a factor of 2
c         if u_0 is to be the max flow speed at surface mid-latitudes;
c         This is corrected here.
c
c   ur     is radial component of circulation
c   uth    is colatitudinal component of circulation
c   durdr  is radial derivative of radial component of circulation
c   durdth is colatitudinal derivative of colatitudinal component of circ.
c===========================================================================
c   Refs: Dikpati & Charbonneau 1999 , ApJ, 518, 508-ff
c         van Ballegooijen & Choudhuri 1988, ApJ, 333, 965
c
c   Last revised 10/07/2001
c===========================================================================
      implicit  none
      include  "arrays.global.f"
c Input
      real    xx,yy
c Output
      real    ur,uth,durdr,duthdt
c Local
      real pi
      parameter(pi=3.1415926536)
c -- these are for the meridional circulation pattern
      real rbot,mm,pp,qq,xib,c1,c2,cc0,cc1,cc2
c -- 5/7/01: rbot=0.7 has u_th peaking at rzc, =0.6 a bit below
c     parameter(rbot=0.7,mm=0.5,pp=0.25,qq=0.,xib=1./rbot-1.)
      parameter(rbot=0.61,mm=0.5,pp=0.25,qq=0.,xib=1./rbot-1.)
c     parameter( c1=(2.*mm+1.)*(mm+pp)/((mm+1.)*pp)*xib**(-mm) )
c     parameter( c2=(2.*mm+pp+1.)*mm/((mm+1.)*pp)*xib**(-mm-pp) )
c     parameter( cc0=1./(mm+1.),cc1=c1/(2.*mm+1.),cc2=c2/(2.*mm+pp+1) )
      real sin1,sin2,cos2,yyp,yym,xir,xirp,xirm,urp,urm,fthur,fruth
      integer ssign
c---------------------------------------------------------------------------
      ssign=1.
c     if(xx.le.0.) ssign=-1.
      if(yy.gt.1..or.yy.le.rbot)then
         ur =0.
	 uth=0.
	 durdr=0.
	 duthdt=0.

      else
         cos2=xx*xx
	 sin1=ssign*sin(acos(xx))
	 sin2=sin1*sin1
         yyp=yy+0.001
         yym=yy-0.001
         xir =1./yy -1.
         xirp=max(0.,1./yyp-1.)
         xirm=1./yym-1.

         c1=(2.*mm+1.)*(mm+pp)/((mm+1.)*pp)*xib**(-mm) 
         c2=(2.*mm+pp+1.)*mm/((mm+1.)*pp)*xib**(-mm-pp) 
         cc0=1./(mm+1.)
	 cc1=c1/(2.*mm+1.)
	 cc2=c2/(2.*mm+pp+1) 

         fthur=(2.*cos2-sin2)
         ur= 2.*(1./yy )**2 *( -cc0 +cc1*xir **mm -cc2*xir **(mm+pp) )
     +       *xir  *fthur
         urp=2.*(1./yyp)**2 *( -cc0 +cc1*xirp**mm -cc2*xirp**(mm+pp) )
     +       *xirp *fthur
         urm=2.*(1./yym)**2 *( -cc0 +cc1*xirm**mm -cc2*xirm**(mm+pp) )
     +       *xirm *fthur
         durdr =(urp -urm )/0.002

         fruth =(1./yy)**3 *( -1.+c1*xir**mm-c2*xir**(mm+pp) )
	 uth   =2.*fruth*sin1*xx
         duthdt=2.*fruth*(cos2-sin2)

      endif

      return
      end