program   transport
	include 'transin'

c  1d simulation program for linear transport using finite differences 
c  1 - the upwind scheme
c  2 - the leapfrog scheme
c  3 - the Lax-Wendroff scheme
c  4 - Crank-Nicholson scheme
c  output of program 
c    sim1.bin - unformated (first nx,x, then each output: time, f)     
c    sim1.asci - formated (first all parameters, then: time, f) 
      open (14,file='sim1.bin',form='unformatted')
      open (16,file='sim1.asci')
c   f,fold,foold are normalized fields for integration   
c   read input parameters
      call pread
c   determine grid:
      call grid
c   calculate runtime params:
      call params
      write(14) nx
      write(14) x
c   write some parameters:      
      call pwrite
      write(14) c
      write(14) alpha
c   set initial conditions:      
       if (init .eq. 1) call initcon1
       if (init .eq. 2) call initcon2
       if (init .eq. 3) call initcon3
       if (init .eq. 4) call initcon4
c   integration: 
c    intftcs - 1d heat conduction with ftcs 
c              (forward time centered space)
       nt = nt+1
       time = time + dt
       if (method .eq. 1) call upwind
       if (method .eq. 2) call upleap
       if (method .eq. 3) call lax
       if (method .eq. 3) call laxwend
       if (method .eq. 4 .and. que .eq. 0.) call crankni
       if (method .eq. 4 .and. que .ne. 0.) call cnupwind
       call bdcon
       if (mod(nt,nout) .eq. 0) call out
       do while( (nt .le. ntmax) .and. (time .le. tmax) )
         nt = nt+1
         time = time + dt
         if (method .eq. 1) call upwind
         if (method .eq. 2) call leap
         if (method .eq. 3) call lax
         if (method .eq. 3) call laxwend
         if (method .eq. 4 .and. que .eq. 0.) call crankni
         if (method .eq. 4 .and. que .ne. 0.) call cnupwind
         call bdcon
c   write output f:      
        if (mod(nt,nout) .eq. 0) call out
       end do
      call out
c   determine and write exact solution to order nterm:  
      call analyt
      close(14) 
      close(16) 
      stop
      end
c ------------ : ------------
      subroutine pread
      include 'transin'
c   input data:
c   parameters explained in input file sim1.dat
       open(11,file='trans1.dat')          
       read(11,1) method
       read(11,1) ntmax
       read(11,3) tmax
       read(11,3) xmin,xmax
       read(11,2) alpha
       read(11,3) u,c,dt
       read(11,1) cordt
       read(11,3) delta,que
       read(11,1) nout
       read(11,3) f0
       read(11,1) init
       read(11,3) lambda
       read(11,1) parbound
       read(11,3) fb(1),fb(2)
       read(11,1) nterm
       close(11) 
       
       if (method.eq.4 .and. que.ne.0. .and. parbound.eq.1) then
         parbound=2
         write(6,*) '!!!! NO PERIODIC BOUNDARY CONDITION IMPLEMENTED',
     +            '      FOR THE CRANK-NICHOLSON SCHEME!!!! '
         write(6,*) '!!!! BOUNDARY CONDITION SET TO DIRICHLET!!!! '
       endif
    1  format(i5)
    2  format(e10.2)
    3  format(f8.2)
      return
      end
c ------------ : ------------
      subroutine grid
      include 'transin'
       integer i
       xsize=xmax-xmin
       dx = (xmax-xmin)/(nx-3.)
       do 10 i = 1,nx
   10    x(i)  = xmin+dx*(i-2.)
      return
      end
c ------------ : ------------
      subroutine params
      include 'transin'
c   compute parameters for schemes
       if (cordt .eq. 0) dt = c*dx/u
       if (cordt .eq. 1) c = u*dt/dx
       s = alpha*dt/dx/dx
       rcell=0.
       if (alpha .ne. 0.) rcell = u*dx/alpha
       nt = 0
       time = 0.0      
       
