      program bldynamo
c**********************************************************************
c                                                                     *
c                       Babcock-Leighton dynamos                      *
c                       ========================                      *
c                                                                     *
c Solves 2D Axisymmetric MHD induction equation in a meridional       *
c plane or quadrant, as nonlinear initial-boundary-value problem      *
c using galerking finite element method and one-step implicit         *
c time-stepping (theta-method)                                        *
c                                                                     *
c Nonlinearity is non-local Babcock-Leighton poloidal source term     *
c Uses solar-like differential rotation, magnetic diffusivity         *
c contrast across core-envelope interface, solar-like single-cell     *
c meridional circulation in convective envelope                       *
c                                                                     *
c Codes solves for toroidal component (B_phi) and toroidal vector     *
c potential (A_phi such that B_p=curl A)                              *
c                                                                     *
c This version has buoyant loss term, as per dario's model            *
c                                                                     *
c Lengths measured in units of solar radius R                         *
c Time measured in units of diffusion time tau=R^2/eta_T, where       *
c   eta_T is envelope magnetic diffusivity                            *
c  mu=cos(theta) used as latitudinal variable.                        *
c=====================================================================*
c References:
c
c Burnett, D.S. 1987, Finite Element Analysis, Addison-Wesley [Burnett]
c
c Dikpati, M., & Charbonneau, P. 1999, ApJ, 518, 508          [DC99]
c
c Golub, & Van Loan, Matrix Computations, 2nd ed., Johns Hopkins
c
c=====================================================================*
c     Version 1.0,    11/7/2001      Paul Charbonneau, HAO/NCAR
c     Last modified   27/11/2002     Paul Charbonneau, Physique/UdeM
c**********************************************************************
      implicit none

      include "size.f"
      include "arrays.global.f"
      include "arrays.mesh.f"
c Local
      integer        it
      character*7    direct
      character*64   infile,outfile,rstfile
      integer        seed
      parameter      (seed=54321)
      real           urand
      external       urand

c---------- I/O control parameters and files
c
c     All output files go to directory [direct]
c
c     all    : output  file [outfile] solution at each time step
c     incr=1 : restart file [infile]  being used as initial condition
c     incw=1 : restart file [rstfile] being written
c
      incr=1
      incw=1
c     direct= '/mephisto/f/paulchar/'
      direct= 'saasfe/'
      infile =direct //'BL_ivp.3030t2b.rst'
      outfile=direct //'BL_ivp.3030t2h.i3e'
      rstfile=direct //'BL_ivp.3030t2h.rst'

      open(8,file=outfile,form='unformatted')
      if(incr.eq.1) open(7,file=infile, form='unformatted')
      if(incw.eq.1) open(9,file=rstfile,form='unformatted')
c---------------------------------------------------------------------------

c---------- Model subroutines
c
c     The code requires various subroutines defining the model's
c     input functionals, such as the meridional flow, differential
c     rotation, etc. These are incorporated via "include" statements
c     at the end of this source code.
c
c     mod.angvel.f : angular velocity as a function of mu and r
c     mod.magdif.f : magnetic diffusivity as a function of r
c     mod.circmd.f : meridional flow components as a function of mu and r
c     mod.srcpol.f : source term as a function of mu and r
c
c---------- Global control parameters 
c
c  1. ipar   : Controls imposed solution parity, or lack thereof
c
c     ipar=0 : full meridional plane solution
c              theta-domain is from -pi to pi, mixed parity initial condition
c     ipar=1 : N-meridional quadrant solution
c              theta-domain is from 0 to pi, dipole-like initial condition
c     ipar=2 : N-meridional quadrant solution
c              theta-domain is from 0 to pi, quadrupole-like initial condition
c
c  2. ialp   : Controls inclusion of mean-field alpha-effect
c
c     ialp=0 : no mean-field alpha-effect
c     ialp=1 : alpha-effect only in poloidal equations (alpha-omega)
c     ialp=2 : alpha-effect in both equations (alpha^2-omega)
c     
c  3. istc   : Controls inclusion of stochastic effects [UNDER CONSTRUCTION]
c
c     istc=0 : steady model
c     istc=1 : fluctuations in source term only            
c     istc=2 : fluctuations in meridional circulation only [NOT THERE YET]
c     istc=3 : fluctuations in source and circulation      [NOT THERE YET]
c
      ipar=1
      ialp=1
      istc=0
c---------------------------------------------------------------------------

c---------- define problem data (see also DC99, sec. 2.4)
c
c    rzc      : core-envelope interface; for sun,      =0.7 
c    ybl      : tachocline half-thickness; for sun,    =0.025
c    calpha   : alpha/source dynamo number:            =s_0 R/eta_T
c    comega   : Differential Rotation dynamo number    =Omega_0 R^2/eta_T
c    reynolds : Reynolds number based on circulation   =u_0 R/eta_T
c    etadf    : envelope-to-core diffusivity ratio     =eta_T/eta_c (>=1)
c    aaa      : Lorentz force parameter in circulation equation
c    bbb      : Decay parameter in circulation equation
c
      ybl=0.025
      rzc=0.7
c -- following for ref model of DC99 (Sec. 3)
      title='B-L Ref.:128x96 Ca=25 deta=300 Rm=840expl B '
c following for BL Ref
      calpha  = 25.
      comega  = 5.e4
      etadf   = 300.
      reynolds= 840.
      reynld2=  0.
      gamfl=0.
      aaa=3.e-1
      bbb=100.
c     aaa=0.
c     bbb=0.

c------ this is the amplitude of the additive noise to vector potential A
c       at the surface
c     addnamp=1.0e+03
      addnamp=0.

c---------------------------------------------------------------------------

c---------- following spatial mesh and time-stepping parameters are set in
c           include file size.f
c     nelx   : number of elements in latitudinal direction
c     nely   : total number of elements in radial direction
c     ityp   : element type (23=bilinear, 24=8-node quadratic)
c     ngp    : number of gaussian quadrature points per dimension
c                (=2 for bilinear elements, =3 for 8-node quadratic)
c     nsteps : number of sub-time steps (=1 for nonlinear problems)
c     theta  : time-stepping parameter (0.5<theta<1.0; 2/3 works well)
c
c     except the following:
c
c     nt:     number of main time steps
c     to:     duration of simulation, in diffusion times
c
      to=0.1
      nt=1

c---------- output some header info
      call header(infile,outfile,rstfile)

c---------- initialize random number generator
      call urand_init(seed)

c---------- mesh generation
      call mesh2d

c---------- some safety checks on flow fields, etc
c     call sfcheck
c     stop ' *TMP*' 

c---------- initial condition 
      call initcnd

c---------- compute energies in initial state; write initial condition
      ijac=1
      call postpro(ijac)
      call writout(ijac-1,0.)

c********** main time-stepping loop begins
      do it=1,nt
         ijac=it+1
c---------- assemble system matrices
         call formm
c        call formf
c---------- apply boundary conditions
         call aplybc
c---------- solver
         call calfun(it)
c---------- compute magnetic energy and such
c        call postpro(it)
c---------- output nodal time series
c        call writout(it)
c---------- additive noise
c        call addnoise
      enddo
c********** main time-stepping loop ends

c---------- produce restart file from final time-step
      if(incw.eq.1) call stlast

c---------- final file write status report
      write(*,'(a64)') title
      if(incr.eq.1) write(*,'(a64)') infile
      write(*,'(a64)') outfile
      write(*,'(a64)') rstfile

      stop ' **END**'
      end

c**************************************************************************
      subroutine header(infile,outfile,rstfile)

c=========================================================================
c  Outputs to screen some basic header info, including run parameters
c=========================================================================
      implicit none
      include "size.f"
      include "arrays.global.f"
      include "arrays.mesh.f"
c Input
      character*64 infile,outfile,rstfile
      character*8 labclock
      character*9 labdate
c-------------------------------------------------------------------------
      call date(labdate)
      call time(labclock)
      write(*,1) title,labdate,labclock
      if(incr.eq.0) write(*,2)
      if(incr.eq.1) write(*,3) infile
      write(*,4) outfile
      write(*,5) rstfile
      write(*,6) calpha,comega,reynolds,etadf,ipar,rzc,ybl
      write(*,7) nelx,nely,ityp,ngp,nt,nsteps,to,theta

      return
c.........................................................................
    1 format(1x,60('*'),/,1x,'*',10x,
     +'NON-LINEAR KINEMATIC BABCOCK-LEIGHTON DYNAMOS',/,
     +1x,60('*'),/,1x,'Run :',a80,/,1x,'Date: ',a9,1x,a8)
    2 format(1x,'Starting from t=0 initial condition')
    3 format(1x,'Extending an earlier run. Restart file: ',/,5x,a80)
    4 format(1x,'Output  file: ',a64)
    5 format(1x,'Restart file: ',a64)
    6 format(/,1x,'Input parameters:',/,
     +         5x,'C_alpha      = ',e12.4,/,
     +         5x,'C_omega      = ',e12.4,/,
     +         5x,'Reynolds     = ',e12.4,/,
     +         5x,'Delta eta    = ',e12.4,/,
     +         5x,'Parity       = ',i3,/,
     +         5x,'R_zc/R       = ',e12.4,/,
     +         5x,'w/R          = ',e12.4)
    7 format(/,1x,'Mesh parameters:',/,
     +         5x,'nelx         = ',i4,/,
     +         5x,'nely         = ',i4,/,
     +         5x,'element type = ',i4,/,
     +         5x,'N_GP         = ',i4,/,
     +         5x,'nt           = ',i5,/,
     +         5x,'nsteps       = ',i4,/,
     +         5x,'to           = ',e12.4,/,
     +         5x,'theta        = ',e12.4)
      end

c***************************************************************************
      subroutine elem23c(ielno,iflag)

