program tacrni
      implicit none
      save    

c..tests the crank-nickolson 1d diffusion solver

c..test problems from 
c..g.d. smith     "numerical solution of partial differential
c..               equations, finite difference methods"
c..               oxford univ. press 1978
c..w.e. schiesser "the numerical method of lines"
c..               academic press 1991


c..declarations for the solver
      external         pdecof
      integer          ngridmax,ngrid,init,j1,jn,ncycle,itemp,iter,ifail
      parameter        (ngridmax=400)
      double precision xgrid(ngridmax),temper(ngridmax),dt,dtnext,
     1                 neu1l,xlpow,neu2l,neu1r,xrpow,neu2r,
     2                 tolx,tempf,flux(ngridmax),term1(ngridmax),
     3                 term2(ngridmax),dttemp

c..math constants
      double precision pi
      parameter        (pi = 3.1415926535897932384d0)


c..local variables
      external         aroot
      integer          i,j,n,njj,icase,niter,nb
      double precision dx,time,xlength,xtime,sum,te(ngridmax),
     1                 errmax,percent,maxper,xwhere,tdum,
     2                 alfan,alfan2,zbrent,aroot,lo,hi,rtol,
     3                 xb1(1000),xb2(1000)
      parameter        (rtol = 1.0d-12)


c..popular format statements
 20   format(1x,a,0p12f8.4,/,(8(3x,0p12f8.4)))
 21   format(1x,a,0p12f8.4)
 22   format(1x,'time=',1pe11.3,'  dt=',1pe11.3,'  iter=',i6)
 23   format(1x,a,1pe12.4,a,1pe12.4)





c..loop over the test cases
      do icase=6,6



c..convergence parameters and initializations
       init   = 1
       tolx   = 1.0e-6
       time   = 0.0
       dtnext = 0.0
       tempf  = 0.05


c..schiesser page 16, simple dirchlett
       if (icase .eq. 1) then
        j1      = 0
        jn      = 0
        neu1l   = 0.0d0  
        xlpow   = 1.0d0
        neu2l   = 0.0d0  
        neu1r   = 0.0d0
        xrpow   = 1.0d0  
        neu2r   = 0.0d0  
        xlength = 1.0d0
        njj     = 1
        ngrid   = 10*njj + 1
        dx      = xlength/float(ngrid-1)
        dt      = 1.0d-6
        do i=1,ngrid
         xgrid(i)  = float(i-1) * dx
         temper(i) = sin(pi*xgrid(i)/xlength)
        enddo
        write(6,*) 
        write(6,*) 'example from schiesser page 16 ngrid=',ngrid


c..schiesser page 28, simple neumann 
       else if (icase .eq. 2) then
        j1      = 1
        jn      = 1
        neu1l   = 0.0d0 
        xlpow   = 1.0d0 
        neu2l   = 0.0d0  
        neu1r   = 0.0d0 
        xrpow   = 1.0d0 
        neu2r   = 0.0d0  
        xlength = 1.0d0
        njj     = 1
        ngrid   = 10*njj + 1
        dx      = xlength/float(ngrid-1)
        dt      = 1.0d-6
        do i=1,ngrid
         xgrid(i) = float(i-1) * dx
         temper(i) = cos(pi*xgrid(i)/xlength) 
        enddo
        write(6,*) 
        write(6,*) 'example from schiesser page 28 ngrid=',ngrid



c..smith page 21, initial discontinuity at x=0.5
       else if (icase .eq. 3) then
        j1      = 0
        jn      = 0
        neu1l   = 0.0d0
        xlpow   = 1.0d0
        neu2l   = 0.0d0
        neu1r   = 0.0d0
        xrpow   = 1.0d0
        neu2r   = 0.0d0
        xlength = 1.0d0
        njj     = 1
        ngrid   = 10*njj + 1
        dx      = xlength/float(ngrid-1)
        dt      = 1.0d-6
        do i=1,ngrid
         xgrid(i)    = float(i-1) * dx
         if (xgrid(i) .le. 0.5) then
          temper(i) = 2.0d0 * xgrid(i)
         else
          temper(i) = 2.0d0 * (1.0d0 - xgrid(i))
         end if
        enddo
        write(6,*) 
        write(6,*) 'example from smith page 21 ngrid=',ngrid