c       delta=1./6.+c**2/12
c        que=0.5+0.25*c*c
       
       star = 0.
       if (method .eq. 1)  then
         s1 = s+c
         s2 = 1. - 2.*s - c
       endif
       if (method .eq. 2)  then
         star = 1.+2.*s
         s1   = (c+2.*s)/star
         s2   = (1.-2.*s)/star
         s3   = (-c+2.*s)/star
       endif
       if (method .eq. 3)  then
         star = s + 0.5*c*c
         s0   = -que*c/3.
         s1   = 0.5*c + star + que*c
         s2   = 1. - 2.*star - que*c
         s3   = -0.5*c + star + que*c/3.
       endif
       if (method .eq. 4)  then
         a0 = que*c/6.
         a1 = delta - 0.25*c - 0.5*s - 0.5*que*c
         a2 = 1. - 2.*delta + s + 0.5*que*c
         a3 = delta + 0.25*c - 0.5*s - 0.5*que*c/3.
         s0 = -que*c/6.
         s1 = delta + 0.25*c + 0.5*s + 0.5*que*c
         s2 = 1. - 2.*delta - s - 0.5*que*c
         s3 = delta - 0.25*c + 0.5*s + 0.5*que*c/3.
       endif
      write(6,*) delta
      return
      end
c ------------ : ------------
      subroutine pwrite
      include 'transin'
       if (method .eq. 1) write(16,*) 'Upwind scheme (explicit)'
       if (method .eq. 2) write(16,*) 'Leapfrog scheme (explicit)'
       if (method .eq. 3) write(16,*) 'Lax-Wendroff scheme (explicit)'
       if (method .eq. 4) then 
         write(16,*) 'Crank-Nicholson scheme (implicit)'
       endif
       write(16,2) nx,ntmax,tmax, xmin, xmax
       write(16,3) dt,dx,nout
       write(16,4) u,alpha,c,s,rcell
       if (method .eq. 3 .or. method .eq. 4) write(16,5) delta, que
       write(16,6) init,lambda,f0
       if (parbound .eq. 1) write(16,*) 'Periodic boundary conditions'
       if (parbound .eq. 2) write(16,*) 
     +              'dirichlet boundary conditions with:'
       if (parbound .eq. 3) write(16,*) 
     +              'neumann boundary conditions with derivatives:'
       if (parbound .ne. 1) write(16,7) fb(1),fb(2)
       write(16,8) nterm
       write(16,*) 'Stability:'
       if (method .eq. 1) write(16,9) c+2.*s, rcell, 2./(1.-c)
       if (method .eq. 2) write(16,10) c, c**2
       if (method .eq. 3) write(16,11) c**2, 2*star, rcell
       if (method .eq. 4) write(16,12) rcell
       
       
       if (method .eq. 1) write(6,*) 'Upwind scheme (explicit)'
       if (method .eq. 2) write(6,*) 'Leapfrog scheme (explicit)'
       if (method .eq. 3) write(6,*) 'Lax-Wendroff scheme (explicit)'
       if (method .eq. 4) then 
         write(6,*) 'Crank-Nicholson scheme (implicit)'
       endif
       write(6,2) nx,ntmax,tmax, xmin, xmax
       write(6,3) dt,dx,nout
       write(6,4) u,alpha,c,s,rcell
       if (method .eq. 3 .or. method .eq. 4) write(6,5) delta, que
       write(6,6) init,lambda,f0
       if (parbound .eq. 1) write(6,*) ' Periodic boundary conditions'
       if (parbound .eq. 2) write(6,*) 
     +              'dirichlet boundary conditions with:'
       if (parbound .eq. 3) write(6,*) 
     +              'neumann boundary conditions with derivatives:'
       if (parbound .ne. 1) write(6,7) fb(1),fb(2)
       write(6,8) nterm
       if (method .eq. 1) write(6,9) c+2.*s, rcell, 2./(1.-c)
       if (method .eq. 2) write(6,10) c, c**2
       if (method .eq. 3) write(6,11) c**2, 2*star, rcell
       if (method .eq. 4) write(6,12) rcell
    2  format(1x,' nx =',i4,'   ntmax =',i5,'     tmax =',f8.2,
     +       '  xmin =',f5.1,'  xmax =',f5.1)
    3  format(1x,' dt =',f8.4,'    dx =',f6.4,' nout =',i5)
    4  format(1x,' u  =',f6.2,'  alpha =',e10.3,
     +           ' c =',f6.4,' s =',f6.4,'  rcell =',e10.3)
    5  format(1x,' delta=',f8.5,'  q(upwind) =',f8.5)
    6  format(1x,' init(cond) =',i3,'  lambda =',f8.4,
     +           ' f0(norm) =',f6.1)
    7  format(1x,' f_bound_min =',f7.3,' f_bound_max =',f7.3)
    8  format(1x,' Terms in exact solution (nterm) =',i4)

    9  format(1x,' Stability: c+2.*s =',f7.3,' <=1; ',/,
     +           '  Accuracy: rcell =',e10.3,',<= 2./(1.-c) =',e10.3)
   10  format(1x,' Stability: c =',f7.3,',<=1; ',
     +           '  Accuracy: c**2 =',e10.3,',<= 1')
   11  format(1x,' Stability: c**2 =',f7.3,',<= ',f8.3,' =2*s_star <=1',
     +         /,'  Oscillations: rcell =',e10.3,',<= 2')
   12  format(1x,'  Oscillations: rcell =',e10.3,',<= 2')
      return
      end