c===========================================================================
c  master routine for element 23c, a quadrilateral bilinear isoparametric    
c  element for a system of two coupled parabolic PDE                        
c
c  ielno  : numbering of element for which computation is carried out
c  iflag=1: computes element matrices   (call from subroutine form)
c  iflag=2: magnetic energy calculation (call from subroutine postpro)
c
c  Called by: subroutines form, post
c  Calls    : subroutines shap23, angvel, magdif, circmd, srcpol
c===========================================================================
      implicit none
      include  "size.f"
      include  "arrays.global.f"
      include  "arrays.mesh.f"
      include  "arrays.femlocal.f"
      include  "arrays.femglobal.f"
c Input
      integer  ielno,iflag
c Local 
      integer  i,j,jj,intx,inty,igp,i1,i2,j1,j2
      real     ybot,tdkill,wgtjac,rr,mu,r2,sin1,sin2,
     +         omg,domgdr,domgdm,remag,drmfac,
     +         ur,uth,durdr,duthdth,
     +         ckat,ckbt,ckarr,ckamm,ckar,cka,
     +         ckbmm,ckbrr,ckbr,ckb,ckb1,ckb3,ckb5,
     +         cadv1,cadv2,cadv3,cadv4,cadv5,cadv6
c---------------------------------------------------------------------------
c********** form elements matrices and vectors
c---------- initialisations
      do i=1,4
         jj=icon(ielno,i)
         xx(i)=x(jj)
         yy(i)=y(jj)
	 do j=1,8
            ke(j,2*i)  =0.
            ce(j,2*i)  =0.
            ke(j,2*i-1)=0.
            ce(j,2*i-1)=0.
         enddo
      enddo
      ybot=y(icon(ielno,1))

c---------- take care of exterior solution
      if(ybot.ge.1.)then
c --> within buffer;
	 tdkill=0.
c        write(*,*) ielno,ybot
      else
         tdkill=1.
      endif

c---------- integrate terms (gaussian quadrature of order ngp X ngp)
      igp=0
      do intx=1,ngp
         xi=gp(ngp,intx)
         do inty=1,ngp
            eta=gp(ngp,inty)
            igp=igp+1
            call shap23(ielno,igp)
            wgtjac=wt(ngp,intx)*wt(ngp,inty)*jac(ielno,igp)
	   
c -- some geometric terms
	    mu  =xgp(ielno,igp)
	    rr  =ygp(ielno,igp)
	    sin2=1.-mu*mu
	    sin1=sin(acos(mu))
	    r2  =rr*rr

c -- compute angular velocity and derivatives at (r,mu)
            call angvel(mu,rr,omg,domgdr,domgdm)

c -- compute magnetic diffusivity and radial derivative at (r,mu)
            call magdif(mu,rr,remag,drmfac)

c -- compute meridional circulation and derivatives at (r,mu)
            call circmd(mu,rr,ur,uth,durdr,duthdth)

c********** THESE ARE THE EQUATIONS BEING SOLVED ************************
c
c    Poloidal field is written as B_pol= curl (0,0,A)
c    Latitudinal dependent variable is mu=cos(theta)
c
c -- Poloidal Equation: (see Tassoul 1978, eq. (15.20); here times r^2)
c
c        dA    d          dA     d        dA    
c   ckat -- = ---( ckamm ---) + --( ckarr --) + cka A   [DIFFUSION]
c        dt   dmu        dmu    dr        dr   
c
c             dA
c      + ckar --                                        [DIAMAGNETISM]
c             dr
c
c              dA                    dA
c      + cadv1 -- + cadv2 A + cadv3 ---                 [ADVECTION BY CM]
c              dr                   dmu
c
c      + ckas                                           [SOURCE TERM]
c
c -- Toroidal Equation: (see Tassoul 1978, eq. (15.18); times r^3 sin theta)
c
c        dB    d          dB     d        dB   
c   ckbt -- = ---( ckbmm ---) + --( ckbrr --) + ckb B   [DIFFUSION]
c        dt   dmu        dmu    dr        dr  
c
c             dB
c      + ckbr --                                        [DIAMAGNETISM]
c             dr
c
c              dB                    dB
c      + cadv4 -- + cadv5 B + cadv6 ---                 [ADVECTION BY CM]
c              dr                   dmu
c
c             dA                  dA
c      + ckb3 -- + ckb5 A + ckb1 ---                    [SHEARING BY DR]
c             dr                 dmu
c
c*************************************************************************

c -- time derivative terms
            ckat= r2*tdkill
            ckbt= r2

c -- diffusion terms, for both equations
c     REDOUBLEDCHECKED 26/8/2002 AND FOUND CORRECT, EXCEPT FOR ckbr, CORRECTED
	    ckarr =  r2      *remag
	    ckar  = -r2      *drmfac
	    ckamm =  sin2    *remag
            cka   = -1./sin2 *remag
	    ckbrr =  ckarr
	    ckbmm =  ckamm
c     NEXT TERM FOUND MISSING 26/8/2002 AND ADDED (need to double check)
c           ckbr  =  rr      *drmfac
            ckb   =  cka

c -- Shearing terms in toroidal field (B) equation
c     REDOUBLEDCHECKED    4/1/2000 AND FOUND CORRECT
c     RE-REDOUBLEDCHECKED 26/8/2002 AND FOUND CORRECT
            ckb1=-r2*sin2*domgdr                   *comega
            ckb3= r2*sin2*domgdm                   *comega
            ckb5= (r2*mu*domgdr+rr*sin2*domgdm)    *comega

c --> terms multiplying A in toroidal equation (alpha^2 contribution;
c     currently commented out in matrix element calculation below)
c     REDOUBLEDCHECKED 4/1/2000 AND FOUND CORRECT
c           ckba1=-alp0*r2                         *calpha*term*ialp
c           ckba2=-alp0*sin2                       *calpha*term*ialp
c           ckab3=(alp0/sin2-rr*dalpdr+dalpdmu*mu )*calpha*term*ialp

c -- advection terms due to meridional circulation, both equations
c     CORRECTED         4/1/2000 for missing r^2 factors
c     REDOUBLEDCHECKED 26/8/2002 AND FOUND CORRECT
            cadv1=- r2*ur                          *reynolds
	    cadv2=-(ur + uth*mu/sin1)*rr           *reynolds
	    cadv3=  uth*sin1*rr                    *reynolds
	    cadv4=- r2*ur                          *reynolds
	    cadv5=-(ur*rr+r2*durdr+rr*duthdth)     *reynolds
	    cadv6=  uth*rr*sin1                    *reynolds
c           cadv1=0.
c           cadv2=0.
c           cadv3=0.
c           cadv4=0.
c           cadv5=0.
c           cadv6=0.

            do i=1,4
	       i1=2*i-1
	       i2=2*i
	       do j=1,4
	          j1=2*j-1
	          j2=2*j

c --> coefficients of poloidal induction equation
c
                  ce(i1,j1)=ce(i1,j1)+wgtjac*ckat*   phi(j)   *phi(i)
                  ke(i1,j1)=ke(i1,j1)+wgtjac*(
     +                                   ckamm    *dphdx(j) *dphdx(i)
     +                                  +ckarr    *dphdy(j) *dphdy(i)
     +                                  -ckar     *dphdy(j)   *phi(i)
     +                                  -cka        *phi(j)   *phi(i)
     +                                  -cadv1    *dphdy(j)   *phi(i)
     +                                  -cadv2      *phi(j)   *phi(i)
     +                                  -cadv3    *dphdx(j)   *phi(i))
c                 ke(i1,j2)=ke(i1,j2)+wgtjac*(
c    +                                  -cka1       *phi(j)   *phi(i))
c --> coefficients of toroidal induction equation
c
                  ce(i2,j2)=ce(i2,j2)+wgtjac*ckbt*   phi(j)   *phi(i)
                  ke(i2,j1)=ke(i2,j1)+wgtjac*(
c    +                                   ckba2    *dphdx(j) *dphdx(i)
c    +                                  +ckba1    *dphdy(j) *dphdy(i)
     +                                  -ckb1     *dphdx(j)   *phi(i)
     +                                  -ckb3     *dphdy(j)   *phi(i)
     +                                  -ckb5       *phi(j)   *phi(i))
c    +                                  -ckab3      *phi(j)   *phi(i))
                  ke(i2,j2)=ke(i2,j2)+wgtjac*(
     +                                   ckbmm    *dphdx(j) *dphdx(i)
     +                                  +ckbrr    *dphdy(j) *dphdy(i)
     +                                  -ckbr     *dphdy(j)   *phi(i)
     +                                  -ckb        *phi(j)   *phi(i)
     +                                  -cadv4    *dphdy(j)   *phi(i)
     +                                  -cadv5      *phi(j)   *phi(i)
     +                                  -cadv6    *dphdx(j)   *phi(i))
               enddo
            enddo
         enddo
      enddo

      return
c-----------------------------------
      end

c***************************************************************************
      subroutine elem23f(ielno,iflag)

c===========================================================================
c  master routine for element 23c, a quadrilateral bilinear isoparametric    
c  element for a system of two coupled parabolic PDE                        
c
c  ielno  : numbering of element for which computation is carried out
c  iflag=1: computes element matrices   (call from subroutine form)
c  iflag=2: magnetic energy calculation (call from subroutine postpro)
c
c  Called by: subroutines form, post
c  Calls    : subroutines shap23, angvel, magdif, circmd, srcpol
c===========================================================================
      implicit none
      include  "size.f"
      include  "arrays.global.f"
      include  "arrays.mesh.f"
      include  "arrays.femlocal.f"
      include  "arrays.femglobal.f"
c Input
      integer  ielno,iflag
c Local 
      integer  i,j,ii,jj,intx,inty,igp,i1,i2,ielno2,ii1,ii2,i3,i4
      real     sum2,term,ybot,wgtjac,rr,mu,r2,sin1,sin2,xbot,
     +         alp0,ckas,ckaf,term2,remag,drmfac,
