;+
; NAME:
;  vt3d_step_impl_nloc
; PURPOSE: (one line)
;  step the temperature one time step
; DESCRIPTION:
;  step the temperature one time step for multiple locations with
;  local mass balance
; CATEGORY:
;  Volatile Transport
; CALLING SEQUENCE:
;  vt3d_step_impl_nloc, alpha_i, beta_0, beta_i, gamma_0, gamma_J, temp_0, temp_i
; INPUTS:
;  alpha_i: dependance of temp[1..J] on temp[0..J-1], unitless, float[J]
;  beta_0:  dependance of temp[0] on temp[1], unitless, float[n_loc]
;  beta_i:  dependance of temp[1..J] on temp[2..J+1], unitless, float[J]
;  gamma_0: net flux at upper boundary, K, float, float[n_loc]
;  gamma_J: net flux at lower boundary, K, float
; INPUT-OUTPUTS:
;  temp_0 : temperature at layer 0, K, float[n_loc]
;  temp_i : temperature at layer 1..J, K, float[n_loc,J]
; COMMON BLOCKS:
;  None
; SIDE EFFECTS:
; RESTRICTIONS:
;  nsurf = 0 or 1
; PROCEDURE:
;  
; MODIFICATION HISTORY:
;  Written 2011 Feb 6, by Leslie Young, SwRI
;-

pro vt3d_step_impl_nloc, alpha_i, beta_0, beta_i, gamma_0, gamma_J, temp_0, temp_i

    method = 'mat'

    if method eq 'loop' then begin
       n_loc = n_elements(beta_0)
       for i_loc = 0, n_loc-1 do begin

          temp_0l = temp_0[i_loc]
          temp_il = reform(temp_i[i_loc,*])

          vt3d_step_impl_1loc, alpha_i, beta_0[i_loc], beta_i, gamma_0[i_loc], gamma_J, temp_0l, temp_il

          temp_0[i_loc] = temp_0l
          temp_i[i_loc,*] = temp_il

       endfor
    endif

    if method eq 'mat' then begin
       n_loc = n_elements(beta_0)
       n_j = n_elements(alpha_i)

       ; make the LU decomposition of Spp
       lower = -alpha_i[1:n_j-1]
       diag = 1 + (alpha_i + beta_i)
       upper = -beta_i[0:n_j-2]
       LA_TRIDC, lower, diag, upper, u2, spp_index       

       ; get y_{L} = Spp x app
       app = replicate(0.,n_j)
       app[0] = - alpha_i[0]
       y = LA_TRISOL(lower, diag, upper, u2, spp_index, app)
       
       ; get Z
       temp_tilda_0 = temp_0 + gamma_0
       temp_tilda_i = temp_i
       temp_tilda_i[*, n_j-1] += gamma_J
       z = LA_TRISOL(lower, diag, upper, u2, spp_index, temp_tilda_i)
       hpp = 1 + beta_0
       c = - beta_0 * y[0]
       d = - beta_0 * z[*,0]
       temp_0 = (temp_tilda_0 - d) / (hpp - c)
       temp_i = z - y ## temp_0

    endif


end