;+
; NAME:
;  vt3d_step_impl_trisol_1loc
; PURPOSE: (one line)
;  step the temperature one time step at a single location
; DESCRIPTION:
;  step the temperature one time step at a single location
; CATEGORY:
;  Volatile Transport
; CALLING SEQUENCE:
;  vt3d_step_impl_trisol_1loc, alpha_1, beta_0, 
;     gamma_0, gamma_J, spp_tridc, spp_indx, temp_0, temp_i
; INPUTS:
;  alpha_1: dependance of temp[1] on temp[0], unitless, float
;  beta_0:  dependance of temp[0] on temp[1], unitless, float
;  gamma_0: net flux at upper boundary, K, float
;  gamma_J: net flux at lower boundary, K, float
;  spp_tridc : LU decomposition of the S'' matrix, made by vt3d_cn_tridc
;  spp_tridc : row permutations for partial pivoting to invert S'', made by vt3d_cn_tridc
; INPUT-OUTPUTS:
;  temp_0 : temperature at layer 0, K, float
;  temp_i : temperature at layer 1..J, K, float[J]
; COMMON BLOCKS:
;  None
; SIDE EFFECTS:
; RESTRICTIONS:
;  nsurf = 0 or 1
; PROCEDURE:
;  Solve implicit equation
;  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
;   Resubmitted to Icarus.
;   This is not used directly in Young 2016, but is used as part of
;   the Crank-Nicholson time step, in vt3d_step_cn_trisol_1loc.  
;   See: Eq. 3.4-9b
;   The banded tridiagonal solution is described in Eq. 3.4-12a-6 and 3.4-13a-b
; MODIFICATION HISTORY:
;  Written 2011 Feb 6, by Leslie Young, SwRI
;  2016 Mar 24 LAY. Modified for inclusion in vt3d library.
;-

pro vt3d_step_impl_trisol_1loc, alpha_1, beta_0, gamma_0, gamma_J, spp_tridc, spp_index, temp_0, temp_i

    n_j = n_elements(temp_i)

    temp_this = [temp_0, temp_i]
    temp_this[0] += gamma_0
    temp_this[n_j] += gamma_J

    ; pull apart the pieces of the
    ; LU decomposition
    lower = spp_tridc[0,0:n_j-2]
    diag = spp_tridc[1,0:n_j-1]
    upper = spp_tridc[2,0:n_j-2]
    u2 = spp_tridc[3,0:n_j-3]

    app = replicate(0.d, n_j)
    app[0] = -alpha_1
    y = LA_TRISOL( lower, diag, upper, u2, spp_index, app,/double)
    z = LA_TRISOL( lower, diag, upper, u2, spp_index, temp_this[1:n_j],/double)
    c = - beta_0 * y[0]
    d = - beta_0 * z[0]
       
    hpp = 1 + beta_0
    
    temp_0 = (temp_this[0] - d)/(hpp - c)
    temp_i = (z - y * temp_0)
end