c    +         sum1,sum3,bgpt,agpt,axgpt,aygpt,br2,bth2,bgpt2,ebt,ebtt,
     +         sum1,sum3,bgpt,agpt,axgpt,aygpt,br2,bth2,bgpt2,
     +         bavrg,alpsfac0,dx,dy,
     +         ur,uth,durdr,duthdth,
     +         cadv1,cadv2,cadv3,cadv4,cadv5,cadv6,
     +         agp,bgp,dadxgp,dadygp,dbdxgp,dbdygp
c---------------------------------------------------------------------------
c********** form elements matrices and vectors
      if(iflag.eq.1)then
c---------- initialisations
      do i=1,4
         jj=icon(ielno,i)
         xx(i)=x(jj)
         yy(i)=y(jj)
         fe(2*i)=0.
         fe(2*i-1)=0.
      enddo
      ybot=y(icon(ielno,1))
      xbot=x(icon(ielno,1))
c     write(*,*) ' ielno,xbot ', ielno,xbot
      alpsfac0=alpsfac
      if(xbot.ge.0.) alpsfac0=alpsfac2

c---------- compute A,B and derivatives in middle of element
c     agp=0.
c     bgp=0.
c     dx=x(icon(ielno,2))-x(icon(ielno,1))
c     dy=y(icon(ielno,4))-y(icon(ielno,1))
c     do i=1,4 
c        i1=2*icon(ielno,i)-1
c        i2=i1+1
c        agp   =agp   +a(i1)/4.
c        bgp   =bgp   +a(i2)/4.
c     enddo
c     i1=2*icon(ielno,1)-1
c     i2=2*icon(ielno,2)-1
c     i3=2*icon(ielno,3)-1
c     i4=2*icon(ielno,4)-1
c     dadxgp=(a(i2)+a(i3)-a(i4)-a(i1))/(2.*dx)
c     dadygp=(a(i4)+a(i3)-a(i2)-a(i1))/(2.*dy)
c     i1=i1+1
c     i2=i2+1
c     i3=i3+1
c     i4=i4+1
c     dbdxgp=(a(i2)+a(i3)-a(i4)-a(i1))/(2.*dx)
c     dbdygp=(a(i4)+a(i3)-a(i2)-a(i1))/(2.*dy)
c---------- compute nonlinear, nonlocal part of source term (DC99, eq. 7)
c           (constant within element!!) Surface source strength is proportional
c           to toroidal field strength at core-envelope inteface at same
c           latitude, with alpha-quenching-like nonlinearity in terms
c           of interface toroidal field:
c
c           B_CZ/( 1.+( B_CZ/B_0 )^2 ) , with B_0=10^5 G
c
      if(ybot.gt.0.5.and.ybot.le.1.0)then
         sum2=0.
c --> first compute B_phi in middle of "cousin" element located
c     beneath core-envelope interface (element numbering ielno2)
	 ielno2=ielnoint+mod(ielno,nelx)
         do i=1,4
c -- following for local alpha-effect
c           jj=2*icon(ielno,i)
c -- following for non-local alpha-effect
            jj=2*icon(ielno2,i)
            sum2=sum2+a(jj)
         enddo
         sum2=sum2/4.
c --> this, finally, is nonlinear nonlocal part of source term,
c     including quenching.
c      following quenching for alpha-Omega
          term=sum2/( 1.+(sum2/5.e3)**2 )
c          term=erf(sum2)*sum2/( 1.+(sum2/3.1)**2 )**2
c      following quenching for BL
c          term=(1.+erf( (abs(sum2)-6.e4)/2.e4 ))
c    +         *(1.-erf( (abs(sum2)-1.e5)/8.e4 ))*sum2
c        if(ielno.eq.2776)then
c           write(*,*) ielno,sum2,term
c        endif
         term2=sum2*sum2*sum2
      else
c --> away from source region, so term=0.
         term=0.
      endif

c---------- integrate terms (gaussian quadrature of order ngp X ngp)
      igp=0
      do intx=1,ngp
         xi=gp(ngp,intx)
         do inty=1,ngp
            eta=gp(ngp,inty)
            igp=igp+1
            call shap23(ielno,igp)
            wgtjac=wt(ngp,intx)*wt(ngp,inty)*jac(ielno,igp)

c -- compute A, B, and derivatives, at gauss points
            agp=0.
            bgp=0.
            dadxgp=0.
            dadygp=0.
            dbdxgp=0.
            dbdygp=0.
            do i=1,4 
               i1=2*icon(ielno,i)-1
               i2=i1+1
               agp   =agp   +a(i1)*phi(i)
               dadxgp=dadxgp+a(i1)*dphdx(i)
               dadygp=dadygp+a(i1)*dphdy(i)
               bgp   =bgp   +a(i2)*phi(i)
               dbdxgp=dbdxgp+a(i2)*dphdx(i)
               dbdygp=dbdygp+a(i2)*dphdy(i)
            enddo
c -- some geometric terms
	    mu  =xgp(ielno,igp)
	    rr  =ygp(ielno,igp)
	    r2  =rr*rr
	    sin1=sin(acos(mu))
c -- compute spatial part of poloidal source term at (r,mu)
            call srcpol(mu,rr,alp0)

c -- compute buoyant loss coefficient (using magnetic diffusivity
c    profile for now...)
c           call magdif(mu,rr,remag,drmfac)

c --> Babcock-Leighton poloidal source term 
	    ckas= r2*alp0                 *alpsfac0*calpha *term

c --> Flux loss term (zero in basic model)
c           ckaf= -r2*gamfl*remag *term2
	    ckaf= 0.

c -- advection terms due to meridional circulation, both equations
c     CORRECTED         4/1/2000 for missing r^2 factors
c     REDOUBLEDCHECKED 26/8/2002 AND FOUND CORRECT

c -- compute meridional circulation and derivatives at (r,mu)
            call circmd(mu,rr,ur,uth,durdr,duthdth)

            cadv1=- r2*ur                          *reynld2
	    cadv2=-(ur + uth*mu/sin1)*rr           *reynld2
	    cadv3=  uth*sin1*rr                    *reynld2
	    cadv4=- r2*ur                          *reynld2
	    cadv5=-(ur*rr+r2*durdr+rr*duthdth)     *reynld2
	    cadv6=  uth*rr*sin1                    *reynld2
c           cadv1=0.
c           cadv2=0.
c           cadv3=0.
c           cadv4=0.
c           cadv5=0.
c           cadv6=0.
            do i=1,4
	       i1=2*i-1
	       i2=2*i
c              ii1=2*icon(ielno,i)-1
c              ii2=ii1+1
c              agp=a(ii1)*phi(i)
c              bgp=a(ii2)*phi(i)
c              dadxgp=a(ii1)*dphdx(i)
c              dbdxgp=a(ii2)*dphdx(i)
c              dadygp=a(ii1)*dphdy(i)
c              dbdygp=a(ii2)*dphdy(i)
c --> source term on RHS of poloidal induction equation
               fe(i1)=fe(i1)+wgtjac*(   ckas*phi(i)
c --> alpha-source term on RHS of toroidal induction equation
c              fe(i1)=fe(i1)+wgtjac*ckas*phi(i)
c --> meridional circulation transport term on RHS of poloidal equation
     +                                  +cadv1*dadygp *phi(i)
     +                                  +cadv2*agp    *phi(i)
     +                                  +cadv3*dadxgp *phi(i))
c --> flux loss term on RHS of toroidal induction equation
               fe(i2)=fe(i2)+wgtjac*(    ckaf*phi(i)
c --> meridional circulation transport term on RHS of toroidal equation
     +                                  +cadv4*dbdygp *phi(i)
     +                                  +cadv5*bgp    *phi(i)
     +                                  +cadv6*dbdxgp *phi(i))
            enddo

         enddo
      enddo

      return

c********** do magnetic energy calculations
      else

c---------- initialisations
      do i=1,4
         j=icon(ielno,i)
         xx(i)=x(j)
         yy(i)=y(j)
      enddo
c---------- integrate terms (gaussian quadrature of order ngp X ngp)
      igp  =0
      sum1=0.
      sum2=0.
      sum3=0.
      do intx=1,ngp
         xi=gp(ngp,intx)
         do inty=1,ngp
            eta=gp(ngp,inty)
            igp=igp+1
            call shap23(ielno,igp)
            wgtjac=wt(ngp,intx)*wt(ngp,inty)*jac(ielno,igp)

	    mu  =xgp(ielno,igp)
	    rr  =ygp(ielno,igp)
	    sin2=1.-mu*mu
	    sin1=sin(acos(mu))
	    r2  =rr*rr

	    bgpt =0.
	    agpt =0.
	    axgpt=0.
	    aygpt=0.
            do i=1,4
	       ii=2*icon(ielno,i)-1
	       agpt =agpt +a(ii  )*phi(i)
	       bgpt =bgpt +a(ii+1)*phi(i)
	       axgpt=axgpt+a(ii)  *dphdx(i)
	       aygpt=aygpt+a(ii)  *dphdy(i)
            enddo

c --> compute field magnitude and energies
c     the forms of B_r and B_theta below suggest B=curl(A), not curl(A/r sin)
	    br2 =(-axgpt*sin1/rr+agpt*mu/rr/sin1)**2
	    bth2=(-aygpt-agpt/rr)**2
	    bgpt2=bgpt**2
	    sum1=sum1+(bgpt2+br2+bth2)     *wgtjac
	    sum2=sum2+bgpt2                *wgtjac
	    sum3=sum3+sqrt(bgpt2+br2+bth2) *wgtjac
         enddo
      enddo

c --> this is total magnetic energy
      ebt=sum1
c --> this is toroidal magnetic energy
      ebtt=sum2
c --> this is average magnetic field strength
      bavrg=sum3

      endif

      return
c-----------------------------------
      end