c ------------ : ------------
      subroutine initcon1
      include 'transin'
       integer i
       real    lhalf
       lhalf=0.5*lambda
       do 10 i = 1,nx-1
   10     f(i) = 0.
       do 20 i = 2,nx-1
   20    if (x(i) .ge. -lhalf .and. x(i).le. lhalf)
     +            f(i) = cos(pi*x(i)/lambda)
       do 25 i=1,nterm
         kf(i)= (pi/2.+(i-1.)*pi)/xmax
         arg1=pi/lambda+kf(i)
         arg2=pi/lambda-kf(i)
         if (arg2 .ne. 0) cf(i)=1./xmax*( 
     +            (sin(arg1*lhalf)/arg1)+(sin(arg2*lhalf)/arg2) )
   25  continue
       do 30 i = 2,nx-1
        fexact(i) = 0
        do 30 l = 1,nterm
        fexact(i)=fexact(i)+cf(l)*cos(kf(l)*x(i))
   30  continue
       call bdcon
       do 40 i = 2,nx-1
        finit(i)=f(i)
   40  continue
       call out
      return
      end
c ------------ : ------------
      subroutine initcon2
      include 'transin'
       integer  i,nwave
       nwave=xsize/(2*lambda)
       kwave=2.*pi*nwave/xsize
       do 10 i = 2,nx-1
         f(i) = sin(kwave*x(i))
         fexact(i) = f(i)
   10  continue
       call bdcon
       do 40 i = 2,nx-1
        finit(i)=f(i)
   40  continue
       call out
      return
      end
c ------------ : ------------
      subroutine initcon3
      include 'transin'
       integer i
       real    lhalf
       lhalf=0.5*lambda
       do 10 i = 1,nx-1
   10     f(i) = 0.
       do 20 i = 2,nx-1
   20    if (x(i) .le. 0.0)   f(i) = 1.
       call bdcon
       do 40 i = 2,nx-1
        finit(i)=f(i)
   40  continue
       call out
      return
      end
c ------------ : ------------
      subroutine initcon4
      include 'transin'
       integer i
       real    lhalf
       lhalf=0.5*lambda
       do 10 i = 1,nx-1
   10     f(i) = 0.
       do 20 i = 2,nx-1
   20    if (x(i) .ge. -lhalf .and. x(i).le. lhalf)
     +            f(i) = 1.
       do 25 i=1,nterm
         kf(i)= (pi/2.+(i-1.)*pi)/xmax
         cf(i)=2./xmax*sin( kf(i)*lhalf )/kf(i)
   25  continue
       do 30 i = 2,nx-1
        fexact(i) = 0
        do 30 l = 1,nterm
        fexact(i)=fexact(i)+cf(l)*cos(kf(l)*x(i))
   30  continue
       call bdcon
       do 40 i = 2,nx-1
        finit(i)=f(i)
   40  continue
       call out
      return
      end
c ------------ : ------------
      subroutine bdcon
      include 'transin'
       if (parbound .eq. 1) then
        f(1)  = f(nx-2)
        f(nx) = f(3)
       endif
       if (parbound .eq. 2) then
        f(2)  = fb(1)
        f(nx-1) = fb(2)
        f(1)  = fb(1)
        f(nx) = fb(2)
       endif
       if (parbound .eq. 3) then
        f(1)  = f(3)-2.*dx*fb(1)
        f(nx) = f(nx-2)+2.*dx*fb(2)
       endif
      return
      end
c ------------ : ------------
      subroutine upwind
      include 'transin'
      integer  i
       dtdx=dt/dx
       s1=1.-2.*s
       do 10 i = 1,nx
        fold(i) = f(i)
        g(i) = 0.5*f(i)*f(i)
   10  continue
       do 20 i = 2,nx-1
   20   f(i) = s1*fold(i) + s*(fold(i+1)+fold(i-1)) 
     +               - dtdx*(g(i)-g(i-1))
      return
      end