c..smith page 35, exercises nonzero neumann boundaries
       else if (icase .eq. 4) then
        j1      = 1
        jn      = 1
        neu1l   = 1.0d0
        xlpow   = 1.0d0
        neu2l   = 0.0d0
        neu1r   = -1.0d0
        xrpow   = 1.0d0
        neu2r   = 0.0d0
        xlength = 1.0d0
        njj     = 1
        ngrid   = 10*njj + 1
        dx      = xlength/float(ngrid-1)
        dt      = 1.0d-6
        do i=1,ngrid
         xgrid(i)    = float(i-1) * dx
         temper(i) = 1.0
        enddo
        write(6,*) 
        write(6,*) 'example from smith page 35 ngrid=',ngrid



c..smith page 62, different neumann conditions on each end
       else if (icase .eq. 5) then
        j1      = 1
        jn      = 1
        neu1l   = 0.0d0
        xlpow   = 1.0d0
        neu2l   = 0.0d0
        neu1r   = 0.0d0
        xrpow   = 1.0d0
        neu2r   = 1.0d0
        xlength = 0.5d0
        njj     = 1
        ngrid   = 5*njj + 1
        dx      = xlength/float(ngrid-1)
        dt      = 1.0d-6
        do i=1,ngrid
         xgrid(i)  = float(i-1) * dx
         temper(i) = 0.0d0
        enddo
        write(6,*) 
        write(6,*) 'example from smith page 62 ngrid=',ngrid


c..scheisser page 54, a nonlinear boundary
       else if (icase .eq. 6) then
        j1      = 0
        jn      = 1
        neu1l   = 0.0d0  
        xlpow   = 1.0d0
        neu2l   = 0.0d0  
        neu1r   = -1.0d0
        xrpow   = 4.0d0
        neu2r   = 1.0d0  
        xlength = 1.0d0
        njj     = 10
        ngrid   = 5*njj + 1
        dx      = xlength/float(ngrid-1)
        dt      = 1.0d-6
        do i=1,ngrid
         xgrid(i)  = float(i-1) * dx
         temper(i) = 0.0d0
        enddo
        temper(1) = 1.0d0
        write(6,*) 
        write(6,*) 'example from schiesser page 54 ngrid=',ngrid


c..end of test problem if
       end if




c..here is the evolution loop
09     continue


c..figure the exact solution
c..scheisser page 16, simple dirchlett
       if (icase .eq. 1) then
        do i=1,ngrid,njj
         te(i) = exp((-pi**2/xlength**2)*time)
     1           * sin(pi*xgrid(i)/xlength)
        enddo


c..scheisser page 28, simple neumann
       else if (icase .eq. 2) then
        do i=1,ngrid,njj
         te(i) = exp((-pi**2/xlength**2)*time)
     1           * cos(pi*xgrid(i)/xlength)
        enddo


c..smith page 21, initial discontinuity at x=0.5
       else if (icase .eq. 3) then
        do i=1,ngrid,njj
         sum = 0.0d0
         do n=1,100
          sum = sum + sin(0.5d0*pi*n)*sin(n*pi*xgrid(i))/n**2
     1                * exp(-(n*pi)**2*time)
         enddo
         te(i) = 8.0d0/pi**2 * sum
        enddo


c..smith page 35, exercises non-zero neumann consitiond
c..get the positive roots of alfa*tan(alfa)=0.5
c..ignore the sign changes at n*pi/2
       else if (icase .eq. 4) then
        lo  = 0.0
        hi  = 1000.0d0
        nb  = 5000
        call zbrak(aroot,lo,hi,5000,xb1,xb2,nb)
        do i=1,ngrid,njj
         sum = 0.0d0
         do n=1,nb/2,2
          alfan  = zbrent(aroot,xb1(n),xb2(n),rtol,niter)
          alfan2 = alfan*alfan
          sum    = sum + 
     1        exp(-4.0d0*alfan2*time)
     2        / (cos(alfan)*(3.0d0+4.0d0*alfan2))
     2          * cos(2.0d0*alfan*(xgrid(i)-0.5d0))
         enddo
         te(i) = 4.0d0 * sum
        enddo


c..smith page 62, different neumann conditions on each end
       else if (icase .eq. 5) then
        do i=1,ngrid,njj
         sum = 0.0d0
         do n=1,100
          sum    = sum + 
     1        (-1.0d0)**n * exp(-4.0d0*(pi*n)**2*time)
     2        * cos(2.0d0*n*pi*xgrid(i))/n**2
         enddo
         te(i) = 2.0d0*time + 0.5d0* ( 
     1           (12.0d0*xgrid(i)**2 - 1.0d0)/6.0d0
     2            - 2.0d0/pi**2 * sum)
        enddo

c..or no analytic solution
       else
        do i=1,ngrid,njj
         te(i) = 0.0d0
        enddo

       end if