c************************************************************************
      subroutine shap23(ielno,igp)

c========================================================================
c  Internal to FEM.
c  Evaluates shape function and derivatives at points (xi,eta)             
c  for 4-node quadrilateral. Jacobian is only computed upon first
c  call (controled by pareameter ijac) and stored for future use
c
c  Called by: subroutine elem23c
c========================================================================
      implicit none
      include "size.f"
      include "arrays.global.f"
      include "arrays.femlocal.f"
c Input
      integer ielno,igp
c Local
      integer i
      real    xip,xim,etp,etm,dxdxi,dydet,dydxi,dxdet,avg
c---------- constants
      xip=1.+xi
      xim=1.-xi
      etp=1.+eta
      etm=1.-eta
c******************** shape functions
      phi(1)=xim*etm/4.
      phi(2)=xip*etm/4.
      phi(3)=xip*etp/4.
      phi(4)=xim*etp/4.
c******************** derivatives of the shape function
c---------- derivatives wrt xi at (xi,eta)
      dphdxi(1)=-etm/4.
      dphdxi(2)= etm/4.
      dphdxi(3)= etp/4.
      dphdxi(4)=-etp/4.
c---------- derivatives wrt eta at (xi,eta)
      dphdet(1)=-xim/4.
      dphdet(2)=-xip/4.
      dphdet(3)= xip/4.
      dphdet(4)= xim/4.
c******************** do jacobian calculations
c---------- compute jacobian at (xi,eta)
      if(ijac.le.1)then
         dxdxi=0.
         dydet=0.
         dydxi=0.
         dxdet=0.
         do i=1,4
            dxdxi=dxdxi+dphdxi(i)*xx(i)
            dydet=dydet+dphdet(i)*yy(i)
            dydxi=dydxi+dphdxi(i)*yy(i)
            dxdet=dxdet+dphdet(i)*xx(i)
c    write(*,*) ielno,xx(i),yy(i)
         enddo
         jac(ielno,igp)=dxdxi*dydet-dxdet*dydxi
c---------- check positiveness of jacobian
         avg=(abs(dxdxi)+abs(dydet))/4.
         if(jac(ielno,igp).lt.avg*1.e-8)then
            write(*,*) ielno,jac(ielno,igp)
            stop '1 shap23'
         end if
c---------- compute inverse of the jacobian at (xi,theta)
         jacinv(ielno,igp,1)=dydet/jac(ielno,igp)
         jacinv(ielno,igp,2)=dxdxi/jac(ielno,igp)
         jacinv(ielno,igp,3)=-dydxi/jac(ielno,igp)
         jacinv(ielno,igp,4)=-dxdet/jac(ielno,igp)
      end if
c---------- derivatives wrt x and y at (xi,eta)
      do i=1,4
         dphdx(i)=jacinv(ielno,igp,1)*dphdxi(i)
     +           +jacinv(ielno,igp,3)*dphdet(i)
         dphdy(i)=jacinv(ielno,igp,2)*dphdet(i)
     +           +jacinv(ielno,igp,4)*dphdxi(i)
      enddo
      return
c----------------------------------------
      end

c******************************************************************************
      subroutine grid1
 
c===========================================================================
c  genere une grid d'elements quadrilateraux bilineaires; calcule
c  coordonnees des noeuds, et matrtice de connectivite
c
c  Called by: subroutine mesh2d                                  
c===========================================================================
      implicit none
      include  "size.f"
      include  "arrays.mesh.f"
c Local
      integer nx,i,j,iel,i2

      nx=nelx+1
c     write(*,*) ' grid: entering' 
c---------- 1. calcul du "fond"
 
      do i=1,nx
         y(i)=yrad(1)
         x(i)=xtht(i)
c        write(*,*) i,x(i),y(i)
      enddo
c---------- 2. nely couches d'elements
c              (nelx elements quadrilateraux bilineaires par couche)
 
      do j=1,nely
         i2=1
         do i=(j-1)*nx+nx+1,j*nx+nx
            y(i)=yrad(j+1)
            x(i)=xtht(i2)
            i2=i2+1
c           write(*,*) i,x(i),y(i)
         enddo
      enddo
      do i=1,nely
         do j=1,nelx
            iel=(i-1)*nelx+j
            icon(iel,1)=(i-1)*nx+j
            icon(iel,2)=(i-1)*nx+j+1
            icon(iel,3)=i*nx+j+1
            icon(iel,4)=i*nx+j
c           ntype(iel)=ityp
c           ngpint(iel)=ngp
         enddo
      enddo

c     write(*,1) nelx,nely,nel,nnd,nm
    1 format(/,1x,'***** GRID:',5x,i3,' X ',i3,' Bilinear Elements',/,
     +         7x,'Total Elements:     ',i6,/,
     +         7x,'Total Nodes:        ',i6,/,
     +         7x,'Total Coefficients: ',i6)

      return
      end

c***************************************************************************
      subroutine gaussp
 
c===========================================================================
c  calcule la position des points de gauss, ainsi que les coefficients de    
c  l'equation differentielle en ces points
c
c  Called by: subroutine mesh2d                                  
c===========================================================================
      implicit none
      include  "size.f"
      include  "arrays.mesh.f"
      include  "arrays.femlocal.f"
c Local
      real     sqr13,sqr35,xl,yl,xmid,ymid
      integer  i

      sqr13=sqrt(1./3.)
      sqr35=sqrt(3./5.)
      do i=1,nel
         ymid=(y(icon(i,4))+y(icon(i,1)))/2.
         xmid=(x(icon(i,2))+x(icon(i,1)))/2.
         xl=(x(icon(i,2))-x(icon(i,1)))/2.
         yl=(y(icon(i,4))-y(icon(i,1)))/2.
c  1 point de gauss
         if(ngp.eq.1)then
            xgp(i,1)=xmid
            ygp(i,1)=ymid
         end if
c  2x2 points de gauss
         if(ngp.eq.2)then
            xgp(i,1)=xmid-xl*sqr13
            xgp(i,3)=xmid+xl*sqr13
            xgp(i,2)=xmid-xl*sqr13
            xgp(i,4)=xmid+xl*sqr13
            ygp(i,1)=ymid-yl*sqr13
            ygp(i,3)=ymid-yl*sqr13
            ygp(i,2)=ymid+yl*sqr13
            ygp(i,4)=ymid+yl*sqr13
         end if
      enddo

      return
      end

c***************************************************************************
      subroutine mesh2d

c===========================================================================
c  Generation du maillage spatial 2-D, maillage temporel, et calcul des    
c                          conditions limites                               
c
c  Called by: program bldynamo
c  Calls    : subroutines  grid1, gaussp
c===========================================================================
      implicit none
      include  "size.f"
      include  "arrays.global.f"
      include  "arrays.mesh.f"
      include  "arrays.femlocal.f"
      include  "arrays.femglobal.f"
c Local
      real     xxi,xxo,yi,yo,tti,tto,pi
      integer  i,nel1,ny1,nel2,ny2,nelbuf,nx,ii,i1,i2
      data     pi/3.1415926536/
c---------------------------------------------------------------------------

c========== set up radial mesh
c---------- radial direction (variable y, --> r/R )
c
c     nely/3  elements from base to C-E interface, concentrated towards latter;
c     nely/10 elements in outer buffer shell, concentrated towards r/R=1;
c     remainder of elements in convective envelope, equally spaced
      nel1=nely/3
      nelbuf=nely/10
      write(*,*) ' nel1,nelbuf', nel1,nelbuf
c --> in radiative core
      yi=0.5
      yo=rzc
      ny1=nel1+1
      do i=1,ny1
         yrad(i)=(yo-yi)*((float(i-1)/(ny1-1))**0.75)+yi
c        write(*,*) i,yrad(i)
      enddo
c     write(*,'()')
      ielnoint=nelx*(nel1-1)
c --> in convective envelope
      yi=rzc
      yo=1.
      nel2=nely-nel1-nelbuf
      ny2=nel2+1
      do i=1,ny2
         yrad(ny1+i)=(yo-yi)*((float(i)/(ny2))**1.)+yi
c        write(*,*) ny1+i,yrad(ny1+i)
      enddo 
c     write(*,'()')
c --> from surface to top of buffer zone
      yi=yo
      yo=3.0
      nelbuf=nelbuf-1
      do i=1,nelbuf
         yrad(ny1+ny2+i)=(yo-yi)*((float(i)/(nelbuf))**1.75)+yi
c        write(*,*) ny1+ny2+i,yrad(ny1+ny2+i)
      enddo
c---------- angular direction (variable x, --> mu=cos(theta) )
c           Mesh is in constant theta-increments
      if(ipar.ne.0)then
c --> solving only in N-hemisphere
	 xxi=pi/2.
      else
c --> solving in full meridional plane
         xxi=pi
      endif
      xxo=0.
      nx=nelx+1
      do i=1,nx
         xtht(i)=(xxo-xxi)*(((float(i-1)/nelx)))+xxi
	 xtht(i)=cos(xtht(i))
c        xxang=xx(i)
c        write(*,*) i,xtht(i)
      enddo
c     xtht(1)=0.
c     stop ' *TMP*'

c---------- generate grid and connectivity matrix
      call grid1

c---------- compute location of gauss points + some PDE coefficients
      call gaussp

c---------- temporal mesh
      tti=0.
      tto=to
      do i=1,nt+1
         times(i)=(tto-tti)*(((float(i-1)/nt))**1.0)+tti
      enddo
      dt=(to-tti)/float(nsteps*nt)

      write(*,15) tti,tto
   15 format(/,1x,' Time in: ',f8.4,/,1x,' Time out: ',f8.4,/)

c---------- conditions limites ------------------------------------------
c
c           1. condition de dirichlet ("essential boundary conditions")
      numebc=0