c ------------ : ------------
      subroutine upleap
      include 'transin'
      integer  i
       dtdx=dt/dx
       s1=1.-2.*s
       do 10 i = 1,nx
        fold(i) = f(i)
        g(i) = 0.5*f(i)*f(i)
   10  continue
       do 20 i = 2,nx-1
   20   f(i) = s1*fold(i) + s*(fold(i+1)+fold(i-1)) 
     +                - dtdx*(g(i)-g(i-1))
      return
      end
c ------------ : ------------
      subroutine leap
      include 'transin'
      integer  i
         dtdx=dt/dx
         star = 1.+2.*s
         s1   = (1.-2.*s)/star
         s2   = 2.*s/star
         s3   = dtdx/star
       do 10 i = 1,nx
   10   foold(i) = fold(i)
       do 15 i = 1,nx
        fold(i) = f(i)
        g(i) = 0.5*f(i)*f(i)
   15  continue
       do 20 i = 2,nx-1
   20   f(i) = s1*foold(i) + s2*(fold(i+1)+fold(i-1))
     +                  - s3*(g(i+1)-g(i-1))
      return
      end
c ------------ : ------------
      subroutine leapvis
      include 'transin'
      integer  i
      real     fhelp(nx),fvis(nx)
         vvv = 0.2
         dtdx=dt/dx
         star = 1.+2.*s
         s1   = (1.-2.*s)/star
         s2   = 2.*s/star
         s3   = dtdx/star
       do 10 i = 1,nx
   10   foold(i) = fold(i)
       do 15 i = 1,nx
        fold(i) = f(i)
        g(i) = 0.5*f(i)*f(i)
   15  continue
       do 20 i = 2,nx-1
   20   f(i) = s1*foold(i) + s2*(fold(i+1)+fold(i-1))
     +                  - s3*(g(i+1)-g(i-1))
       do 30 i = 2,nx-1
        fvis(i) =  vvv*0.5*(f(i+1)+f(i-1))
   30   fhelp(i) = (f(i+1)-f(i))*(f(i)-f(i-1))
       do 40 i = 4,nx-3
   40   if ( fhelp(i) .le. 0 .and. 
     +      ( fhelp(i+1) .le. 0 .or. fhelp(i-1) .le. 0  
     +        .or. fhelp(i+2) .le. 0 .or. fhelp(i-2) .le. 0) )
     +     f(i) = (1.0-vvv)*f(i) + fvis(i)
      return
      end
c ------------ : ------------
      subroutine lax
      include 'transin'
      integer  i
      real     fhmin,fhmax
         dtdx=dt/dx
       do 10 i = 1,nx
        g(i) = 0.5*f(i)*f(i)
        flw(i)=0.
   10  continue
        do 30 i = 2,nx-2
         flw(i) = 0.5*( f(i)+f(i+1) ) - 0.5*dtdx*( g(i+1)-g(i) )
     +          + 0.25*s*( (f(i-1)+f(i+1)-2.*f(i))
     +                    + (f(i)+f(i+2)-2.*f(i+1)) )
   30  continue
        nx1=nx-1
        nx2=nx-2
        nx3=nx-3
        if (parbound .eq. 1) then
         fhmin = f(nx3)
         fhmax = f(4)
        elseif (parbound .eq. 2) then
         fhmin = fb(1)
         fhmax = fb(2)
        elseif (parbound .eq. 3) then
         fhmin = f(2)-2.*dx*fb(1)
         fhmax = f(nx1)+2.*dx*fb(2)
        endif
         flw(1) = 0.5*(f(1)+f(2)) - 0.5*dtdx*(g(2)-g(1))
     +          + 0.25*s*( (fhmin+f(2)-2.*f(1))
     +                    + (f(1)+f(3)-2.*f(2)) )
         flw(nx1) = 0.5*(f(nx1)+f(nx)) - 0.5*dtdx*(g(nx)-g(nx1))
     +          + 0.25*s*( (f(nx2)+f(nx)-2.*f(nx1))
     +                    + (f(nx1)+fhmax-2.*f(nx)) )
      return
      end