c..max percentage difference between the numerical and analytic solution
       errmax = -1.0e30
       do i=1,ngrid,njj
        if (abs(te(i)) .gt. 1.0d-12) then
         percent = ( (te(i) - temper(i))/te(i) ) * 100.0d0
         errmax = max(errmax,abs(percent))
         if (errmax .eq. abs(percent)) then
          xwhere = xgrid(i)
          maxper = percent
         end if
        end if
       enddo



c..write out the grid, numerical solution and exact solution
       write(6,22) time,dt,iter
       if (ngrid .le. 12) then
        write(6,21) 'x =',(xgrid(j),j=1,ngrid,njj)
        write(6,21) 't =',(temper(j),j=1,ngrid,njj)
        write(6,21) 'te=',(te(j),j=1,ngrid,njj)
       else
        write(6,20) 'x =',(xgrid(j),j=1,ngrid,njj)
        write(6,20) 't =',(temper(j),j=1,ngrid,njj)
        write(6,21) 'te=',(te(j),j=1,ngrid,njj)
       end if
       write(6,23) 'max % error=',maxper,' occurs at x=',xwhere



c..get how many time steps to take
       write(6,*)
       write(6,*) 'time to evolve (0 next case; -1 to stop) =>'
       read(5,*) xtime
       if (xtime .eq. 0.0) goto 200
       if (xtime .eq. -1.0) stop 'normal termination'



c..evolve for the given amount of xtime
       tdum = 0.0d0
       do ncycle=1,100000

        call acrni(ngrid,xgrid,temper,init,dt,dtnext,
     1            j1,jn,neu1l,xlpow,neu2l,neu1r,xrpow,neu2r,
     2            ncycle,tolx,tempf,
     3            flux,term1,term2,dttemp,itemp,iter,ifail,
     4            pdecof)

        time   = time + dt
        tdum   = tdum + dt
        dtnext = dttemp
        if (dtnext .eq. 0.0) dtnext = dt
        if (tdum + dtnext .gt. xtime) dtnext = xtime - tdum
        if (dtnext .le. 1.0d-14) goto 100
       enddo
 100   continue


c..back for another evolution
       goto 09



c..or start another test case
200   continue
      enddo
      stop 'normal termination'
      end







      subroutine pdecof(nn,temp,
     1                den,drhodt,
     2                ener,denerdt,
     3                sigma,dsigmadt,
     4                sdot,dsdotdt,
     5                flim,dflimdt)

      implicit none  
      save

c..this routine sets the coeficients for the diffusion solver

c..input:
c..nn       = number of grid points
c..temp     = temperature at each gridpoint

c..output:
c..den      = density
c..drhodt   = d(rho)/d(temperature)
c..ener     = specific internal energy
c..denerdt  = d(ener)/d(temperature) 
c..sigma    = thermal diffusion coefficient
c..dsigmadt = d(sigma)/d(temperature)
c..sdot     = energy generation rate
c..dsdotdt  = d(sdot)/d(temperature)
c..flim     = flux limiter coefficient
c..dflimdt  = d(flim)/d(temperature)


c..declare the pass
      integer          nn
      double precision temp(nn),
     1                 den(nn),drhodt(nn),
     2                 ener(nn),denerdt(nn),
     3                 sigma(nn),dsigmadt(nn),
     4                 sdot(nn),dsdotdt(nn),
     5                 flim(nn),dflimdt(nn)

c..local variables
      integer          i


c..set the parameters for the test problems 
      do i=1,nn
       den(i)      = 1.0d0
       drhodt(i)   = 0.0d0 
       ener(i)     = temp(i)
       denerdt(i)  = 1.0d0
       sigma(i)    = 1.0d0
       dsigmadt(i) = 0.0d0
       sdot(i)     = 0.0d0
       dsdotdt(i)  = 0.0d0
       flim(i)     = 1.0d99
       dflimdt(i)  = 0.0d0
      enddo
      return
      end





      subroutine acrni(ngrid,x,tnew,init,dt,dtnext,
     1                 j1,jn,neu1l,xlpow,neu2l,neu1r,xrpow,neu2r,
     2                 ncycle,tolx,tempf,
     3                 flux,term1,term2,dttemp,itemp,iter,ifail,
     4                 func)
      implicit none
      save


c..this routine takes a single time step for the 1d diffusion equation:
c..
c..                de(T)/dt = 1/rho * d/dx (sigma(T) * dT/dx) + s(T)
c..
c..where 
c..T        = temperature
c..e(T)     = specific internal energy 
c..rho(T)   = density 
c..sigma(T) = thermal diffusion coefficient
c..s(T)     = energy generation rate