c --> conditions at r=r_b: B_phi=0, A=0
      ii=numebc
      do i=1,nelx
         i1=2*icon(i,1)-1
         i2=2*icon(i,1)
         iuebc(ii+1)=i1
         iuebc(ii+2)=i2
c        uebc (ii+1)=(100000./rzc**2)*sin(acos(xx(i)))
         uebc (ii+1)=0.
         uebc (ii+2)=0.
         ii=ii+2
      enddo
      write(*,*) 'numebc = ',ii,
     +           ' ***** conditions along bottom boundary'

c --> conditions along top boundary: A=0
      do i=nnd,nnd-nelx+1,-1
c        i1=2*icon(i,1)-1
         i1=2*i-1
         iuebc(ii+1)=i1
         uebc (ii+1)=0.
         ii=ii+1
      enddo
      write(*,*) 'numebc = ',ii,
     +           ' ***** conditions along top boundary'

c --> conditions within buffer zone: B_phi=0
      do i=nnd,nnd-((nelbuf+1)*(nelx+1))+1,-1
         i2=2*i
         iuebc(ii+1)=i2
         uebc (ii+1)=0.
         ii=ii+1
c write(*,*) i,x(i),y(i)
      enddo
      write(*,*) 'numebc = ',ii,
     +           ' ***** conditions within buffer zone'

c --> conditions along N-pole boundary: B_phi=0, A=0
      do i=nelx+1,nnd,nelx+1
         i1=2*i-1
         i2=2*i
         iuebc(ii+1)=i1
         iuebc(ii+2)=i2
         uebc (ii+1)=0.
         uebc (ii+2)=0.
         ii=ii+2
      enddo
      write(*,*) 'numebc = ',ii,
     +           ' ***** conditions along N-pole'

c --> conditions along equator/S-pole boundary: set by parity
      do i=1,nnd,nelx+1
c following for full meridional plane model: A=0, B=0 along S-polar axis
         if(ipar.eq.0)then
            i1=2*i-1
            i2=2*i
            iuebc(ii+1)=i1
            iuebc(ii+2)=i2
            uebc (ii+1)=0.
            uebc (ii+2)=0.
            ii=ii+2
         else
c following for dipole-like, single quadrant: B=0, dA/dmu=0  on equator
            if(ipar.eq.1)then
               i1=2*i
               iuebc(ii+1)=i1
               uebc (ii+1)=0.
               ii=ii+1
            else
c following for quadrupole-like, single quadrant: A=0, dB/dmu=0 on equator
               i1=2*i-1
               iuebc(ii+1)=i1
               uebc (ii+1)=0.
               ii=ii+1
            endif
         endif
      enddo
      numebc=ii
      write(*,*) 'numebc = ',numebc,
     +           ' ***** conditions along right boundary'

c -- safety check on boundary conditions
      if(numebc.gt.nbc)then
 	write(*,19) numebc,nbc
 	stop ' *MESH*'
      endif
   19 format(/,1x,'Number of boundary conditions    : ',i6,/
     +         1x,'Size of boundary condition arrays: ',i6)
      return
c...........................................................................
      end

c***************************************************************************
      subroutine sfcheck

c===========================================================================
c     Safety checks on the flow fields and equation coefficients
c     Produces a lot of tabular output, basically profiles as a 
c     function of depth and latitude
c
c     Called by : program bldynamo
c===========================================================================
      implicit none
      include  "size.f"
      include  "arrays.global.f"
      include  "arrays.mesh.f"
      include  "arrays.femglobal.f"
c Local
      integer  i,ir1,ir2,ir3,ir4,it1,it2,it3
      real     mu,rr,omg,domgdr,domgdm,remag,drmfac,ur,uth,durdr,
     +         duthdth,alp0,sin1,r1sin1,qdiv,dlnrhodr
c---------------------------------------------------------------------------

c first search for nodal positions
      do i=1,nely
 	if(yrad(i).le. (0.65) .and. yrad(i+1).gt. (0.65)) ir1=i
 	if(yrad(i).le. (0.70) .and. yrad(i+1).gt. (0.70)) ir2=i
 	if(yrad(i).le. (0.85) .and. yrad(i+1).gt. (0.85)) ir3=i
 	if(yrad(i).le. (1.00) .and. yrad(i+1).gt. (1.00)) ir4=i
      enddo
      it1=1
      it2=nelx/2+1
      it3=nelx+1

c -- radial profiles at equator
      write(*,1) it1,acos(xtht(it1)),xtht(it1)
      write(*,10)
      mu=xtht(it1)
      do i=1,ir4+2
 	rr=yrad(i)
        call angvel(mu,rr,omg,domgdr,domgdm)
        call magdif(mu,rr,remag,drmfac)
        call circmd(mu,rr,ur,uth,durdr,duthdth)
        call srcpol(mu,rr,alp0)
	 dlnrhodr=(-0.5/rr**2)/(1./rr-1.)
	 sin1    =sin(acos(mu))
	 r1sin1  =rr*sin1
	 qdiv=r1sin1*dlnrhodr*ur+r1sin1*durdr+2.*sin1*ur
     +       +sin1*duthdth+uth*mu
c       alp0=qdiv
        write(*,11) i,yrad(i),alp0,remag,drmfac,omg,domgdr,domgdm
      enddo

c -- radial profiles at mid-lat
      write(*,2) it2,acos(xtht(it2)),xtht(it2)
      write(*,10)
      mu=xtht(it2)
      do i=1,ir4+2
 	rr=yrad(i)
        call angvel(mu,rr,omg,domgdr,domgdm)
        call magdif(mu,rr,remag,drmfac)
        call circmd(mu,rr,ur,uth,durdr,duthdth)
        call srcpol(mu,rr,alp0)
	 dlnrhodr=(-0.5/rr**2)/(1./rr-1.)
	 sin1    =sin(acos(mu))
	 r1sin1  =rr*sin1
	 qdiv=r1sin1*dlnrhodr*ur+r1sin1*durdr+2.*sin1*ur
     +       +sin1*duthdth+uth*mu
c alp0=qdiv
        write(*,11) i,yrad(i),alp0,remag,drmfac,omg,domgdr,domgdm
      enddo

c -- radial profiles at N-pole
      write(*,3) it3,acos(xtht(it3)),xtht(it3)
      write(*,10)
      mu=xtht(it3)-0.0001
      do i=1,ir4+2
 	rr=yrad(i)
        call angvel(mu,rr,omg,domgdr,domgdm)
        call magdif(mu,rr,remag,drmfac)
        call circmd(mu,rr,ur,uth,durdr,duthdth)
        call srcpol(mu,rr,alp0)
	 dlnrhodr=(-0.5/rr**2)/(1./rr-1.)
	 sin1    =sin(acos(mu))
	 r1sin1  =rr*sin1
	 qdiv=r1sin1*dlnrhodr*ur+r1sin1*durdr+2.*sin1*ur
     +       +sin1*duthdth+uth*mu
	 alp0=qdiv
        write(*,11) i,yrad(i),alp0,remag,drmfac,omg,domgdr,domgdm
      enddo

c -- latitudinal profiles below shear layer (r=0.65)
      write(*,4) ir1,yrad(ir1)
      write(*,12)
      rr=yrad(ir1)
      do i=1,nelx+1
        mu=xtht(i)
        call angvel(mu,rr,omg,domgdr,domgdm)
        call magdif(mu,rr,remag,drmfac)
        call circmd(mu,rr,ur,uth,durdr,duthdth)
        call srcpol(mu,rr,alp0)
	 dlnrhodr=(-0.5/rr**2)/(1./rr-1.)
	 sin1    =sin(acos(mu))
	 r1sin1  =rr*sin1
	 qdiv=r1sin1*dlnrhodr*ur+r1sin1*durdr+2.*sin1*ur
     +       +sin1*duthdth+uth*mu
	 alp0=qdiv
        write(*,11) i,xtht(i),alp0,remag,drmfac,omg,domgdr,domgdm
      enddo

c -- latitudinal profiles at shear layer (r=0.7)
      write(*,5) ir2,yrad(ir2)
      write(*,12)
      rr=yrad(ir2)
      do i=1,nelx+1
        mu=xtht(i)
        call angvel(mu,rr,omg,domgdr,domgdm)
        call magdif(mu,rr,remag,drmfac)
        call circmd(mu,rr,ur,uth,durdr,duthdth)
        call srcpol(mu,rr,alp0)
	 dlnrhodr=(-0.5/rr**2)/(1./rr-1.)
	 sin1    =sin(acos(mu))
	 r1sin1  =rr*sin1
	 qdiv=r1sin1*dlnrhodr*ur+r1sin1*durdr+2.*sin1*ur
     +       +sin1*duthdth+uth*mu
	 alp0=qdiv
        write(*,11) i,xtht(i),alp0,remag,drmfac,omg,domgdr,domgdm
      enddo

c -- latitudinal profiles at mid-envelope (r=0.85)
      write(*,6) ir3,yrad(ir3)
      write(*,12)
      rr=yrad(ir3)
      do i=1,nelx+1
        mu=xtht(i)
        call angvel(mu,rr,omg,domgdr,domgdm)
        call magdif(mu,rr,remag,drmfac)
        call circmd(mu,rr,ur,uth,durdr,duthdth)
        call srcpol(mu,rr,alp0)
	 dlnrhodr=(-0.5/rr**2)/(1./rr-1.)
	 sin1    =sin(acos(mu))
	 r1sin1  =rr*sin1
	 qdiv=r1sin1*dlnrhodr*ur+r1sin1*durdr+2.*sin1*ur
     +       +sin1*duthdth+uth*mu
	 alp0=qdiv
        write(*,11) i,xtht(i),alp0,remag,drmfac,omg,domgdr,domgdm
      enddo

