;***********************************************************************
pro yder0
;
;  Periodic y-derivatives
;  1st derivative, 6th order compact (cf. Lele, JCP 103, 16-42)
;  truncation error: f' = (f')_exact + 7.9e-4*dy^6*d^7f/dy^7
;
;  alpha*[f'(i+1)+f'(i-1)] + f'(i) = a*[f(i+1)-f(i-1)] + b*[f(i+2)-f(i-2)]
;  where alpha=1/3, a=14/(9*2*dy), b=1/(9*4*dy)
;
;  rescale whole equation by (9/7)*dy, such that a=1:
;  alpha*[f'(i+1)+f'(i-1)] + beta*f'(i) = [f(i+1)-f(i-1)] + b*[f(i+2)-f[i-2)]
;  alpha=3*dy/7, beta=9*dy/7, and
;  e = (1/28)*[f(i+2)-f(i-2)] + [f(i+1)-f(i-1)]
;
;  Calculates the matrix elements and an LU decomposition of the problem
;  We have an equation system of the form:
;
;     | w2 w3        w1 |   | d |     | e | 
;     | w1 w2 w3        |   | d |     | e |
;     |     .  .  .     | * | . |  =  | . |
;     |        .  .  .  |   | . |     | . |
;     | w3        w1 w2 |   | d |     | e |
;
;  This may be solved by adding a fraction (w4) of each row to the
;  following row, and a fraction (w5) to the last row.
;
;  This corresponds to an LU decomposition of the problem:
;
;     | w2 w3        w1 |   | d |     | 1               |   | e |
;     |    w2 w3     w1 |   | d |     | w4  1           |   | e |
;     |        .  .  .  | * | . |  =  |    w4  1        | * | . |
;     |           .  .  |   | . |     |        .  .     |   | . |
;     |              w2 |   | d |     | w5 w5  .  w4  1 |   | e |
;
;  In this startup routine, we calculate the fractions w4 and w5.
;  These are the same for all periodic derivatives with the same period,
;  and only have to be calculated once for all nx*nz derivatives.
; 
;  In principle, this startup routine only has to be called once,
;  but to make the yder entry selfcontained, yder0 is called from
;  within yder.  In practice, the yder0 time is insignificant.
; 
;  Update history:
;                      .
;  20-aug-87: coded by Ake Nordlund, Copenhagen University Observatory
;  14-sep-87: performance with alternative y-/z-loop ordering tested
;  04-nov-87: inverted w2[], to turn divides into multiplications
;  23-sep-94/db: changed from spline to 6th-order compact derivative
;
;***********************************************************************
;
common cdat,x,y,z,nx,ny,nz,nw,ntmax,date0,time0
common cyder,w1,w2,w3,w4,w5
;
;  set coefficients on lhs: (1, 3, 1)*3*dy/7
;
      dy=y[1]-y[0]
      cy=dy*3./7.
;
      w1=fltarr(ny)
      w2=w1
      w3=w1
      w4=w1
      w5=w1
      w1[*]=cy
      w2[*]=3.*cy
      w3[*]=cy
;
;  Eliminate subdiagonal elements of rows 2 to ny-2
;
      for k=1,ny-3 do begin
        w4[k]=-w1[k]/w2[k-1]
        w2[k]=w2[k]+w4[k]*w3[k-1]
        w1[k]=w4[k]*w1[k-1]
        w5[k]=-w3[ny-1]/w2[k-1]
        w2[ny-1]=w2[ny-1]+w5[k]*w1[k-1]
        w3[ny-1]=w5[k]*w3[k-1]
      end
;
;  We now have a 3*3 system on the last three rows:
;
;     |        w2 w3 w1 | * | . |  =  | . |
;     |        w1 w2 w3 |   | . |     | . |
;     |        w3 w1 w2 |   | d |     | e |
;
      w4[ny-2]=-w1[ny-2]/w2[ny-3]
      w2[ny-2]=w2[ny-2]+w4[ny-2]*w3[ny-3]
      w3[ny-2]=w3[ny-2]+w4[ny-2]*w1[ny-3]
      w5[ny-2]=-w3[ny-1]/w2[ny-3]
      w1[ny-1]=w1[ny-1]+w5[ny-2]*w3[ny-3]
      w2[ny-1]=w2[ny-1]+w5[ny-2]*w1[ny-3]