c..alternatively, e(t) could be the enthalpy or just the temperature,
c..in which case the meaning of rho(T) may change. 


c..a newton-raphson iteration on the crank-nicolson algorithm is used.
c..for constant coefficients on a uniform mesh, this solver is second-order
c..accurate in space and time. for variable coefficients or nonuniform grids,
c..i find that the solutions are less than second order,
c..but better than first order, in space.

c..input:
c..ngrid   = number of grid points 
c..x       = the spatial grid 
c..tnew    = temperature on the grid
c..init    = initialization flag (= 1 for a complete initialization)
c..dt      = last timestep succesfully completed, or first timestep to try  
c..dtnext  = timestep to try
c..j1 & jn = boundary condition flags. dirchlett boundary are used at 
c..          u(j=1) if j1=0 and at u(j=n) if jn=0. neumann boundary are used at
c..          u(j=1) if j1=1 and at u(j=n) if jn=1. acceptable forms for the 
c..          neumann boundaries are 
c..          du/dx = neu1l * u**xlpow + neu2l at the j=1 boundary
c..          du/dx = neu1r * u**xrpow + neu2r at the j=n boundary 
c..          where neu1l, neu2l, neu1r and neu2r are constants.
c..ncycle  = cycle number, only used for error message printing
c..tolx    = convergence criteria, maximum allowd relative difference
c..tempf   = maximum temperature change allowed between time steps. 
c..          only used to form the next suggested timestep
c..func    = external routine, returns the diffusion equation coefficients 
c..          e(T), rho(T),sigma(T), s(T) and their derivatives with T.


c..output:
c..dt      = timestep taken
c..tnew    = temperature on the grid
c..flux    = value of the sigma(T) * dT/dt term
c..term1   = value of de(T)/dt term
c..term2   = value of 1/rho * d/dx (sigma(T) * dT/dt) term
c..dttemp  = suggestion for the next timestep based on tempf
c..itemp   = zone that is setting the value of dttemp
c..iter    = number of newton-raphson iterations to get the new solution


c..declare the pass
      external         func
      integer          ngrid,init,j1,jn,ncycle,itemp,ifail
      double precision x(ngrid),tnew(ngrid),dt,dtnext,
     1                 neu1l,xlpow,neu2l,neu1r,xrpow,neu2r,tolx,tempf,
     2                 flux(ngrid),term1(ngrid),term2(ngrid),dttemp 


c..local variables
      integer          nmax
      parameter        (nmax=800)
      integer          itmax,maxbak,i,j,iter,ibakup,irrmax
      parameter        (itmax=20, maxbak=5)
      double precision oper(nmax),told(nmax),
     1                 dold(nmax),dnew(nmax),dnewdt(nmax),
     2                 eold(nmax),enew(nmax),denewdt(nmax),
     3                 aold(nmax),anew(nmax),danewdt(nmax),
     4                 sold(nmax),snew(nmax),dsnewdt(nmax),
     5                 lold(nmax),lnew(nmax),dlnewdt(nmax),
     6                 flux0,dflux01,dflux00,dflux0m,fluxm,dfluxm1,
     7                 dfluxm0,dfluxmm,u0,du00,du01,up,dupm,dup0,dxp,
     8                 dxm,dx2,dto2,err,errmax,rat,xx,yy,si,rrmax,lim1,
     9                 lim2,lim3,tiny,alfa
      parameter        (tiny=1.0e-12, alfa = 1.0d0)




c..for the solution to the tridiagonal matrix
      double precision aa(nmax),bb(nmax),cc(nmax),rr(nmax),corr(nmax)



c..initialize the coefficients and store the solution
      if (init .eq. 1) then
       init    = 0
       if (ngrid .gt. nmax) stop 'ngrid > nmax in routine crni'
        call func(ngrid,tnew,
     1            dnew,dnewdt,
     2            enew,denewdt,
     3            anew,danewdt,
     4            snew,dsnewdt,
     5            lnew,dlnewdt)
       do i=1,ngrid
        told(i) = tnew(i)
        eold(i) = enew(i)
        dold(i) = dnew(i)
        aold(i) = anew(i)
        sold(i) = snew(i)
        lold(i) = lnew(i)
       enddo
      end if



c..backup loop; hopefully we do this loop once, that is, the newton-raphson 
c..loop convergences the first time 

      ifail = 0
      iter  = 0
      do ibakup=1,maxbak