c -- latitudinal profiles at surface (r=0.1)
      write(*,7) ir4,yrad(ir4)
      write(*,12)
      rr=yrad(ir4)
      do i=1,nelx+1
        mu=xtht(i)
        call angvel(mu,rr,omg,domgdr,domgdm)
        call magdif(mu,rr,remag,drmfac)
        call circmd(mu,rr,ur,uth,durdr,duthdth)
        call srcpol(mu,rr,alp0)
	 dlnrhodr=(-0.5/rr**2)/(1./rr-1.)
	 sin1    =sin(acos(mu))
	 r1sin1  =rr*sin1
	 qdiv=r1sin1*dlnrhodr*ur+r1sin1*durdr+2.*sin1*ur
     +       +sin1*duthdth+uth*mu
c        alp0=qdiv
        write(*,11) i,xtht(i),alp0,remag,drmfac,omg,domgdr,domgdm
      enddo

c -- radial profiles at equator
      write(*,1) it2,acos(xtht(it2)),xtht(it2)
      write(*,13)
      mu=xtht(it2)
      do i=1,ir4+2
 	rr=yrad(i)
        call circmd(mu,rr,ur,uth,durdr,duthdth)
	 dlnrhodr=(-0.5/rr**2)/(1./rr-1.)
	 sin1    =sin(acos(mu))
	 r1sin1  =rr*sin1
	 qdiv=r1sin1*dlnrhodr*ur+r1sin1*durdr+2.*sin1*ur
     +       +sin1*duthdth+uth*mu
	 alp0=qdiv
        write(*,11) i,yrad(i),ur,uth,durdr,duthdth,dlnrhodr,qdiv
      enddo

      stop ' *TMP*'
  
      return
c...........................................................................
    1 format(/,1x,'Radial profiles at equator, x-node= ',i4,
     +         ' theta= ',f6.4,2x,'cos(theta)= ',f6.4)
    2 format(/,1x,'Radial profiles at mid-lat, x-node= ',i4,
     +         ' theta= ',f6.4,2x,'cos(theta)= ',f6.4)
    3 format(/,1x,'Radial profiles at N-pole, x-node= ',i4,
     +         ' theta= ',f6.4,2x,'cos(theta)= ',f6.4)
    4 format(/,1x,'Latitudinal profiles below shear layer, y-node= ',i4,
     +         ' r/R= ',f6.4)
    5 format(/,1x,'Latitudinal profiles at shear layer, y-node= ',i4,
     +         ' r/R= ',f6.4)
    6 format(/,1x,'Latitudinal profiles mid-envelope, y-node= ',i4,
     +         ' r/R= ',f6.4)
    7 format(/,1x,'Latitudinal profiles at surface, y-node= ',i4,
     +         ' r/R= ',f6.4)
   10 format(/,7x,
     +       'r/R',7x,'S(r)',6x,'eta(r)',4x,'deta/dr',3x,'Omega',5x,
     +       'dOmg/dr',3x,'dOmg/dmu')
   13 format(/,7x,
     +       'r/R',7x,'u_r',7x,'u_tht',5x,'dur/dr',4x,'duthdth',3x,
     +       'dlnrho/dr',2x,'div rho u')
   12 format(/,7x,
     +       'mu',8x,'S(th)',5x,'eta(th)',3x,'deta/dr',3x,'Omega',5x,
     +       'dOmg/dr',3x,'dOmg/dmu')
   11 format(i5,7e10.2)
      end

c***************************************************************************
      subroutine initcnd

c===========================================================================
c     Calcule ou lit une condition initiale pour les composantes
c     poloidale et toroidale du champ magnetique
c
c     Initialisation de la machinerie FEM
c
c     Called by : program bldynamo
c===========================================================================
      implicit none
      include  "size.f"
      include  "arrays.global.f"
      include  "arrays.mesh.f"
      include  "arrays.femglobal.f"
c Local
      real     fr,q1,q2,tht
      integer  i,i1,i2
      integer  nm0,nel0,nelx0,nely0,nt0,nstep0
      real     comega0,calpha0,reynolds0,etadf0,ybl0,rzc0,to0,fac
      character*80 title0
c---------------------------------------------------------------------------

c---------- initialisations de la machinerie FEM
      call initdt

c---------- now initial condition on magnetic field
c
      if(incr.eq.0)then
c---------- Starting from scratch
	 timein=0.
         do i=1,nnd
            i1=2*i-1
            i2=2*i
            tht=acos(x(i))
c           write(*,*) i,x(i),tht
c --> seed initial poloidal field within domain
	    a(i1)=0.
c --> seed initial toroidal field within domain
	    fr=exp(-400.*(y(i)-rzc-ybl)**2)
            q1=1.e2*fr *sin(2.*tht)
            q2=q1*x(i)
c following for full meridional plane
            if(ipar.eq.0)then
               a(i2)=q1+q2
            else
c following for dipolar
               if(ipar.eq.1)then
                  a(i2)=q1
               else
c following for quadrupolar
                  a(i2)=q2
               endif
            endif
         enddo

      else
c---------- Starting from previous run
         read(7) title0
         read(7) nm0,nel0,nelx0,nely0,nt0,nstep0
         read(7) comega0,calpha0,reynolds0,etadf0,ybl0,rzc0,to0
         timein=to0
c -- some safety check
         if( nelx0.ne.nelx .or. nely0.ne.nely )then
            write(*,1) nelx0,nelx,nely0,nely
            stop ' *INIT*'
         endif
c        if( comega0  .ne.comega   .or. calpha0.ne.calpha .or. 
c    +       reynolds0.ne.reynolds .or. etadf0 .ne.etadf  .or.
c    +       ybl0     .ne.ybl                                  )then
c           write(*,2) comega0,comega,calpha0,calpha,reynolds0,
c    +                 reynolds,etadf,etadf0,ybl0,ybl
c           stop ' *INIT*'
c        endif
         read(7) (xtht(i),i=1,nelx0+1)
         read(7) (yrad(i),i=1,nely0+1)
         read(7) a

c -- facteur multiplicatif pour condition initiale
         fac=1.
         do i=1,nnd
            i1=2*i-1
            i2=2*i
c  if(y(i).le.0.7)then
               a(i1)=fac*a(i1)
               a(i2)=fac*a(i2)
c  endif
         enddo
      endif

      return
c...........................................................................
    1 format(/,1x,'Incompatible Spatial Meshes:',
     +       /,5x,'nelx= ',i3,'/',i3,
     +       /,5x,'nely= ',i3,'/',i3)
    2 format(/,1x,'Incompatible Problem parameters:',
     +       /,5x,'comega  = ',e12.4,'/',e12.4,
     +       /,5x,'calpha  = ',e12.4,'/',e12.4,
     +       /,5x,'reynolds= ',e12.4,'/',e12.4,
     +       /,5x,'etadf   = ',e12.4,'/',e12.4,
     +       /,5x,'ybl     = ',e12.4,'/',e12.4)
c...........................................................................
      end

c***************************************************************************
      subroutine formm

c===========================================================================
c  form the system matrices from element matrices contributions            
c  (See Burnett, chap. 5)                                          
c                                                                          
c  This is an adaptation of equivalent routine in Burnett's code          
c  modified to handle non-symmetric system matrices                        
c
c  Called by: program bldynamo
c  Calls    : subroutine elem23c
c===========================================================================
      implicit none
      include  "size.f"
      include  "arrays.mesh.f"
      include  "arrays.femglobal.f"
      include  "arrays.femlocal.f"
c Local
      integer  i,j,ielno,i1,i2,j1,j2,ii1,ii2,jj1,jj2,ish,ishm,ishp,
     +         iflag
      real     onemth
c---------------------------------------------------------------------------
c---------- initialize arrays for the system matrices
c     write(*,*) ' FORM: entering '
      onemth=1.-theta

      do j=1,nm
         do i=1,mfb
            k(i,j)=0.
            ck(i,j)=0.
         enddo
      enddo
c******************** loop over all elements
      iflag=1
      do ielno=1,nel
c     write(*,*) ' FORM: working on element #',ielno

c---------- form element matrices
      call elem23c(ielno,iflag)

c---------- assemble the system matrices from the element matrices
      do i=1,nd
	 i1=2*i-1
	 i2=2*i
	 ii1=2*icon(ielno,i)-1
	 ii2=ii1+1
         do j=1,nd
	    j1=2*j-1
	    j2=2*j
	    jj1=2*icon(ielno,j)-1
	    jj2=jj1+1
	    ish=ii1+mhb-jj1
	    ishm =ish-1
	    ishp =ish+1
	    k(ish  ,jj1) = k(ish  ,jj1)+ce(i1,j1)/dt +theta *ke(i1,j1)
	    k(ishm ,jj2) = k(ishm ,jj2)+ce(i1,j2)/dt +theta *ke(i1,j2)
	    k(ishp ,jj1) = k(ishp ,jj1)+ce(i2,j1)/dt +theta *ke(i2,j1)
	    k(ish  ,jj2) = k(ish  ,jj2)+ce(i2,j2)/dt +theta *ke(i2,j2)
	    ck(ish  ,jj1)=ck(ish  ,jj1)+ce(i1,j1)/dt -onemth*ke(i1,j1)
	    ck(ishm ,jj2)=ck(ishm ,jj2)+ce(i1,j2)/dt -onemth*ke(i1,j2)
	    ck(ishp ,jj1)=ck(ishp ,jj1)+ce(i2,j1)/dt -onemth*ke(i2,j1)
	    ck(ish  ,jj2)=ck(ish  ,jj2)+ce(i2,j2)/dt -onemth*ke(i2,j2)
         enddo
      enddo

      enddo

c     write(*,*) ' FORM: exiting '

      return
      end

c***************************************************************************
      subroutine formf

c===========================================================================
c  form the RHS vector by assembling individual element's contributions
c  (See Burnett, chap. 5)                                          
c                                                                          
c  Called by: program bldynamo, subroutine calfun
c  Calls    : subroutine elem23f
c===========================================================================
      implicit none
      include  "size.f"
      include  "arrays.mesh.f"
      include  "arrays.femglobal.f"
      include  "arrays.femlocal.f"