c ------------ : ------------
      subroutine laxwend
      include 'transin'
      integer  i
         dtdx=dt/dx
       do 10 i = 1,nx
        fold(i) = f(i)
   10  continue
       do 20 i = 1,nx-1
        g(i) = 0.5*flw(i)*flw(i)
   20  continue
        do 30 i = 2,nx-1
         f(i) = fold(i) - dtdx*(g(i)-g(i-1)) 
     +          + s*(fold(i+1)+fold(i-1)-2.*fold(i))
   30   continue
      return
      end
c ------------ : ------------
      subroutine crankni
      include 'transin'
      integer  j
      real     a(5,nx),r(nx),td(nx)
       do 10 j = 2,nx-1
        jm=j-1
        a(1,jm)=0.
        a(2,jm)=a1
        a(3,jm)=a2
        a(4,jm)=a3
        r(jm) = s1*f(jm)+s2*f(j)+s3*f(j+1)
   10  continue
       r(1)=r(1)-a(2,1)*f(1)
       a(2,1)=0.
       a(4,nx-2)=0.
       
       call banfac(a,nx-2,1)
       call bansol(r,td,a,nx-2,1)

       do 30 i = 2,nx-1
   30   f(i) = td(i-1)
      return
      end
c ------------ : ------------
      subroutine cnupwind
      include 'transin'
      integer  j
      real     a(5,nx),r(nx),td(nx)
       do 10 j = 2,nx-1
        jm=j-1
        a(1,jm)=a0
        a(2,jm)=a1
        a(3,jm)=a2
        a(4,jm)=a3
        r(jm) = s1*f(jm)+s2*f(j)+s3*f(j+1)
   10  continue
       r(1)=r(1)-a(2,1)*f(1)
       if (que .ne. 0.) then
        do 20 j = 2,nx-2
   20     r(j) = r(j)+s0*f(j-1)
        r(1)=r(1)-a(1,1)*f(1)
        r(2)=r(2)-a(1,2)*f(1)
c  reduce to tridiagonal:   
        do 25 j = 3,nx-2
         jm = j-1
         dum=a(1,j)/a(2,jm)
         a(2,j) = a(2,j) - a(3,jm)*dum
         a(3,j) = a(3,j) - a(4,jm)*dum
         a(1,j) = 0.
         r(j) = r(j) - r(jm)*dum
   25   continue
        a(1,1) = 0.
        a(1,2) = 0.
       endif
       a(2,1)=0.
       a(4,nx-2)=0.
       
       call banfac(a,nx-2,1)
       call bansol(r,td,a,nx-2,1)

       do 30 i = 2,nx-1
   30   f(i) = td(i-1)
      return
      end
c ------------ : ------------
      subroutine out
      include 'transin'
      integer  i
      real     fout(nx),fexout(nx)
       do 10 i = 1,nx
        fexout(i)=f0*fexact(i)
   10   fout(i) = f0*f(i)
       write(16,1) time, (fout(i),i=1,nx)
       call analyt
       write(14) time,rms,fout,fexout
       write(16,*) 'time step:', nt
    1  format(' time = ',f5.0,'  T =', 11f7.2)
      return
      end
c ------------ : ------------
      subroutine analyt
      include 'transin'
      integer  i
      real     sum,xst,piquart
       piquart = sqrt(pi/4.)
       xst = -xmax
       do 10 i = 1,nx
        fexact(i)=0.
   10  continue
        
       do 20 i = 2,nx-1
        xb = -x(i)
        xa = x(i)-time
        pp = sqrt(pi*time*alpha)
        xa = xa*piquart/pp
        xb = xb*piquart/pp
        call errfc(xa,exa)
        call errfc(xb,exb)
        suma = pp*exa
        sumb = pp*exb
        dum = 0.5*(x(i)-0.5*time)/alpha
        if(dum .gt.  20.) dum =  20.
        if(dum .lt. -20.) dum = -20.
        sumc=exp(dum)
        fexact(i) = suma/(suma + sumc*sumb)
   20  continue

       sum=0.
       do 40 i = 2,nx-1
  40    sum = sum + ( fexact(i)-f0*f(i) )**2
       if (method .eq. 2) then 
         do 50 i = 2-mod(nt,2),nx,2
   50        sum = sum + ( fexact(i)-f0*f(i) )**2
       endif
c   rms is the rms error
       rms = sqrt(sum/nx)
       write(16,1) nterm
       write(16,2) time, (fexact(i),i=1,nx)
       write(16,3) rms