c..copy old solution
c..and make a linear extrapolation guess for the new solution
       rat = alfa * dtnext/dt
       do i=1,ngrid
        told(i) = tnew(i)
        eold(i) = enew(i)
        dold(i) = dnew(i)
        aold(i) = anew(i)
        sold(i) = snew(i)       
        lold(i) = lnew(i)
        tnew(i) = tnew(i) + rat * (tnew(i) - told(i))
       enddo
       if (dtnext .ne. 0.0) dt = dtnext


c..this section computes the old operator
c..for the first grid point, compute flux, limit flux, store the operator
       i       = 1
       dxp     = x(i+1) - x(i)
       dxm     = dxp
       dx2     = dxp + dxm
       flux0   = aold(i) * (told(i+1) - told(i))/dxp
       u0      = told(i+1) - dx2*(neu1l*told(i)**xlpow + neu2l)
       fluxm   = aold(i) * (told(i) - u0)/dxm
       xx      = flux0/lold(i)
       lim1    = 1.0d0/(1.0d0 + abs(xx))
       flux0   = flux0 * lim1
       oper(i) = (flux0 - fluxm)/(dxm * dold(i)) + sold(i)
       dxm     = dxp
       fluxm   = flux0

c..same but for the middle of the grid
       do i=2,ngrid-1
        dxp     = x(i+1) - x(i)
        dx2     = dxp + dxm
        flux0   = aold(i) * (told(i+1) - told(i))/dxp
        xx      = flux0/lold(i)
        lim1    = 1.0d0/(1.0d0 + abs(xx))
        flux0   = flux0 * lim1
        oper(i) = (flux0 - fluxm)/(dxm * dold(i)) + sold(i)
        dxm     = dxp
        fluxm   = flux0
       enddo

c..same, but for the last grid point
       i       = ngrid
       dxp     = dxm
       dx2     = dxp + dxm
       up      = told(i-1) + dx2*(neu1r*told(i)**xrpow + neu2r)
       flux0   = aold(i) * (up - told(i))/dxp
       xx      = flux0/lold(i)
       lim1    = 1.0d0/(1.0d0 + abs(xx))
       flux0   = flux0 * lim1
       oper(i) = (flux0 - fluxm)/(dxm * dold(i)) + sold(i)
       dxm     = dxp
       fluxm   = flux0





c..start the newton-raphson iteration
       dto2 = 0.5d0 * dt
       do j=1,itmax

c..get the coefficients for the whole grid
         call func(ngrid,tnew,
     1             dnew,dnewdt,
     2             enew,denewdt,
     3             anew,danewdt,
     4             snew,dsnewdt,
     5             lnew,dlnewdt)


c..form the crank-nicholson tridiagonal correction matrix 
c..for the first grid point
        i       = 1
        dxp     = x(i+1) - x(i)
        dxm     = dxp
        dx2     = dxp + dxm

        flux0   = anew(i) * (tnew(i+1) - tnew(i))/dxp
        dflux01 = anew(i)/dxp
        dflux00 = (danewdt(i)*(tnew(i+1) - tnew(i)) - anew(i))/dxp
        dflux0m = 0.0d0

        u0      = tnew(i+1) - dx2*(neu1l*tnew(i)**xlpow + neu2l)
        du01    = 1.0d0
        du00    = -dx2*neu1l*xlpow*tnew(i)**(xlpow-1.0d0)

        fluxm   = anew(i) * (tnew(i) - u0)/dxm
        dfluxm1 = -anew(i)/dxm * du01
        dfluxm0 = (danewdt(i)*(tnew(i)-u0) + anew(i)*(1.0d0-du00))/dxm
        dfluxmm = 0.0d0

        lim1    = 1.0d0/(1.0d0 + abs(flux0/lnew(i)))
        flux0   = flux0 * lim1
        si      = sign(1.0d0,flux0)
        lim2    = lim1 * (1.0d0 - si * flux0/lnew(i))
        lim3    = si * (flux0/lnew(i))**2
        dflux01 = dflux01 * lim2
        dflux00 = dflux00 * lim2 + dlnewdt(i) * lim3


c..dirchlett on the first grid point. 
c..rr(i) is the difference equation.
c..aa(i) is derivative of rr with respect to temp in zone  i - 1 
c..bb(i) is derivative of rr with respect to temp in zone  i
c..cc(i) is derivative of rr with respect to temp in zone  i + 1
        if (j1 .eq. 0) then
         rr(i)    = tnew(i) - told(i)
         aa(i)    = 0.0d0
         bb(i)    = 1.0d0
         cc(i)    = 0.0d0
         flux(i)  = 0.0d0
         term1(i) = 0.0d0
         term2(i) = 0.0d0