c Local
      integer  i,ielno,i1,i2,ii1,ii2,iflag
      real     onemth
c---------------------------------------------------------------------------
c---------- initialize arrays for the system matrices
c     write(*,*) ' FORM: entering '
      onemth=1.-theta

      do i=1,nm
         f(i)=0.
      enddo
c******************** loop over all elements
      iflag=1
      do ielno=1,nel
c     write(*,*) ' FORM: working on element #',ielno

c---------- form element matrices
      call elem23f(ielno,iflag)

c---------- assemble the system matrices from the element matrices
      do i=1,nd
	 i1=2*i-1
	 i2=2*i
	 ii1=2*icon(ielno,i)-1
	 ii2=ii1+1
         f(ii1)=f(ii1)+fe(i1)
         f(ii2)=f(ii2)+fe(i2)
      enddo

      enddo

c     write(*,*) ' FORM: exiting '

      return
      end

c***************************************************************************
      subroutine aplybc
 
c===========================================================================
c  Modifies system matrices to incorporate Dirichlet conditions            
c  (See Burnett, chap. 7)                                                  
c                                                                          
c  This is an adaptation of equivalent routine in Burnett's code           
c  modified to handle non-symmetric system matrices                        
c
c  Called by: program bldynamo
c===========================================================================
      implicit none
      include  "size.f"
      include  "arrays.mesh.f"
      include  "arrays.femglobal.f"
      include  "arrays.femlocal.f"
c Local
      integer  i,j,jj,ibc,jpmb,jshift,ishp,ishm,jrp,jrm
      real     sum1,sum2
 
c******************** apply mixed boundary conditions, by modifying
c                     diagonal of stiffness matrix
c     write(*,*) ' APLYBC: entering '
c     if(nummbc.le.0) goto 400
c     do 300 i=1,nm
c        k(mband,i)=k(mband,i)-umbc(i)
c 300 continue
c******************** apply the essential boundary conditions, if any,
c                     by imposing constraint equations
  400 if(numebc.le.0) return
c---------- calculate step change in essential boundary conditions
c              (for ramping in calfun)
      do i=1,nm
         u0(i)=0.
         delu(i)=0.
         keffdu(i)=0.
         keffu0(i)=0.
      enddo
      do ibc=1,numebc
         i=iuebc(ibc)
         u0(i)=a(i)
         delu(i)=(uebc(ibc)-a(i))/float(nsteps)
c        write(*,*) ibc,' node # ',jj, ' BC= ',uebc(ibc),u0(i)
      enddo
c********** enforce dirichlet B.C. by applying constraint equations to system
c---------- 1. calculate modified load vector [see Burnett, eq.(11.116)]

      do j=1,nm
	 sum1=0.
	 sum2=0.
	 jpmb=j+mhb
	 do i=max(1,jpmb-nm),min(mfb,jpmb-1)
	    jshift=jpmb-i
	    sum1=sum1+k(i,jshift)*delu(jshift)
	    sum2=sum2+k(i,jshift)*  u0(jshift)
         enddo
	 keffdu(j)=sum1
	 keffu0(j)=sum2
      enddo

c---------- 2. modify system matrix
      do ibc=1,numebc
         jj=iuebc(ibc)
c        write(*,*) ' working on BC # ',ibc,' node # ',jj
c --> zero matrix column
         do i=1,mfb
            k(i,jj)=0.
         enddo
c --> zero matrix "line" (shifted version)
         do i=2,mhb
            ishp=mhb-i+1
            ishm=mhb+i-1
            jrp=jj+i-1
            jrm=jj-i+1
            if(jrm.ge.1)  k(ishm,jrm)=0.
            if(jrp.le.nm) k(ishp,jrp)=0.
         enddo
c --> set corresponding diagonal element of system matrix to one
	 k(mhb,jj)=1.
      enddo
c     write(*,*) ' APLYBC: exiting '
      return
c------------------------------------
      end

c******************************************************************************
      subroutine calfun(icall)
 
c===========================================================================
c  Solve the system equations for the function values, using
c  the Theta-method, an implicit one-step time-stepping scheme
c  (See Burnett, sec. 11.3, 11.4)                                          *
c                                                                          *
c  This is an adaptation of equivalent routine in Burnett's code           *
c  modified to handle non-symmetric system matrices                        *
c
c  Called by: program bldynamo
c  Calls    : subroutines  tribu, rhsbu
c===========================================================================
      implicit none
      include   "size.f"
      include   "arrays.global.f"
      include   "arrays.mesh.f"
      include   "arrays.femglobal.f"
c Input
      integer   icall
c Local
      integer   i,j,jpmb,jshift,ieb,n
      real      onemth,rns,rn,cka,time,feff,quant
      dimension feff(nm)
      real           urand
      external       urand
 
c-----------------------------------
c******************** implicit initial-boundary-value problem
c     write(*,*) ' CALFUN: entering '
c     write(*,*) 'nm = ',nm
c---------- save load vector and function at start of interval
      time=times(icall)
      do i=1,nm
c        fnm1(i)=f0(i)
         fnm1(i)=f(i)
         anm1(i)=a(i)
      enddo
      onemth=1.-theta
      rns=float(nsteps)

c---------- forward reduction of effective stiffness matrix (keff)
      call tribu(k)

c---------- solve recurrence relation over all time steps
      do n=1,nsteps

         time=time+dt
         rn=float(n)
         write(*,*) ' CALFUN: n,time = ',n,time
c---------- calculate nat+int loads at step n by ramping
c        do 420 i=1,nm
c           fn(i)=f0(i)+(rn/rns)*(f(i)-f0(i))
c 420    continue

c---------- compute stochastic multiplicative factor for alpha-effect
c        quant=float(n)/2.-float(n/2)
         quant=float(n)/5.-float(n/5)
c        if(quant.le.0.) alpsfac=4.*urand()-1.
c        if(quant.le.0.) alpsfac=2.*urand()
c        if(quant.le.0.) alpsfac=0.5+urand()
c        if(quant.le.0.) alpsfac=0.75+0.5*urand()
         if(quant.le.0.)then
            alpsfac=0.0+2.0*urand()
            alpsfac2=0.0+2.0*urand()
         endif
c        if(quant.le.0.) alpsfac=2.*urand()
c        alpsfac=0.5+urand()
            write(*,*) alpsfac,alpsfac2
         alpsfac=1.
         alpsfac2=1.
c---------- calculate RHS vector at this sub-timestep
	 call formf

         do i=1,nm
            fn(i)=f(i)
         enddo

c---------- calculate effective load vector (feff) at step n
      do j=1,nm
	 cka=0.
	 jpmb=j+mhb
	 do i=max(1,jpmb-nm),min(mfb,jpmb-1)
	    jshift=jpmb-i
	    cka=cka+ck(i,jshift)*anm1(jshift)
         enddo
	 feff(j)=onemth*fnm1(j)+theta*fn(j)+cka
      enddo
c---------- modify effective load vector for ramped essential bc
c --> insertion du "if" pour eviter keffu0 indefini si numebc=0  7/4/88
c        if(numebc.gt.0) then
            do i=1,nm
               feff(i)=feff(i)-keffu0(i)-rn*keffdu(i)
            enddo
            do ieb=1,numebc
               i=iuebc(ieb)
               feff(i)=anm1(i)+delu(i)
            enddo
c        end if

c---------- reduction of load vector: backsubstitution
         call rhsbu(k,feff,a)

c---------- compute magnetic energy and such
         call postpro(n)

c---------- adjust meridional flow speed
         reynld2=reynld2+dt*(aaa*lorentz-bbb*reynld2)
         write(*,*) ' reynld2= ',reynld2

c---------- output nodal time series
         call writout(n,time)

c---------- additive noise
         call addnoise

c---------- save values for next step
         do i=1,nm
            anm1(i)=a (i)
            fnm1(i)=fn(i)
         enddo

      enddo

c---------- save nat+int loads at end of interval
      do i=1,nm
         f0(i)=f(i)
      enddo
c     write(*,*) ' CALFUN: exiting '
      return
c-------------------------------------
      end

c**************************************************************************
      subroutine initdt
c========================================================================= 
c  reads initial data, sets limits and integration points (Internal to FEM)
c
c  Called by : subroutine initcnd
c========================================================================= 
      implicit none
      include "size.f"
      include "arrays.femlocal.f"
c Local
      integer n,l
c-------------------------------------
c     do i=1,30
c        intdef(i)=0
c        ioptau(i)=0
c     enddo
c     nnd(23)=4
c     nnd(24)=8
c     intdef(23)=2
c     intdef(24)=3
c     ioptau(23)=1
c     ioptau(24)=4
 
c---------- gauss points and weights for gauss-legendre quadrature
      do n=1,5
         do l=1,5
            gp(n,l)=0.
            wt(n,l)=0.
         enddo
      enddo
c---------- for one gauss point
      gp(1,1)=0.
      wt(1,1)=2.
c---------- for two gauss points
      gp(2,1)=-1./sqrt(3.)
      gp(2,2)=-gp(2,1)
      wt(2,1)=1.
      wt(2,2)=1.
c---------- for three gauss points
      gp(3,1)=-sqrt(3./5.)
      gp(3,2)=0.
      gp(3,3)=-gp(3,1)
      wt(3,1)=5./9.
      wt(3,2)=8./9.
      wt(3,3)=wt(3,1)
c---------- for four gauss points
      gp(4,1)=-sqrt((15.+2.*sqrt(30.))/35.)
      gp(4,2)=-sqrt((15.-2.*sqrt(30.))/35.)
      gp(4,3)=-gp(4,2)
      gp(4,4)=-gp(4,1)
      wt(4,1)=49./(6.*(18.+sqrt(30.)))
      wt(4,2)=49./(6.*(18.-sqrt(30.)))
      wt(4,3)=wt(4,2)
      wt(4,4)=wt(4,1)

      return