;
;  Row ny
;
      w4[ny-1]=-w1[ny-1]/w2[ny-2]
      w2[ny-1]=w2[ny-1]+w4[ny-1]*w3[ny-2]
;
;  Invert all the w2[]'s, turns divides into mults
;
      w2[*]=1./w2[*]
;
      end
;***********************************************************************
function yder,f
;
;  The following LU decomposition of the problem was calculated by the
;  startup routine:
;
;  | w2[1] w3[1]        w1[1] |   |d|   |  1                    |   |e|
;  |       w2[2] w3[2]  w1[2] |   |d|   | w4[2]  1              |   |e|
;  |          .     .      .  | * |.| = |       w4[3]  1        | * |.|
;  |                .      .  |   |.|   |              .  .     |   |.|
;  |                    w2[N] |   |d|   | w5[2] w5[3]  .  w4  1 |   |e|
;
;  The solutions of all nx*nz equations are carried out in parallel.
;
;  Update history:
;                      .
;  20-aug-87: coded by Ake Nordlund, Copenhagen University Observatory
;  14-sep-87: performance with alternative y-/z-loop ordering tested
;  04-nov-87: inverted w2[], to turn divides into multiplications
;  25-feb-88: removed use of swapped ff[), speeded things up 10 %
;  04-mar-88: padded blank common with one chunk, for curl, curl2, nodiv
;  07-mar-88: coalesced m,n loops -> mn loop
;  12-apr-88: removed blank common padding for hd (no curl, curl2)
;  23-sep-94/db: changed spline to 6th-order compact derivative
;
;***********************************************************************
;
common cdat,x,y,z,nx,ny,nz,nw,ntmax,date0,time0
common cyder,w1,w2,w3,w4,w5
;
;  Get dimensions, assume 3-dimensional
;
;      s=size(f)
;        nx=s(1) & ny=s(2) & nz=s(3)
;  Check if we have a degenerate case (no y-extension)
;
      d=f
      if ny eq 1 then begin
        d[*,*,*]=0.0
        return,d
      endif
;
;  Do LU decomposition (this is 1-D only)
;
;  Right hand sides, e = [f[i+1]-f[i-1]] + (1/28)*[f[i+2]-f[i-2]]
;
      yder0
      c=1./28.
      ny1=ny-1 & ny2=ny-2 & ny3=ny-3 & ny4=ny-4 & ny5=ny-5
      d[*,0,*]=(f[*,1,*]-f[*,ny1,*])+c*(f[*,2,*]-f[*,ny2,*])
      d[*,1,*]=(f[*,2,*]-f[*,0,*])+c*(f[*,3,*]-f[*,ny1,*])
      d[*,2:ny3,*]=(f[*,3:ny2,*]-f[*,1:ny4,*]) $
                   + c*(f[*,4:ny1,*]-f[*,0:ny5,*])
      d[*,ny2,*]=(f[*,ny1,*]-f[*,ny3,*])+c*(f[*,0,*]-f[*,ny4,*])
      d[*,ny1,*]=(f[*,0,*]-f[*,ny2,*])+c*(f[*,1,*]-f[*,ny3,*])
;
;  Do the forward substitution [2m+2a = 4 flop]
;
      for k=1,ny2 do begin
        d[*,k,*]=d[*,k,*]+w4[k]*d[*,k-1,*]
        d[*,ny1,*]=d[*,ny1,*]+w5[k]*d[*,k-1,*]
      end
      d[*,ny1,*]=d[*,ny1,*]+w4[ny1]*d[*,ny2,*]
;
;  Do the backward substitution [3m+2s = 5 flop]
;
      d[*,ny1,*]=d[*,ny1,*]*w2[ny1]
      d[*,ny2,*]=(d[*,ny2,*]-w3[ny2]*d[*,ny1,*])*w2[ny2]
      for k=ny-3,0,-1 do begin
        d[*,k,*]=(d[*,k,*]-w3[k]*d[*,k+1,*]-w1[k]*d[*,ny1,*])*w2[k]
      end
;
      return,d
      end