c..or neumann on the first grid point
        else if (j1 .eq. 1) then
         flux(i)  = flux0
         term1(i) = (enew(i) - eold(i))/dto2
         term2(i) = (flux0 - fluxm)/(dxm * dnew(i))
         rr(i)    = enew(i) - eold(i) - dto2*(term2(i)+snew(i)+oper(i))
         aa(i)    = 0.0d0
         bb(i)    = denewdt(i) - dto2 *((dflux00-dfluxm0)/(dxm*dnew(i))
     1              - term2(i)/dnew(i)*dnewdt(i)
     2              + dsnewdt(i))
         cc(i)    = -dto2 * (dflux01 - dfluxm1)/(dxm * dnew(i))
        end if


c..swap fluxes, grid differences and get the first max difference
        dxm     = dxp
        fluxm   = flux0
        dfluxm1 = 0.0d0
        dfluxm0 = dflux01
        dfluxmm = dflux00
        rrmax   = abs(rr(i))
        irrmax  = i



c..for the middle of the grid
        do i=2,ngrid-1
         dxp     = x(i+1) - x(i)
         dx2     = dxp + dxm

         flux0   = anew(i) * (tnew(i+1) - tnew(i))/dxp
         dflux01 = anew(i)/dxp
         dflux00 = (danewdt(i)*(tnew(i+1) - tnew(i)) - anew(i))/dxp
         dflux0m = 0.0d0

c..limit the flux; wilson perscription
         lim1    = 1.0d0/(1.0d0 + abs(flux0/lnew(i)))
         flux0   = flux0 * lim1
         si      = sign(1.0d0,flux0)
         lim2    = lim1 * (1.0d0 - si * flux0/lnew(i))
         lim3    = si * (flux0/lnew(i))**2
         dflux01 = dflux01 * lim2
         dflux00 = dflux00 * lim2 + dlnewdt(i) * lim3


c..store come quatities
         flux(i)  = flux0
         term1(i) = (enew(i) - eold(i))/dto2
         term2(i) = (flux0 - fluxm)/(dxm * dnew(i))


c..rr(i) is the difference equation.
c..aa(i) is derivative of rr with respect to temp in zone  i - 1 
c..bb(i) is derivative of rr with respect to temp in zone  i
c..cc(i) is derivative of rr with respect to temp in zone  i + 1

         rr(i)    = enew(i) - eold(i) - dto2*(term2(i)+snew(i)+oper(i))
         aa(i)    = -dto2 * ((dflux0m - dfluxmm)/(dxm * dnew(i)))
         bb(i)    = denewdt(i) - dto2 *((dflux00-dfluxm0)/(dxm*dnew(i))
     1              - term2(i)/dnew(i)*dnewdt(i)
     2              + dsnewdt(i))
         cc(i)    = -dto2 * ((dflux01 - dfluxm1)/(dxm * dnew(i)))



c..push the forward spatial step, rotate the flux derivatives
         dxm     = dxp
         fluxm   = flux0
         dfluxm1 = 0.0d0
         dfluxm0 = dflux01
         dfluxmm = dflux00


c..get the maximum rr and where it occurs, then back for another grid point
         err = abs(rr(i))
         if (err .gt. rrmax) then
          rrmax  = err
          irrmax = i
         end if
        enddo



c..for the last grid point
        i       = ngrid
        dxp     = dxm
        dx2     = dxp + dxm
        up      = tnew(i-1) + dx2*(neu1r*tnew(i)**xrpow + neu2r)
        dup0    = dx2*neu1r*xrpow*tnew(i)**(xrpow-1.0d0)
        dupm    = 1.0d0

        flux0   = anew(i) * (up - tnew(i))/dxp
        dflux01 = 0.0d0
        dflux00 = (danewdt(i)*(up-tnew(i)) + anew(i)*(dup0-1.0d0))/dxp
        dflux0m = anew(i)/dxp * dupm

        lim1    = 1.0d0/(1.0d0 + abs(flux0/lnew(i)))
        flux0   = flux0 * lim1
        si      = sign(1.0d0,flux0)
        lim2    = lim1 * (1.0d0 - si * flux0/lnew(i))
        lim3    = si * (flux0/lnew(i))**2
        dflux01 = dflux01 * lim2
        dflux00 = dflux00 * lim2 + dlnewdt(i) * lim3