c------------------------------------------------------------------------------
      end

c***************************************************************************
      subroutine writout(icall,time)

c***************************************************************************
c   outputs results of selected time step                                  *
c         unit8 --> formatted for plotting                                 *
c         unit* --> tabular form                                           *
c***************************************************************************
      implicit none
      include  "size.f"
      include  "arrays.global.f"
      include  "arrays.mesh.f"
      include  "arrays.femglobal.f"
c Input
      integer   icall
      real      time
c Local
      integer   i,j
      real      au,ab,omgn,urn,uthn,etan,dum1,dum2
      dimension au(nnd),ab(nnd),omgn(nnd),urn(nnd),uthn(nnd),etan(nnd)
c---------------------------------------------------------------------------
c---------------------------------------------------------------------------

c---------- on first time step, output header info
      if(icall.eq.0)then
         write(8) title
         write(8) nm,nel,nelx,nely,nt,nsteps
         write(8) comega,calpha,reynolds,etadf,ybl,rzc,to
	 write(8) (xtht(i),i=1,nelx+1)
         write(8) (yrad(i),i=1,nely+1)
         do i=1,nnd
c -- compute angular velocity, diffusivity, and circulation at nodes
            call angvel(x(i),y(i),omgn(i),        dum1,dum2)
            call magdif(x(i),y(i),etan(i),             dum1)
            call circmd(x(i),y(i), urn(i),uthn(i),dum1,dum2)
	 enddo
	 write(8) (omgn(i),i=1,nnd)
	 write(8) ( urn(i),i=1,nnd)
	 write(8) (uthn(i),i=1,nnd)
	 write(8) (etan(i),i=1,nnd)
      end if

c---------- separate poloidal and toroidal fields
      do i=1,nnd
	 au(i)=a(2*i-1)
         ab(i)=a(2*i)
      enddo

      write(8) icall
      write(8) time,ebtot,ebttot,b1max,b2max,reynld2,lorentz
      write(8) (au(j),j=1,nnd)
      write(8) (ab(j),j=1,nnd)

      return
c..........................................................................
      end

c***************************************************************************
      subroutine addnoise

c***************************************************************************
c   adds stochastic noise on poloidla field at surface
c***************************************************************************
      implicit none
      include  "size.f"
      include  "arrays.global.f"
      include  "arrays.mesh.f"
      include  "arrays.femglobal.f"
c Local
      integer   ii,iel,nelnoise
      real      urand
      external  urand
c---------------------------------------------------------------------------

c---------- introduce noise at surface
      nelnoise=nelx*nely/4
      do iel=ielnoint,ielnoint+nelnoise-1
         ii=icon(iel,2)
         if(x(ii).gt.0.and.x(ii).lt.1.0)then
	 a(2*ii)=a(2*ii)+addnamp*(2.*urand()-1.)
         a(2*ii-1)=a(2*ii-1)+addnamp*(2.*urand()-1.)/nelx
         endif
c        a(2*i)=a(2*i)+addnamp*2.*(urand()-1.)
      enddo

      return
c..........................................................................
      end

c***************************************************************************
      subroutine stlast

c***************************************************************************
c   stores solutions of last step, for future usage as initial condition   *
c***************************************************************************
      implicit none
      include  "size.f"
      include  "arrays.global.f"
      include  "arrays.mesh.f"
      include  "arrays.femglobal.f"
c Local
      integer i
c---------------------------------------------------------------------------
c---------- output grid values
      write(9) title
      write(9) nm,nel,nelx,nely,nt,nsteps
      write(9) comega,calpha,reynolds,etadf,ybl,rzc,to
      write(9) (xtht(i),i=1,nelx+1)
      write(9) (yrad(i),i=1,nely+1)

      write(9) a

      return
      end
     
c***************************************************************************
      subroutine postpro(it)

c===========================================================================
c     Postprocessing: computes magnetic energy and field strength in domain  
c     and portions thereof, at timestep it
c
c     ebttot: Magnetic energy based on B_phi only, in whole domain
c     ebtot : Magnetic energy based on full B, in whole domain
c     b1max : Maximal toroidal field strength in r/R<rzc+ybl
c     b2max : Maximal toroidal field strength in r/R>rzc+ybl
c
c     Called by: program bldynamo
c===========================================================================
      implicit none

      include "size.f"
      include "arrays.global.f"
      include "arrays.mesh.f"
      include "arrays.femglobal.f"

c Input
      integer        it
c Local
      integer        iflag,ielno,i,ii,iell
      real           pi,bphi
      data           pi/3.1415926536/
c-----------------------------------
c     write(*,*) ' POST: entering '
c********** loop over all elements to compute magnetic energy
      iflag=2
      ebtot=0.
      ebttot=0.
      
      do ielno=1,nel
         call elem23f(ielno,iflag)
         ebtot=ebtot+ebt
         ebttot=ebttot+ebtt
      enddo
      ebttot=ebttot/(8.*pi)
      ebtot=ebtot/(8.*pi)

c********** now look for maximum toroidal field strength
c --> at interface
      iell=2*icon(ielnoint+nelx/2,1)
      lorentz=a(iell-1)*a(iell)

      b1max=-1.e10
      b2max=-1.e10
      do i=1,nnd
c --> in core and shear layer
	 if(y(i).ge.0..and.y(i).lt.rzc+ybl)then
	    ii=2*i
	    bphi=abs(a(ii))
	    if(bphi.gt.b1max) b1max=bphi
	 endif
c --> in envelope
	 if(y(i).gt.(rzc+ybl).and.y(i).le.1.0)then
	    ii=2*i
	    bphi=abs(a(ii))
	    if(bphi.gt.b2max) b2max=bphi
	 endif
      enddo

c********** some output
      write(*,210)
      write(*,211) it,times(ijac),ebtot,ebttot,b2max,b1max

      return
c...................................................................
  210 format(/,'SOLUTION CHARACTERISTICS AT THIS STEP:',/)
  211 format(5x,'Time step #                  :',i6,/,
     +       5x,'Simulation time              :',e12.4,/,
     +       5x,'Total Mag. Energy Content    :',e12.4,/,
     +       5x,'Mag. Energy Content (B_y)    :',e12.4,/,
     +       5x,'Max(B_phi) envelope          :',e12.4,/,
     +       5x,'Max(B_phi) shear layer       :',e12.4)

      end

c***************************************************************************
      subroutine tribu(gk)

c===========================================================================
c     triangularizes shifted matrix
c     (modified alg., after Golub & Van Loan, p. 152)
c===========================================================================
      implicit none
      include "size.f"
c Local
      integer i,j,k,kpj,imj
      real      gk,c
      dimension gk(mfb,nm)

      do k=1,nm-1
	 if(gk(mhb,k).eq.0.)then
            write(*,100) k,gk(mhb,k)
            stop ' *TRIBU*'
	 endif
	 c=1./gk(mhb,k)
	 do i=mhb+1,mfb
            gk(i,k)=c*gk(i,k)
         enddo
         do j=1,min(mhb-1,nm-k)
	    kpj=k+j
	    c=gk(mhb-j,kpj)
	    if(c.ne.0.)then
CDIR$ IVDEP
	    do i=mhb+1,min(mfb,nm-k+mhb)
	       imj=i-j
	       gk(imj,kpj)=gk(imj,kpj) -c*gk(i,k)
            enddo
	    end if
         enddo
      enddo
      return
c...........................................................................
  100 format(//,'  set of equations may be singular..',/,
     +'   diagonal term of equation',i5,'  at tribu is equal to',
     +1pe15.8)
      end

c***************************************************************************
      subroutine rhsbu(gk,gf,u)

c===========================================================================
c     forward reduction & backsubstitution of triangularized system matrix
c     (modified alg., after Golub & Van Loan, p. 152)
c===========================================================================
      implicit none
      include "size.f"
c Local
      integer i,j,mhbmj,ishift
      real      gk,gf,u
      dimension gk(mfb,nm),gf(nm),u(nm)

c---------- forward reduction
      do j=1,nm
	 if(gk(mhb,j).eq.0.)then
            write(*,100) j,gk(mhb,j)
            stop ' *RHSBU*'
	 endif
	 mhbmj=mhb-j
	 do i=j+1,min(nm,j+mhb-1)
	    ishift=i+mhbmj
	    gf(i)=gf(i)-gk(ishift,j)*gf(j)
         enddo
      enddo

c---------- backsubstitution 

      do j=nm,1,-1
c        gf(j)=gf(j)/gk(mhb,j)
	 u(j)=gf(j)/gk(mhb,j)
	 mhbmj=mhb-j
	 do i=max(1,j-mhb+1),j-1
	    ishift=i+mhbmj
c           gf(i)=gf(i)-gk(ishift,j)*gf(j)
            gf(i)=gf(i)-gk(ishift,j)*u(j)
         enddo
c        u(j)=gf(j)
      enddo

      return
c.........................................................................
  100 format(//,'  set of equations are singular.',/,
     +'   diagonal term of equation',i5,'  at tribu is equal to',
     +1pe15.8)
      end

c********************************************************************
c      INCLUDE FILES
c
c      mod.angvel.f    --> Solar-like differential rotation
c      mod.magdif.f    --> Depth-dependent magnetic diffusivity profile
c      mod.circmd.f    --> Solar-like envelope meridional circulation
c      mod.srcpol.f    --> Babcock-Leighton source term
c      numrecp.f       --> needed Numerical Recipes routines
c      random.f        --> needed for stochastically forced version
c********************************************************************

      include "random.f"
      include "numrecp.f"
      include "mod.angvel.f"
      include "mod.magdif.f"
c     include "mod.srcpol.ao.f"
      include "mod.srcpol.f"
      include "mod.circmd.f"
