;+
; NAME:
;  vt3d_zdelta
; PURPOSE: (one line)
;  Return depths and layer widths for volatile transport
; DESCRIPTION:
;  Return depths and layer widths for volatile transport
; CATEGORY:
;  Volatile Transport
; CALLING SEQUENCE:
;  vt3d_zdelta, z, z_delta, z_delta_top, z_delta_bot, k=k, z_delta_z0 = z_delta_z0
; INPUTS:
;  z - cm. depth of grid points (z <= 0)
;        1-D, length n_z, if skindepth is a scalar, 
;        2-D, [n_loc, n_z] if skindepth is an array
; OPTIONAL INPUT PARAMETERS:
;  k - thermal conductivity
;  z_delta0 - width of the top layer
; KEYWORD INPUT PARAMETERS:
; KEYWORD OUTPUT PARAMETERS:
;  none
; OUTPUTS:
; z_delta[i] (cm) is the spacing between the top and bottom of layer i.
; z_delta_top[i] (cm) is s.t.
;   flux = - k[i] (T[i-1] - T[i]) / z_delta_top[i]
;   is continuous at the boundary.
; z_delta_bot[i] (cm) is s.t.
;   flux = - k[i] (T[i] - T[i+1]) / z_delta_bot[i]
;   is continuous at the boundary.
; COMMON BLOCKS:
;  None
; SIDE EFFECTS:
; RESTRICTIONS:
;  nsurf = 0 or 1
; PROCEDURE:
;  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
;   Resubmitted to Icarus.
;   Section 3.3
; MODIFICATION HISTORY:
;  Written 2011 Feb 6, by Leslie Young, SwRI
;  2014 May 18 Modify for z being 1D.  LAY.
;-

pro vt3d_zdelta, z, z_delta, z_delta_top, z_delta_bot, k=k, z_delta_z0 = z_delta_z0

    if size(z,/n_dim) eq 1 then begin
       z2d = reform(z,1,n_elements(z))
       if keyword_set(k) then begin
          if size(k,/n_dim) eq 1 then begin
             k2d = reform(k,1,n_elements(z))
          endif else begin
             k2d = replicate(k,1,n_elements(z))
          endelse
          zdelta, z2d, z_delta, z_delta_top, z_delta_bot, k=k2d, z_delta_z0 = z_delta_z0
       endif else begin
          zdelta, z2d, z_delta, z_delta_top, z_delta_bot, z_delta_z0 = z_delta_z0
       endelse
       z_delta = reform(z_delta)
       z_delta_top = reform(z_delta_top)
       z_delta_bot = reform(z_delta_bot)
       return
    endif

    nxnyn, z, n_loc, n_z
    
    z_delta = fltarr(n_loc,n_z)
    
    if keyword_set(z_delta_z0) then begin
       z_delta[*,0] = z_delta_z0
    endif else begin
       z_delta[*,0] = (z[*,0] - z[*,1])/2.
    endelse

    ; top layer is special
    z_delta[*,1] = 2 * (z[*,0] - z_delta[*,0] - z[*,1])

;    z_delta[*,0] = 0.
    for i = 2, n_z - 1 do begin
       ; z(i-1) - dz(i-1)/2 = z(i) + dz(i)/2
       z_delta[*,i] = 2 * (z[*,i-1] - z[*,i] - z_delta[*,i-1]/2)
    endfor

    z_delta_top = fltarr(n_loc,n_z)
    z_delta_bot = fltarr(n_loc,n_z)
    z_delta_top = z_delta/2
    z_delta_bot = z_delta/2
    z_delta_top[0] = 0
    z_delta_bot[0] = z_delta[0]
    if keyword_set(k) then begin
       z_delta_top[*,1] += (k[*,1]/k[*,0]) * z_delta[*,0]       
       z_delta_top[*,2:n_z-1] += 0.5*(k[*,2:n_z-1]/k[*,0:n_z-2]) * z_delta[*,1:n_z-2]       
       z_delta_bot[*,0:n_z-2] += 0.5*(k[*,0:n_z-2]/k[*,1:n_z-1]) * z_delta[*,1:n_z-1]       
    endif else begin
       z_delta_top[*,1] += z_delta[*,0] 
       z_delta_top[*,2:n_z-1] += 0.5 * z_delta[*,1:n_z-2] 
       z_delta_bot[*,0:n_z-2] += 0.5 * z_delta[*,1:n_z-1]
    endelse

end