c..aa(i) is derivative of rr with respect to temp in zone  i - 1 
c..bb(i) is derivative of rr with respect to temp in zone  i
c..cc(i) is derivative of rr with respect to temp in zone  i + 1
c..dirichlet on the outer boundary
        if (jn .eq. 0) then
         rr(i)     = tnew(i) - told(i)
         aa(i)     = 0.0d0
         bb(i)     = 1.0d0
         cc(i)     = 0.0d0
         flux(i)   = 0.0d0
         term1(i)  = 0.0d0
         term2(i)  = 0.0d0

c..or neumann on the outer boundary
        else if (jn .eq. 1) then
         flux(i)   = flux0
         term1(i)  = (enew(i) - eold(i))/dto2
         term2(i)  = (flux0-fluxm)/(dxm*dnew(i))
         rr(i)     = enew(i)-eold(i) - dto2*(term2(i)+snew(i)+oper(i))
         aa(i)     = -dto2 * ((dflux0m - dfluxmm)/(dxm * dnew(i)))
         bb(i)     = denewdt(i) - dto2*((dflux00-dfluxm0)/(dxm*dnew(i))
     1              - term2(i)/dnew(i)*dnewdt(i)
     2               + dsnewdt(i))
         cc(i)     = 0.0d0
        end if


c..the last maximum difference
        err = abs(rr(i))
        if (err .gt. rrmax) then
         rrmax  = err
         irrmax = i
        end if




c..flip the correction matrix; the adjustments are in the corr array
        call tridag(aa,bb,cc,rr,corr,ngrid)



c..record the maximum temperature correction
        errmax = 0.0d0
        do i = 1,ngrid
         if (abs(tnew(i)) .gt. tiny) then
          err = abs(corr(i)/(tnew(i) + tiny))
          if (err .gt. errmax) then
           errmax  = err
          end if
         end if
        enddo


c..apply the corrections to the solution
        do i=1,ngrid
         tnew(i) = tnew(i) - corr(i)
        enddo
        iter = iter + 1


c..this is the normal exit point if we have converged
c..form a suggestion for the next timestep dttemp based on tempf

        if (errmax .lt. tolx) then
         xx = 1.0d99
         do i=1,ngrid
          yy = abs(tnew(i)-told(i)) 
          if (yy*xx .gt. told(i)) then
           xx    = told(i)/yy
           itemp = i
          end if
         enddo
         dttemp = xx*dt*tempf
         return  
        end if

c..end of newton-raphson loop
       enddo



c..land here if newton-raphson failed 
c..cut the timestep in half and reset the new solution

       dtnext = 0.5d0 * dtnext
       if (ibakup .ne. maxbak) then
        do i=1,ngrid
         tnew(i) = told(i)
         enew(i) = eold(i)
         dnew(i) = dold(i)
         anew(i) = aold(i)
         snew(i) = sold(i)       
         lnew(i) = lold(i)
        enddo
       end if


c..end of backup loop
      enddo
      write(6,*) 'error: exceeded newton loop at cycle',ncycle
      write(6,*) 'setting ifail to 1 and returning'
      ifail = 1
      return
      end





      subroutine tridag(a,b,c,r,u,n) 
      implicit none
      save

c..solves for a vector u of length n in the tridiagonal linear system 
c.. a_i u_(i-1) + b_i u_i + c_i u_(i+1) = r_i 
c..input are the a, b, c, and r vectors. they are not modified on output.

c..in its clearest incarnation, this algorithm uses two storage arrays,
c..say p and q. here, the solution vector u is used for q, cutting 
c..the extra storage down to one array. 

c..declare 
      integer          n,nmax,i 
      parameter        (nmax=800) 
      double precision a(n),b(n),c(n),r(n),u(n),p(nmax),bet,xx 


      if (b(1) .eq. 0.0)  stop 'eliminate u2 trivially' 
      xx = 1.0d0/b(1) 
      p(1) = c(1)*xx
      u(1) = r(1)*xx
      do i=2,n 
       bet = b(i) - a(i) * p(i-1) 
       if (bet .eq. 0.0) stop 'singular in tridag3' 
       xx = 1.0d0/bet
       p(i) = c(i)*xx
       u(i) = (r(i) - a(i) * u(i-1) )*xx
      enddo
      do i=n-1,1,-1 
       u(i) = u(i) - p(i) * u(i+1) 
      enddo
      return 
      end 




      double precision function aroot(a)
      implicit none
      save

c..function for the problem of smith page 35
      double precision a

      aroot = a*tan(a) - 0.5d0
      return
      end




      double precision function zbrent(func,x1,x2,tol,niter)   
      implicit none
      save