c       write(6,3) rms
    1  format(1x,'No of terms in exact solution (nterm) =',i5)
    2  format(/,' time = ',f5.0,'  T =', 11f7.2)
    3  format(1x,' rms dif = ', e11.4,//)
      return
      end
c ------------ : ------------
      subroutine errfc(xin,erc)
      include 'transin'
       b = abs(xin)
       if (b .lt. 4) goto 1
       d = 1.0
       goto 2
    1  c = exp(-b*b)
       t = 1./(1. + 0.3275911*b)
       d = 0.254829592*t - 0.284496736*t*t + 1.421413741*t*t*t
     +     - 1.453152027*t*t*t*t + 1.061405429*t*t*t*t
       d = 1. - d
    2 if (xin .lt. 0.) d = -d
      erc = 1.-d
      return
      end
c ------------ : ------------
      subroutine banfac(b,n,int)
      include 'transin'
c
c  factorizes band matrix arising from linear or quadratic elements
c  into l.u
c  
      real     b(5,nx)
      
c
c  int = 1, linear elements = tridiagonal system
c  int = 2, qudratic elements = pentadiagonal system
c  assumes that first equation formed at midside node
c  
      if (int .eq. 1) then
       np = n-1
       do 1 j = 1,np
        jp = j+1
        b(2,jp) = b(2,jp)/b(3,j)
        b(3,jp) = b(3,jp) - b(2,jp)*b(4,j)
    1  continue
       return
      else
       nh = n/2
       do 5 i = 1,2
        js = 3-i
        do 4 j = js,nh
         ja = 2*(j-1)
c  i = 1, first pass, reduction to tridiagonal
c  i = 2, second pass, reduction to upper triangular
         if (i .eq. 1) then
          jb = ja+2
          b(1,jb) = b(1,jb)/b(2,jb-1)
          b(2,jb) = b(2,jb) - b(1,jb)*b(3,jb-1)
          b(3,jb) = b(3,jb) - b(1,jb)*b(4,jb-1)
         else
          jb = ja+3
          b(2,jb-1) = b(2,jb-1)/b(3,jb-2)
          b(3,jb-1) = b(3,jb-1) - b(2,jb-1)*b(4,jb-2)
          if (jb .le. n) then
           b(2,jb) = b(2,jb)/b(3,jb-1)
           b(3,jb) = b(3,jb) - b(2,jb)*b(4,jb-1)
           b(4,jb) = b(4,jb) - b(2,jb)*b(5,jb-1)
          endif
         endif
    4   continue
    5  continue
      endif
      return
      end
c ------------ : ------------
      subroutine bansol(r,xx,b,n,int)
      include 'transin'
c
c  uses l.u factorization to solve for x, given r
c  
      real     r(nx),xx(nx),b(5,nx)

c
c  int = 1, tridiagonal system
c  int = 2, pentadiagonal system
c   assumes that first equation formed at midside node
c   ibc=0 if last equation formed at midside node
c   ibc=1 if last equation formed at corner node
c  
      if (int .eq. 1) then
       np = n-1
       do 1 j = 1,np
        jp = j+1
        r(jp) = r(jp) - b(2,jp)*r(j)
    1  continue
       xx(n)=r(n)/b(3,n)
       do 2 j = 1,np
        ja = n - j
        xx(ja) = ( r(ja) - b(4,ja)*xx(ja+1) )/b(3,ja)
    2  continue
       return
      else
    3  ibc = 1
       nh = n/2
       if (2*nh .eq. n) r(n+1) = 0.
       if (2*nh .eq. n) b(2,n+1) = 0.
       do 6 i = 1,2
        do 5 j = 1,nh
         ja = 2*j
         do 4 k = 1,i
          jb = ja -1 + k
          r(jb) = r(jb) - b(i,jb)*r(jb-1)
    4    continue
    5   continue
    6  continue
       nen = nh -ibc
       xx(n) = r(n)/b(3,n)
       if ( ibc .eq. 1) xx(n-1) = ( r(n-1) - b(4,n-1)*xx(n) )/b(3,n-1)
       do 7 j = 1,nen
        ja = n -2*j + 1 - ibc
        xx(ja) = ( r(ja) - b(4,ja)*xx(ja+1) - b(5,ja)*xx(ja+2) )/b(3,ja)
        xx(ja-1) = ( r(ja-1) - b(4,ja-1)*xx(ja) )/b(3,ja-1)
    7  continue
      endif
      return
      end