program DERIV                                                            
c                                                                               
c this simple program computes the first derivative of a function  
c using the second and fourth order accurate schemes. the following
c driver estimates the derivatives of  p(z) = po*exp(-a*z)                                         
c                                                                               
      parameter(l=10)                                                          
      real z(l),p(l),p2(l),p4(l),anal(l)
c                                       
      data z /0000.,1000.,2000.,3000.,4000.,                                       
     &         5000.,6000.,7000.,8000.,9000./                                    
      data a,p0,dz/0.000125062,1000.,1000./
c
c  initialize the work arrays
c
      do 2100 k = 1, l
         p2(k)  = 0.
         p4(k)  = 0.
2100  continue
c                                  
      do 2102 k = 1, l
c                                                             
c  construct the function p(z)
c                                                 
      p(k)      =  p0*exp(-a*z(k))
c                                                   
c  compute the analytical derivative
c                                           
      anal(k)   = -a*p0*exp(-a*z(k))
c                                             
 2102 continue 
c                                                                
c  compute the second order estimate                                           
c
      call DDX2 (p,p2,l,dz,1)                                                   
c
c  compute the fourth order estimate                                           
c
      call DDX4 (p,p4,l,dz)
c
c  write output . The 4th order derivative is omitted at the first 
c  and last 2 points  since it is not defined as 4th order.
c                                                     
      write (6,1000)                                                           
      write (6,1001)
c                                                           
      do 2104 k = 1, l 
         if (k.le.2.or.k.ge.(l-1)) then                                                            
      write(6,1002) anal(k),   p2(k)
         else                                   
      write(6,1003) anal(k),   p2(k),  p4(k)
         endif                                   
 2104 continue
c                                                                 
1000  format (5x,'analytical',15x,'second order',13x,'fourth order')           
1001  format (5x,' solution ',15x,'  estimate  ',13x,'  estimate  ',/)         
1002  format (5x,f10.6,15x,f10.6)                                    
1003  format (5x,f10.6,15x,f10.6,15x,f10.6)                                    
      stop                                                                     
      end                                                                      
      subroutine DDX2 (a,b,l,dx,index)                                           
c                                                                               
c  This subroutine performs the second 
c  order finite difference                       
c                                                                               
c  input      : a,l,dx,index                                                         
c  output     : b                                                                    
c  index=0    : cyclic boundary conditions                                       
c  index=1    : noncyclic boundary conditions                                    
c                                                                               
      real a(l),b(l)                                                       
      l1                   = l-1                                                                    
      do 2110 i = 2, l1                                                              
 2110    b(i)   =  (a(i+1) - a(i-1))/(2.*dx)
      if (index.eq.0) b(1) = (a(2)-a(l1))/(2.*dx)                                   
      if (index.eq.0) b(l) = b(1)                                                   
      if (index.eq.1) b(1) = (a(2)-a(1))/dx                                         
      if (index.eq.1) b(l) = (a(l)-a(l1))/dx                                        
      return                                                                    
      end
      subroutine DDX4 (a,b,l,dx)                                                 
c                                                                               
c  This subroutine performs the fourth order 
c  finite difference .The computation is cyclic .                                                    
c                                                                               
c  input  : a,l,dx                                                               
c  output : b                                                                    
      real a(l),b(l)                                                       
      l1        = l-1                                                                    
      l2        = l-2                                                                    
      l3        = l-3                                                                    
      c1        = 4./3.                                                                  
      c2        = 1./3.                                                                  
      dx2       = 1./(2.*dx)                                                            
      dx4       = 1./(4.*dx)
      do 2120 i = 3, l2                                                              
          b(i)  = c1*((a(i+1)-a(i-1))*dx2)
     &          - c2*((a(i+2)-a(i-2))*dx4)                    
 2120 continue                                                                  
      b(2)      = c1*((a(3)-a(1))*dx2)-c2*((a(4)-a(l1))*dx4)                           
      b(1)      = c1*((a(2)-a(l1))*dx2)-c2*((a(3)-a(l2))*dx4)                          
      b(l)      = b(1)                                                                 
      b(l1)     = c1*((a(l)-a(l2))*dx2)-c2*((a(2)-a(l3))*dx4)                         
      return                                                                    
      end