c..using brent's method this routine finds the root of a function func  
c..between the limits x1 and x2. the root is when accuracy is less than tol. 
c.. 
c..note: eps the the machine floating point precision 
c.. 
c..declare 
      external          func 
      integer           niter,itmax,iter   
      parameter         (itmax = 100)   
      double precision  func,x1,x2,tol,a,b,c,d,e,fa, 
     1                  fb,fc,xm,tol1,p,q,r,s,eps    
      parameter         (eps=3.0d-15)   

c..initialize 
      niter = 0
      a     = x1 
      b     = x2 
      fa    = func(a)   
      fb    = func(b)   
      if ( (fa .gt. 0.0  .and. fb .gt. 0.0)  .or. 
     1     (fa .lt. 0.0  .and. fb .lt. 0.0)       ) then 
       write(6,100) x1,fa,x2,fb 
100    format(1x,' x1=',1pe11.3,' f(x1)=',1pe11.3,/, 
     1        1x,' x2=',1pe11.3,' f(x2)=',1pe11.3) 
       stop 'root not bracketed in routine zbrent'    
      end if 
      c  = b 
      fc = fb    

c..rename a,b,c and adjusting bound interval d   
      do iter =1,itmax  
       niter = niter + 1  
       if ( (fb .gt. 0.0  .and. fc .gt. 0.0)  .or. 
     1      (fb .lt. 0.0  .and. fc .lt. 0.0)      ) then 
        c  = a    
        fc = fa  
        d  = b-a  
        e  = d    
       end if    
       if (abs(fc) .lt. abs(fb)) then    
        a  = b    
        b  = c    
        c  = a    
        fa = fb  
        fb = fc  
        fc = fa  
       end if    
       tol1 = 2.0d0 * eps * abs(b) + 0.5d0 * tol 
       xm   = 0.5d0 * (c-b)  

c..convergence check 
       if (abs(xm) .le. tol1 .or. fb .eq. 0.0) then  
        zbrent = b   
        return   
       end if    

c..attempt quadratic interpolation   
       if (abs(e) .ge. tol1 .and. abs(fa) .gt. abs(fb)) then 
        s = fb/fa    
        if (a .eq. c) then   
         p = 2.0d0 * xm * s    
         q = 1.0d0 - s  
        else 
         q = fa/fc   
         r = fb/fc   
         p = s * (2.0d0 * xm * q *(q-r) - (b-a)*(r - 1.0d0))   
         q = (q - 1.0d0) * (r - 1.0d0) * (s - 1.0d0) 
        end if   

c..check if in bounds    
        if (p .gt. 0.0) q = -q    
        p = abs(p)   

c..accept interpolation  
        if (2.0d0*p .lt. min(3.0d0*xm*q - abs(tol1*q),abs(e*q))) then    
         e = d   
         d = p/q 

c..or bisect
        else 
         d = xm  
         e = d   
        end if   

c..bounds decreasing to slowly use bisection 
       else  
        d = xm   
        e = d    
       end if    

c..move best guess to a  
       a  = b 
       fa = fb   
       if (abs(d) .gt. tol1) then    
        b = b + d    
       else  
        b = b + sign(tol1,xm)    
       end if    
       fb = func(b)  
      enddo
      stop 'too many iterations in routine zbrent'   
      end    





      subroutine zbrak(fx,x1,x2,n,xb1,xb2,nb)    
      implicit none
      save

c..given a function fx defined on the interval x1 to x2, subdivide the  
c..interval into n equally spaced segments, and search for zero crossings  
c..of the function. nb is input as the maximum number of roots sought and  
c..is reset to the number of bracketing pairs xb1,xb2 that are found.    
c..  
c..declare 
      external           fx 
      integer            n,nb,nbb,i  
      double precision   fx,x1,x2,xb1(nb),xb2(nb),x,dx,fp,fc   

c..initialize and determine the appropriate mesh spacing  
      nbb = 0 
      x   = x1   
      dx  = (x2-x1)/n    
      fp  = fx(x) 

c..loop over all intervals 
      do i=1,n    
       x  = x + dx   
       fc = fx(x)    

c..if a sign change occurs then record the values for the bounds 
       if (fc*fp .le. 0.0) then  
        nbb      = nbb + 1 
        xb1(nbb) = x - dx   
        xb2(nbb) = x  
        if (nbb .eq. nb) goto 10 
       end if    
       fp = fc   
      enddo
10    nb = nbb 
      return 
      end