;***********************************************************************
pro zder0
;
;  Calculates the matrix elements and an LU decomposition of the problem
;  of finding the derivatives of a cubic spline with continuous second
;  derivatives.
;
;  This particular version (to be used with the hd/mhd code with fixed
;  boundary values at a fiducial extra boundary layer) has a one-sided
;  first order derivative at the top boundary.
;
;  To make a "natural spline" (zero second derivative at the end pts),
;  change the weights to
;
;     w2(1)=2./3.
;     w3(1)=1./3.
;     w2(nz)=2./3.
;     w1(nz)=1./3.
;
;  For a general cubic spline, we have an equation system of the form:
;
;     | w2 w3          |   | d |     | e |
;     | w1 w2 w3       |   | d |     | e |
;     |       .  .     | * | . |  =  | . |
;     |       .  .  .  |   | . |     | . |
;     |          w1 w2 |   | d |     | e |
;
;  This may be solved by adding a fraction (w1) of each row to the
;  following row.
;
;  This corresponds to an LU decomposition of the problem:
;
;     | w2 w3          |   | d |     | 1               |   | e |
;     |    w2 w3       |   | d |     | w1  1           |   | e |
;     |       .  .     | * | . |  =  |    w1  1        | * | . |
;     |          .  .  |   | . |     |        .  .     |   | . |
;     |             w2 |   | d |     |           w1  1 |   | e |
;
;  In this startup routine, we calculate the fractions w1.
;  These are the same for all splines with the same z-scale,
;  and only have to be calculated once for all nx*ny splines.
;
;  In principle, this startup routine only has to be called once,
;  but to make the zder entry selfcontained, zder0 is called from
;  within zder.  In practice, the zder0 time is insignificant.
; 
;  Update history:
;                      .
;  20-aug-87: coded by Ake Nordlund, Copenhagen University Observatory
;  27-sep-87: Removed bdry[123] parameters.
;  02-oct-87: Moved factor 3 from zder loop to zder0 loop.
;  04-nov-87: inverted w2(), to turn divides into multiplications
;  11-apr-96: first-order one-sided derivative at top
;
;***********************************************************************
;
;
common cdat,x,y,z,nx,ny,nz,nw,ntmax,date0,time0 
common czder,w1,w2,w3
    
	w1=fltarr(nz)
	w2=fltarr(nz)
	w3=fltarr(nz)
;
;  Interior points
;
      a3=1./3.
      w1(1:*)=a3/(z(1:*)-z)
      w3(0:nz-2)=w1(1:nz-1)
      w2=2.*(w1+w3)
;
;  First point
;
      w1(0)=0.
      w2(0)=1.
      w3(0)=0.

;
;  Last point
;
      cza=1./(z(nz-2)-z(nz-3))
      czb=1./(z(nz-1)-z(nz-2))
      w1(nz-1)=cza*cza
      w3(nz-1)=-czb*czb
      w2(nz-1)=w1(nz-1)+w3(nz-1)
;
;  Eliminate at last point
;
      c=-w1(nz-1)/w1(nz-2)
      w2(nz-1)=w2(nz-1)+c*w2(nz-2)
      w3(nz-1)=w3(nz-1)+c*w3(nz-2)
      w1(nz-1)=w2(nz-1)
      w2(nz-1)=w3(nz-1)
      w3(nz-1)=c
;
;  Eliminate subdiagonal elements of rows 2 to nz:
;
;     |     w2 w3       | * | . |  =  | . |
;     |     w1 w2 w3    |   | . |     | . |
;
      for k=1,nz-1 do begin
        w1(k)=-w1(k)/w2(k-1)
        w2(k)=w2(k)+w1(k)*w3(k-1)
	end
;
;  Invert all the w2()'s, turns divides into mults
;
	w2=1./w2
;
      end
;***********************************************************************
function zder,f
;
;  The following LU decomposition of the problem was calculated by the
;  startup routine:
;
;  | w2(1) w3(1)              |   |d|   |  1                    |   |e|
;  |       w2(2) w3(2)        |   |d|   | w1(2)  1              |   |e|
;  |          .     .         | * |.| = |       w1(3)  1        | * |.|
;  |                .      .  |   |.|   |              .  .     |   |.|
;  |                    w2(N) |   |d|   |                 w1  1 |   |e|
;
;  The solutions of all nx*ny equations are carried out in paralell.
;
;  Operation count:  4m+1d+4a = 9 flops per grid point
;
;  Timing:
;     14.2 Mfl = 1.6 Mpts/sec on the Univ of Colo Alliant FX-8.
; 
;  Update history:
;                      .
;  20-aug-87: Coded by Ake Nordlund, Copenhagen University Observatory
;  27-sep-87: Removed bdry[123] parameters.
;  02-oct-87: Moved factor 3 from zder loop to zder0 loop.
;  15-oct-87: Factors a and b moved out of lm-loops 121. 13->9 flops.
;  04-nov-87: inverted w2(), to turn divides into multiplications
;
;***********************************************************************
;
common cdat,x,y,z,nx,ny,nz,nw,ntmax,date0,time0 
common czder,w1,w2,w3
;
;  Check if we have a degenerate case (no z-extension)
;
	d=fltarr(nx,ny,nz)
      if (nz eq 1) then begin
	  d(*,*,*)=0.0
        return,d
      endif
      zder0
;
;  First point
;
      d(*,*,0)=(f(*,*,1)-f(*,*,0))/(z(1)-z(0))
;
;  Interior points [2m+2s = 4 flop]
;
      for k=1,nz-2 do begin
        a=1./(z(k+1)-z(k))^2
        b=1./(z(k)-z(k-1))^2
        d(*,*,k)=(f(*,*,k+1)-f(*,*,k))*a+(f(*,*,k)-f(*,*,k-1))*b
	end
;
;  Last point
;
      cza=1./(z(nz-2)-z(nz-3))
      czb=1./(z(nz-1)-z(nz-2))
      dfa=f(*,*,nz-2)-f(*,*,nz-3)
      dfb=f(*,*,nz-1)-f(*,*,nz-2)
      d(*,*,nz-1)=2.*(dfa*cza*cza*cza-dfb*czb*czb*czb)
;
;  Eliminate at last point
;
      d(*,*,nz-1)=d(*,*,nz-1)+w3(nz-1)*d(*,*,nz-2)
;
;  Do the forward substitution [1m+1a = 2 flop]
;
      for k=1,nz-1 do begin
        d(*,*,k)=d(*,*,k)+w1(k)*d(*,*,k-1)
	end
;
;  Do the backward substitution [1m+1d+1s = 3 flop]
;
      d(*,*,nz-1)=d(*,*,nz-1)*w2(nz-1)
      for k=nz-2,0,-1 do begin
        d(*,*,k)=(d(*,*,k)-w3(k)*d(*,*,k+1))*w2(k)
	end
;
      return,d
	end