Extended IDL Help

This page was created by the IDL library routine mk_html_help. For more information on this routine, refer to the IDL Online Help Navigator or type:

     ? mk_html_help

at the IDL command line prompt.

Last modified: Thu Jul 7 23:39:43 2016.


List of Routines


Routine Descriptions

AAVT_README

[Next Routine] [List of Routines]
 NAME:
   AAvt_README
 PURPOSE
   Print the file 'AAvt_README.pro'
 CATEGORY:
  Volatile transport and thermophysical modeling (vt)
 CALLING SEQUENCE:
  AAvt_README
 
 ---------------------------------------------------------------------------
 FILES AND ENVIRONMENT VARIABLES
 ---------------------------------------------------------------------------
 
 None

 ---------------------------------------------------------------------------
 FUNCTIONS
 ---------------------------------------------------------------------------

 A help file was  

 mk_html_help,['.','vty16'],'AAvt.html'

 ---------------------------------------------------------------------------
 VARIABLES
 ---------------------------------------------------------------------------
 
 VT16
 ----
 In vty16, there are some common variable names used in multiple
 functions.
 ARRAY LENGTHS
 n_z : number of layers

 ---------------------------------------------------------------------------
 FUNCTIONS FROM ELSEWHERE IN THE LAYOUNG LIBRARY
 ---------------------------------------------------------------------------
 general/physconstants
 ---------------------------------------------------------------------------

(See AAvt_README.pro)


BZ1980_VP

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  bz1980_vp
 PURPOSE: (one line)
  Return vapor pressure, optionally latent heat
 DESCRIPTION:
  Return vapor pressure, optionally latent heat from Brown & Ziegler 1980
 CATEGORY:
  Volatile Transport
 CALLING SEQUENCE:
  vp = bz1980_vp(species, t, latentheat, energyratio, phase)
 INPUTS:
  species: string 
    Valid species are: p-H2, O2, N2, F2, CO, CO2, Ne, Ar, Kr, Xe, CH4, C2H4, C2H6
  t : temperature, K
 OUTPUTS:
  vp : vapor pressure in microbar
 OPTIONAL OUTPUT
  latentheat : latent heat of sublimation or vaporization in erg/g
  energyratio : ratio of latent to kinetic energy, unitless.
    latent heat per molecule / (k*T)
  phase : a string representing the phase
  key: string giving designation in  bz1980_vp_table3.txt
  ap: the polynomial coefficients A0..A6 in bz1980_vp_table3.txt
 RESTRICTIONS:
  The files "bz1980_vp_table1.txt" and "bz1980_vp_table3.txt"
   must be in the same directory as bz1980_vp.pro
 PROCEDURE:
  This function implements the tables of 
   Brown, G.N., Ziegler, W.T., 1980.
   Vapor pressure and heats of vaporization and
   sublimation of liquids and solids of interest in cryogenics below 1atm
    pressure. Adv. Cryog. Eng. 25, 662–670.
   (note there is confustion about the year of this publication, as
   the conference was 1979 and the procedings were published in
   1980).
  Latent heat is calculated from the derivative of the vapor
  pressure, as John Stansberry pointed out that the BZ latent heat
  tables are incorrect.
 MODIFICATION HISTORY:
  Written 2011 Apr 10, by Leslie Young, SwRI
  2016 Mar 22 Cleaned up documentation, LAY 
  2016 June 25 enforce output to be same dimensions as t

(See bz1980_vp.pro)


CH4LATENTHEAT

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  ch4latentheat
 PURPOSE: (one line)
  Return latent heat of sublimation of CH4 in erg/g from Brown & Ziegler 1980
 CATEGORY:
  Volatile Transport
 CALLING SEQUENCE:
  latentheat = ch4latentheat(t)
 INPUTS:
  t : temperature, K
 KEYWORD INPUTS
  None
 OUTPUTS:
  latentheat : latent heat of sublimation in erg/g
 OPTIONAL OUTPUT
  None
 RESTRICTIONS:
 PROCEDURE:
  This function implements the tables of Brown and Zeigler 
 MODIFICATION HISTORY:
  Written 2011 Apr 10, by Leslie Young, SwRI

(See ch4latentheat.pro)


CH4T

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  ch4t
 PURPOSE: (one line)
  Return vapor pressure of CH4 from Brown & Ziegler 1980
 CATEGORY:
  Volatile Transport
 CALLING SEQUENCE:
  t = ch4t(vp)
 INPUTS:
  vp : vapor pressure in microbar
 KEYWORD INPUTS
  None
 OUTPUTS:
  t : temperature, K
 OPTIONAL OUTPUT
  None
 RESTRICTIONS:
 PROCEDURE:
  This function implements the tables of Brown and Zeigler 
 MODIFICATION HISTORY:
  Written 2013 Jun 16, by Leslie Young, SwRI

(See ch4t.pro)


CH4VP

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  ch4vp
 PURPOSE: (one line)
  Return vapor pressure of CH4 from Brown & Ziegler 1980
 CATEGORY:
  Volatile Transport
 CALLING SEQUENCE:
  vp = ch4vp(t)
 INPUTS:
  t : temperature, K
 KEYWORD INPUTS
  None
 OUTPUTS:
  vp : vapor pressure in microbar
 OPTIONAL OUTPUT
  None
 RESTRICTIONS:
 PROCEDURE:
  This function implements the tables of Brown and Zeigler 
 MODIFICATION HISTORY:
  Written 2013 Jun 13, by Leslie Young, SwRI

(See ch4vp.pro)


COLATENTHEAT

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  colatentheat
 PURPOSE: (one line)
  Return latent heat of sublimation of CO in erg/g from Brown & Ziegler 1980
 CATEGORY:
  Volatile Transport
 CALLING SEQUENCE:
  latentheat = colatentheat(t)
 INPUTS:
  t : temperature, K
 KEYWORD INPUTS
  None
 OUTPUTS:
  latentheat : latent heat of sublimation in erg/g
 OPTIONAL OUTPUT
  None
 RESTRICTIONS:
 PROCEDURE:
  This function implements the tables of Brown and Zeigler 
 MODIFICATION HISTORY:
  Written 2011 Apr 10, by Leslie Young, SwRI

(See colatentheat.pro)


COVP

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  covp
 PURPOSE: (one line)
  Return vapor pressure of CO from Brown & Ziegler 1980
 CATEGORY:
  Volatile Transport
 CALLING SEQUENCE:
  vp = covp(t)
 INPUTS:
  t : temperature, K
 KEYWORD INPUTS
  None
 OUTPUTS:
  vp : vapor pressure in microbar
 OPTIONAL OUTPUT
  None
 RESTRICTIONS:
 PROCEDURE:
  This function implements the tables of Brown and Zeigler 
 MODIFICATION HISTORY:
  Written 2009 Dec 7, by Leslie Young, SwRI

(See covp.pro)


N2LATENTHEAT

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  n2latentheat
 PURPOSE: (one line)
  Return latent heat of sublimation of N2 in erg/g from Brown & Ziegler 1980
 CATEGORY:
  Volatile Transport
 CALLING SEQUENCE:
  latentheat = n2latentheat(t)
 INPUTS:
  t : temperature, K
 KEYWORD INPUTS
  None
 OUTPUTS:
  latentheat : latent heat of sublimation in erg/g
 OPTIONAL OUTPUT
  None
 RESTRICTIONS:
 PROCEDURE:
  This function implements the tables of Brown and Zeigler 
 MODIFICATION HISTORY:
  Written 2011 Jan 25, by Leslie Young, SwRI

(See n2latentheat.pro)


N2T

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  n2t
 PURPOSE: (one line)
  Return vapor pressure of N2 from Brown & Ziegler 1980
 CATEGORY:
  Volatile Transport
 CALLING SEQUENCE:
  t = n2t(vp)
 INPUTS:
  vp : vapor pressure in microbar
 KEYWORD INPUTS
  None
 OUTPUTS:
  t : temperature, K
 OPTIONAL OUTPUT
  None
 RESTRICTIONS:
 PROCEDURE:
  This function implements the tables of Brown and Zeigler 
 MODIFICATION HISTORY:
  Written 2013 Jun 16, by Leslie Young, SwRI

(See n2t.pro)


N2VP

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  n2vp
 PURPOSE: (one line)
  Return vapor pressure of N2 from Brown & Ziegler 1980
 CATEGORY:
  Volatile Transport
 CALLING SEQUENCE:
  vp = n2vp(t)
 INPUTS:
  t : temperature, K
 KEYWORD INPUTS
  None
 OUTPUTS:
  vp : vapor pressure in microbar
 OPTIONAL OUTPUT
  None
 RESTRICTIONS:
 PROCEDURE:
  This function implements the tables of Brown and Zeigler 
 MODIFICATION HISTORY:
  Written 2011 Aug 24, by Leslie Young, SwRI

(See n2vp.pro)


VT3D_1LOC_DIURNAL_LOCAL

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vt3d_1loc_diurnal_local
 PURPOSE: (one line)
  Calculate the diurnal variation at a single non-interacting location
 DESCRIPTION:
  Calculate the diurnal variation at a single non-interacting location
 CATEGORY:
  Volatile Transport 
 CALLING SEQUENCE:
  vt3d_1loc_diurnal_bare, constants, input, grid, program, output
 INPUTS:
 constants: constants structure, with
    i                   ; sqrt(-1)
    deg                 ; radian per deg
    hour                ; s per hour
    year                ; s per year
    km                  ; cm per km
    ti_SI               ; [cgs thermal inertia] per [SI thermal inertia]
    m_u                 ; atomic mass constant, g
    N_a                 ; Avogadro constant, mol^-1
    k                   ; Boltzmann constant, erg K^-1
    sigma               ; Stefan-Boltzmann constant erg cm^-2 K^-4 s^-1
    G                   ; gravitational constant, cm^3 g^-1 s^-2
    sol_norm_1au        ; solar flux at 1 AU, erg/cm^2/s
 input : inputs defining the problem
    dist_sol_au         ; heliocentric distance in AU
    albedo              ; bond albedo, assumed equal to hemispheric
    emis                ; emissivity, scalar
    flux_int            ; internal heat flux, erg cm^-2 s^-1
    therminertia        ; thermal inertia in erg cm^-2 s^1/2 K^-1, scalar
    dens                ; density, g cm^-3, scalar
    specheat            ; specific heat in erg g^-1 K^-1, scalar
    thermcond           ; thermal conductivity in erg cm^-1 s^-1 K^-1, scalar
    lon_sol_ref         ; reference sub-solar longitude, radian
    t_lon_sol_ref       ; time of lon_sol_ref, s
    dlon_sol_dt         ; rate of change of lon_sol, radian/s
                        ; (negative for rotation N pole, E longitude)
    lat_sol             ; sub-solar latitude, radian
    is_volatile         ; 1 if has a volatile, 0 if bare. flag. scalar
    mass_volatile       ; mass of the volatile slab, g cm^-2, scalar
    specheat_volatile   ; specific heat of the volatile slab, erg g^-1 K^-1, scalar
    gravacc             ; gravitational acceleration, cm s^-2
    name_species        ; name of the species.  String. "N2", "CO", or "CH4"
    mflux_esc           ; escape rate in g cm^-2 s^-1, constant with time
 grid : grids 
    lat                 ; latitude of spot, radian, scalar
    lon                 ; longitude of spot, radian, scalar
    t_start             ; time of calculation start, s, scalar
    t_end               ; time of calculation end, s, scalar
    tau                 ; unitless timestep, (delta time) * (frequency)
    zeta                ; unitless depths of layers, cm, [nz]
 program : program control
    n_terms             ; number of terms in the sinusoidal expansion
    iterate_temp_terms  ; 1 to iterate the explicit temp terms
    iterate_dfluxdtemp_e ; iterate the flux-per-temperature for emission
    explicit            ; 1 if explicit timestep, 0 if implicit (Crank-Nicholson)
    debug               ; 1 to print debug statements, 0 otherwise
 OUTPUT:
 output: output
    sol_terms : terms in solar expansion, erg cm^-2 s^-1, complex[nterm+1]
    temp_terms : terms in temperature expansion, K, complex[nterm+1]
    s : array of states
      t : time
      temp_0 : temperature of layer 0, K, scalar
      temp_1J : temperature of layers 1..J, K, [nz]
 SIDE EFFECTS:
 RESTRICTIONS:
 PROCEDURE:
  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
   Resubmitted to Icarus.
   Section 3 (is_volatile=0) & 4 (is_volatile=1)

(See vt3d/vt3d_1loc_diurnal_local.pro)


VT3D_ALPHA_INTERIOR

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vt3d_alpha_interior
 PURPOSE: (one line)
  Return alpha[1..J] for volatile transport timestepping
 DESCRIPTION:
  Return alpha[1..J] for volatile transport timestepping
 CATEGORY:
  Volatile Transport
 CALLING SEQUENCE:
  alpha_i = vt3d_alpha_interior(tau, del, del_a)
 INPUTS:
  tau : non-dimentionalized time step, scalar
  del : non-dimentionalized layer width, float[J+1]
 OPTIONAL INPUTS
  del_a : non-dimentionalized distance to layer above, float[J+1]
 OUTPUTS:
  alpha_i : dependance of temp[1..J] on temp[0..J-1], float[J]
 RESTRICTIONS:
 PROCEDURE:
  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
   Resubmitted to Icarus.
   Table 4 & 5
 MODIFICATION HISTORY:
  Written 2011 Dec 30, by Leslie Young, SwRI
  2016 Mar 24 LAY. Modified for inclusion in vt3d library.

(See vt3d/vt3d_alpha_interior.pro)


VT3D_BETA_INTERIOR

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vt3d_beta_interior
 PURPOSE: (one line)
  Return beta[1..J] for volatile transport timestepping
 DESCRIPTION:
  Return beta[1..J] for volatile transport timestepping
 CATEGORY:
  Volatile Transport
 CALLING SEQUENCE:
  beta_i = vt3d_beta_interior(tau, del, del_b)
 INPUTS:
  tau : non-dimentionalized time step, scalar
  del : non-dimentionalized layer width, float[J+1]
 OPTIONAL INPUTS
  del_b : non-dimentionalized distance to layer below, float[J+1]
 OUTPUTS:
  beta_i : dependance of temp[1..J] on temp[2..J+1], float[J]
  beta_[0:n_z-3] are defined; final element is zero.
 RESTRICTIONS:
 PROCEDURE:
  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
   Resubmitted to Icarus.
   Table 4 & 5
 MODIFICATION HISTORY:
  Written 2011 Dec 30, by Leslie Young, SwRI
  2016 Mar 24 LAY. Modified for inclusion in vt3d library.

(See vt3d/vt3d_beta_interior.pro)


VT3D_CN_TRIDC

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vt3d_cn_tridc
 PURPOSE: (one line)
  Upper-Lower decomposition of the substrate matrix, S''
 DESCRIPTION:
  Upper-Lower decomposition of the substrate matrix, S''
 CATEGORY:
  Volatile Transport
 CALLING SEQUENCE:
  spp_tridc = vt3d_cn_tridc(alpha_i, beta_i, spp_indx)
 INPUTS:
  alpha_i: dependance of temp[1..J] on temp[0..J-1], unitless, float[n_j]
  beta_i:  dependance of temp[1..J] on temp[2..J+1], unitless, float[n_j]
 OUTPUTS:
   spp_tridc - the tridiagonal LU decomposition of the substrate matrix, S.
    spp_tridc[0,0:n_j-2] = spp_al 
      lower bidiagonal array.  float[n_j]
    spp_tridc[1,0:n_j-1] = spp_a  
      diagonal elements of the upper array. float[n_j+1]
    spp_tridc[2,0:n_j-2] = spp_au 
      superdiagonal elements of the upper array. float[n_j]
    spp_tridc[3,0:n_j-3] = u2 
      Second superdiagonal of the upper array. float(n_j-1)
   spp_index An output vector that records the row permutations which occurred as a result of partial pivoting. 
 PROCEDURE:
  Solve Crank Nickleson implicit  equation
  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
   Resubmitted to Icarus.
   This is used as part of the Crank-Nicholson time step.  See:
   Table 4 & 5, Figure 3-6, Figure 3-8
  The LU decomposition is done using LA_TRIDC
 MODIFICATION HISTORY:
  Written 2011 Feb 6, by Leslie Young, SwRI
  2016 Mar 25 LAY. Modified for inclusion in vt3d library.

(See vt3d/vt3d_cn_tridc.pro)


VT3D_DFLUXDTEMP_ATM

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vt3d_dfluxdtemp_atm
 PURPOSE: (one line)
  Return (d flux / d T) associated with atmosphere
 DESCRIPTION:
  Return (d flux / d T) associated with atmosphere, latent heat
 CATEGORY:
  Volatile Transport
 CALLING SEQUENCE:
  phi_a = vt3d_dfluxdtemp_atm(freq, temp_v, frac_varea, gravacc, name_species)
 INPUTS:
  freq : frequency in s^-1
  temp_v : volatile frost temperature, 
  frav_varea : fraction of area covered by volatiles, unitless
   Always 1 for local mass conservation (isolated volatile-covered areas)
  gravacc : effective grav. accelleration, cm/s^2, defined as mass_atm / p_surface
   gravacc = (GM/R^2) * (1 - 2 H/R)
   where G is the gravitational constant
   M is the mass of the body
   R is the surface radius
   H is the scale height
  name_species : string. "N2", "CO", or "CH4"
 OUTPUTS:
  phi_a : d flux/d T in erg cm^-2 s^-1 K^-1
 OPTIONAL OUTPUTS:
  latheat : latent heat (erg/g)
  a : latent heat / kinetic energy ratio, ( (erg/molecule) / (k T) )
 RESTRICTIONS:
  name_species must be one of 'N2', 'CO', 'CH4'
 PROCEDURE:
  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
   Resubmitted to Icarus.
   Eq. 4.1-5 and 4.2-4b for isolated areas (frac_varea = 1)
   Eq. 4.1-5 and 5.2-2d for interactive areas (frac_varea = 0 to 1)
 MODIFICATION HISTORY:
  Written 2011 Apr 10, by Leslie Young, SwRI
  2016 Mar 24 LAY. Modified for inclusion in vt3d library.

(See vt3d/vt3d_dfluxdtemp_atm.pro)


VT3D_DFLUXDTEMP_EMIT

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vt3d_dfluxdtemp_emit
 PURPOSE: (one line)
  Return (d flux / d T) associated with thermal emission
 DESCRIPTION:
  Return (d flux / d T) associated with thermal emission
 CATEGORY:
  Volatile Transport
 CALLING SEQUENCE:
  phi_e = vt3d_dfluxdtemp_emit(emis, temp)
 INPUTS:
  emis : emissivity, unitless
  temp : temperature, K
 OUTPUTS:
  phi_e : d flux/d T in erg cm^-2 s^-1 K^-1
    where flux is the thermal emission
 RESTRICTIONS:
 PROCEDURE:
  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
   Resubmitted to Icarus.
   Eq. 3.2-9b
 MODIFICATION HISTORY:
  Written 2011 Apr 10, by Leslie Young, SwRI
  2011-05-18 LAY.  change from e*s*T^3 to 4*e*s*T^3
  2011-12-30 LAY.  change from temp_0 to temp
  2016 Mar 24 LAY. Modified for inclusion in vt3d library.

(See vt3d/vt3d_dfluxdtemp_emit.pro)


VT3D_DFLUXDTEMP_SLAB

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vt3d_dfluxdtemp_slab
 PURPOSE: (one line)
  Return (d flux / d T) associated with slab
 DESCRIPTION:
  Return (d flux / d T) associated with slab specific heat
 CATEGORY:
  Volatile Transport
 CALLING SEQUENCE:
  phi_v = vt3d_dfluxdtemp_slab(freq, mass_0, specheat)
 INPUTS:
  freq : frequency in s^-1
  mass_0 : time-averaged slab mass, averaged over the volatiles. g cm^-2
  specheat : specific heat, averaged over the volatiles, erg g^-1 K^-1
    specheat is defined as
      total( specheat * mass_0 * angarea_delta ) / total( mass_0 * angarea_delta )
 OUTPUTS:
  phi_v : d flux/d T in erg cm^-2 s^-1 K^-1
 RESTRICTIONS:
 PROCEDURE:
  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
   Resubmitted to Icarus.
   Eq. 4.2-4a
 MODIFICATION HISTORY:
  Written 2011 Apr 10, by Leslie Young, SwRI
  2011-11-05 LAY slightly changed documentation
  2016 Jul 3 LAY. Modified for inclusion in vt3d library.

(See vt3d/vt3d_dfluxdtemp_slab.pro)


VT3D_DFLUXDTEMP_SUBSTRATE

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vt3d_dfluxdtemp_substrate
 PURPOSE: (one line)
  Return (d flux / d T) associated with thermal inertia
 DESCRIPTION:
  Return (d flux / d T) associated with thermal inertia
 CATEGORY:
  Volatile Transport
 CALLING SEQUENCE:
  phi_s = vt3d_dfluxdtemp_substrate(freq, therminertia)
 INPUTS:
  freq : frequency in s^-1
  therminertia : thermal inertia in erg cm^-2 s^-0.5 K^-1
 OUTPUTS:
  phi_gamma : d flux/d T in erg cm^-2 s^-1 K^-1
    where flux is the thermal conduction
 RESTRICTIONS:
 PROCEDURE:
  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
   Resubmitted to Icarus.
   Eq. 3.2-9a
 MODIFICATION HISTORY:
  Written 2011 Apr 10, by Leslie Young, SwRI
  2016 Mar 24 LAY. Modified for inclusion in vt3d library.

(See vt3d/vt3d_dfluxdtemp_substrate.pro)


VT3D_EFLUX_TERMS_BARE

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vt3d_eflux_terms_bare
 PURPOSE: (one line)
  Return temperature in balance with mean absorbed solar flux
 DESCRIPTION:
  Return temperature in balance with mean absorbed solar flux
 CATEGORY:
  Volatile Transport
 CALLING SEQUENCE:
  flux_terms = vt3d_eflux_terms_bare(sol_terms,flux_int,emis,freq,therminertia)
 INPUTS:
  sol_terms : absorbed solar flux, complex or real, erg/cm^2/s:
   complex[n_term+1] or complex[n_loc,n_term+1]
     sol = sol_terms[0] + sol_terms[1]*exp(_i * freq * t) 
                    + sol_terms[2]*exp(2 *_i freq t) + ..
     
  flux_int : internal hear flux, erg/cm^2/s: float or float[n_loc]
  emis : thermal emissivity, unitless : float or float[n_loc]
  freq : frequency of forcing
  therminertia : thermal inertia, erg cm-2 s-1/2 K-1 : float or float[n_loc]
 OUTPUTS:
  flux_terms: terms of approximate emitted flux, erg cm^-2 s^-1, s.t.
     flux = flux_terms[0] + flux_terms[1]*exp(   _i * freq * t + sqrt(    _i * zeta) ) 
                      + flux_terms[2]*exp(2 *_i * freq * t + sqrt(2 * _i * zeta) ) + ..
  flux_terms is complex[n_term+1] or complex[n_loc,n_term+1]
 OPTIONAL OUTPUTS:
  temp_terms: terms of approximate temperture, K, s.t.
     temp = temp_terms[0] + temp_terms[1]*exp(   _i * freq * t + sqrt(    _i * zeta) ) 
                      + temp_terms[2]*exp(2 *_i * freq * t + sqrt(2 * _i * zeta) ) + ..
  temp_terms is complex[n_term+1] or complex[n_loc,n_term+1]
 RESTRICTIONS:
 PROCEDURE:
  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
   Resubmitted to Icarus.
   Eq 3.2-17
 MODIFICATION HISTORY:
  Written 2011 Dec 31, by Leslie Young, SwRI

(See vt3d/vt3d_eflux_terms_bare.pro)


VT3D_EFLUX_TERMS_LOCAL

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vt3d_eflux_terms_local
 PURPOSE: (one line)
  Return temperature in balance with mean absorbed solar flux
 DESCRIPTION:
  Return temperature in balance with mean absorbed solar flux
 CATEGORY:
  Volatile Transport
 CALLING SEQUENCE:
   flux_terms = vt3d_eflux_terms_local(sol_terms,flux_int,emis,freq,therminertia, $
    is_volatile, mass_volatile, specheat_volatile, gravacc, name_species)
 INPUTS:
  sol_terms : absorbed solar flux, complex or real, erg/cm^2/s:
   complex[n_term+1] or complex[n_loc,n_term+1]
     sol = sol_terms[0] + sol_terms[1]*exp(_i * freq * t) 
                    + sol_terms[2]*exp(2 *_i freq t) + ..
     
  flux_int : internal hear flux, erg/cm^2/s: float or float[n_loc]
  emis : thermal emissivity, unitless : float or float[n_loc]
  freq : frequency of forcing
  therminertia : thermal inertia, erg cm-2 s-1/2 K-1 : float or float[n_loc]
  is_volatile : 1 if has a volatile, 0 if bare. flag.  int or int[n_loc]
  mass_volatile : mass of the volatile slab, g cm^-2 : float or float[n_loc]
  specheat_volatile : specific heat of the volatile slab, erg g^-1 K^-1, scalar
  gravacc : gravitational acceleration, cm s^-2
  name_species : name of the species.  String. "N2", "CO", or "CH4"
 OPTIONAL INPUTS:
  mflux_esc_terms     ; escape rate in g cm^-2 s^-1
   complex[n_term+1] or complex[n_loc,n_term+1]
     defined so that (if _i = sqrt(-1) )
     mflux_esc = mflux_esc_terms[0] + mflux_esc_terms[1]*exp(   _i * freq * t) 
                                    + mflux_esc_terms[2]*exp(2 *_i * freq * t) + ..
 OUTPUTS:
  flux_terms: terms of approximate emitted flux, erg cm^-2 s^-1, s.t.
     flux = flux_terms[0] + flux_terms[1]*exp(   _i * freq * t + sqrt(    _i * zeta) ) 
                      + flux_terms[2]*exp(2 *_i * freq * t + sqrt(2 * _i * zeta) ) + ..
  flux_terms is complex[n_term+1] or complex[n_loc,n_term+1]
 OPTIONAL OUTPUTS:
  temp_terms: terms of approximate temperture, K, s.t.
     temp = temp_terms[0] + temp_terms[1]*exp(   _i * freq * t + sqrt(    _i * zeta) ) 
                      + temp_terms[2]*exp(2 *_i * freq * t + sqrt(2 * _i * zeta) ) + ..
  temp_terms is complex[n_term+1] or complex[n_loc,n_term+1]
 RESTRICTIONS:
  If defined, mflux_esc_terms must be the same dimensions as sol_terms
 PROCEDURE:
  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
   Resubmitted to Icarus.
   Eq 4.2-10
 MODIFICATION HISTORY:
  Written 2011 Dec 31, by Leslie Young, SwRI

(See vt3d/vt3d_eflux_terms_local.pro)


VT3D_LOCAVG

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vt3d_locavg
 PURPOSE: (one line)
  Return spatial average over location
 DESCRIPTION:
  Return spatial average over location
 CATEGORY:
  Util
 CALLING SEQUENCE:
  avg = vt3d_locavg(val, weight, indx=indx, default=default)
 INPUTS:
  val = a scalar, array[n_loc], or array[n_loc,n_time]
  weight = array float[n_loc] 
  - usually angular area of surface elements (ster)
 OPTIONAL INPUT PARAMETER
  indx = array long[n_indx], 0 <= indx < nloc.  
    If not passed, all indices are used.
    If indx = -1, use "default"
  default = returned value if indx = -1 (e.g., do not choose any elements)
 OUTPUTS:
  average of val or val[indx], weighted by area (weight)
 PROCEDURE
  At its core this is just total(val * weight)/ total(weight)
  This is very similar to wmean, except:
  1. it accepts weights instead of errors
  2. it can accpet the indx keyword
  3. the names are taylored for weighting a map
  4. If val is 2-D, then average over the first dimension only
  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
   Resubmitted to Icarus.
   Eq. 5.3-3
 MODIFICATION HISTORY:
  Written 2011 Aug 3, by Leslie Young, SwRI
  2016 Apr 21, moved into layoung IDL library

(See vt3d/vt3d_locavg.pro)


VT3D_SKINDEPTH

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
 vt3d_skindepth
 PURPOSE: (one line)
  Return vt3d_skindepth
 DESCRIPTION:
  Return vt3d_skindepth
 CATEGORY:
  Volatile Transport
 CALLING SEQUENCE:
  z_skin = vt3d_skindepth(dens, specheat, thermcond, freq)
 INPUTS:
  dens : density of the substrate (g/cm^3)
  specheat : specific heat (erg/(g K))
  thermcond : thermal conductivity of the substrate (erg/(cm s K))
  freq : frequency of the solar forcing, 2 pi/period (1/(s))
 OUTPUTS:
  z_skin : skin depth (cm)
 RESTRICTIONS:
 PROCEDURE:
  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
   Resubmitted to Icarus.
   Eq 3.2-6
 MODIFICATION HISTORY:
  Written 2011 Aug 10, by Leslie Young, SwRI
  2016 Mar 24 LAY. Modified for inclusion in vt3d library.

(See vt3d/vt3d_skindepth.pro)


VT3D_SOLAR_MU

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vt3d_solar_mu
 PURPOSE: (one line)
  Return cosine of the solar incidence angle, mu >= 0
 DESCRIPTION:
  Return cosine of the solar incidence angle, mu >= 0
 CATEGORY:
  Volatile Transport
 CALLING SEQUENCE:
  mu0 = vt3d_solar_mu(lat, lon, lat0, lon0)
 INPUTS:
  lat - scalar or array of latitudes (radian)
  lon - scalar or array of longitudes (radian)
  lat0 - scalar or array of sub-solar latitude (radian)
  lon0 - scalar or array of sub-solar longitude (radian)
 KEYWORD INPUT PARAMETERS:
  lonavg - return the longitudinal average for the given lat and lat0
  coslat, sinlat, coslon, sinlon : can pass these if they've
  been already computed
 OUTPUTS:
  mu0 = cos(theta0), where theta0 is angle of the sun with respect to the zenith
  If 
 COMMON BLOCKS:
  None
 SIDE EFFECTS:
 RESTRICTIONS:
  lat and lon must be scalars or both 1-D arrays of length n_loc
  lat0 and lon0 must be scalars or both 1-D arrays of length n_time
 PROCEDURE:
  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
  Resubmitted to Icarus.
  Eq. 3.1-5 (Eq 3.2-4a for longitudinal averaged mu0)
 MODIFICATION HISTORY:
  Written 2011 March 6, by Leslie Young, SwRI
  2016 Mar 22 LAY. Modified for inclusion in vt3d library.

(See vt3d/vt3d_solar_mu.pro)


VT3D_SOLWAVE

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vt3d_solwave
 PURPOSE: (one line)
  Return solar wave
 DESCRIPTION:
  Return solar wave
 CATEGORY:
  Volatile Transport
 CALLING SEQUENCE:
  sol = vt3d_solwave(sol_terms, phase)
 INPUTS:
  sol_terms = [sol_0, sol_1, ...], erg cm^-2 s^-1 : complex[n_term+1] or complex[n_loc,n_term+1]
   s.t., at the surface, sol(phase) = sol_0 + sol_1 exp(i*phase) + temp_2 exp(2*i*phase) + ...
  phase: float or float[n_t], where n_t is the number of time steps. frequency * time, radian
 OUTPUTS:
  sol: solar flux, erg cm^-2 s^-1. float, float[n_t], float[n_loc], or float[n_t,n_loc]

   sol_t           phase  sol  
   n_term+1        scalar scaler
   n_loc,n_term+1  scalar n_loc
   n_term+1        n_t    n_t
   n_loc,n_term+1  n_t    n_t,n_loc
   
 RESTRICTIONS:
 PROCEDURE:
  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
   Resubmitted to Icarus.
   Eq 3.2-1
 MODIFICATION HISTORY:
  Written 2012 Jan 11, by Leslie Young, SwRI
  2016 Mar 22 LAY. Modified for inclusion in vt3d library.

(See vt3d/vt3d_solwave.pro)


VT3D_SOL_TERMS_DIURNAL

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vt3d_sol_terms_diurnal
 PURPOSE: (one line)
  return the diurnal solar terms in sinusoidal expansion
 DESCRIPTION:
  return the diurnal solar terms such that
    s[i,lon] = sol_terms[i,0] + Re( sol_terms[i,1]*exp(_i*  p) )$
                          + Re( sol_terms[i,2]*exp(_i*2*p) ) + 
 CATEGORY:
  Volatile Transport
 CALLING SEQUENCE:
   sol_terms = vt3d_sol_terms_diurnal(dist_sol_au, albedo, lat, h_phase0, lat_sol, n_terms)
 INPUTS:
  dist_sol_au : sub-solar distance, in AU, float
  albedo : hemispheric albedo, unitless, float or float[n_loc]
  lat : latitude, radian, float or float[n_loc]
  h_phase0 : hour angle at zero phase
  lat_sol : subsolar latitude, radian, float
  n_terms : number of terms in the expansion
 OPTIONAL INPUTS:
  sol_norm_1au - Normal incident solar flux at 1 AU (scalar, erg/cm^2/s)
    Default 1367.6e3
 OUTPUTS:
  sol_terms : the complex diurnal solar terms in a sinusoidal expansion, 
    ergcm^2/s, complex[n_loc, n_term+1] or complex[n_term+1]
 OPTIONAL OUTPUTS:
  h_max - Hour angle at sunset, radian, float or float[n_loc]
 PROCEDURE:
  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
  Resubmitted to Icarus.
  Eq 3.2-4a,  3.2-4b,  3.2-4c
 MODIFICATION HISTORY:
  Written 2012 Jan 8, by Leslie Young, SwRI
  2012 Apr 19 LAY, change offset to h_phase0
  2016 Mar 22 LAY. Modified for inclusion in vt3d library.

(See vt3d/vt3d_sol_terms_diurnal.pro)


VT3D_STEP_CN_1LOC

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vt3d_step_cn_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_cn_1loc, 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
  beta_i:  dependance of temp[1..J] on temp[2..J+1], unitless, float[J]
  gamma_0: net flux at upper boundary, K, float
  gamma_J: net flux at lower boundary, K, float
 INPUT-OUTPUTS:
  temp_0 : temperature at layer 0, K, float
  temp_i : temperature at layer 1..J, K, float[J]
 PROCEDURE:
  Solve Crank-Nicholson equation
  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
   Resubmitted to Icarus.
   Eq. 3.4-9a, 3.4-9b
 MODIFICATION HISTORY:
  Written 2011 Feb 6, by Leslie Young, SwRI
  2016 Mar 24 LAY. Modified for inclusion in vt3d library.

(See vt3d/vt3d_step_cn_1loc.pro)


VT3D_STEP_CN_NLOC

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vt3d_step_cn_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_cnl_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:
  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
   Resubmitted to Icarus.
   Eq. 3.4-18
 MODIFICATION HISTORY:
  Written 2011 Feb 6, by Leslie Young, SwRI

(See vt3d/vt3d_step_cn_nloc.pro)


VT3D_STEP_CN_TRISOL_1LOC

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vt3d_step_cn_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_cn_trisol_1loc, alpha_i, beta_0, beta_i, gamma_0, gamma_J, spp_tridc, spp_indx, 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
  beta_i:  dependance of temp[1..J] on temp[2..J+1], unitless, float[J]
  gamma_0: net flux at upper boundary, K, float
  gamma_J: net flux at lower boundary, K, float
 INPUT-OUTPUTS:
  temp_0 : temperature at layer 0, K, float
  temp_i : temperature at layer 1..J, K, float[J]
 PROCEDURE:
  Solve Crank-Nicholson equation
  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
   Resubmitted to Icarus.
   Eq. 3.4-9a, 3.4-9b
 MODIFICATION HISTORY:
  Written 2011 Feb 6, by Leslie Young, SwRI
  2016 Mar 24 LAY. Modified for inclusion in vt3d library.

(See vt3d/vt3d_step_cn_trisol_1loc.pro)


VT3D_STEP_CN_TRISOL_NLOC

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vt3d_step_cn_trisol_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_cn_trisol_nloc, alpha_i, beta_0, beta_i, gamma_0, gamma_J, spp_tridc, spp_index, 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]
  spp_tridc : LU decomposition of the S'' ("spp") matrix, made by vt3d_cn_tridc
  spp_tridc : row permutations for partial pivoting to invert S'' ("spp"), made by vt3d_cn_tridc
  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:
  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
   Resubmitted to Icarus.
   Eq. 3.4-18
 MODIFICATION HISTORY:
  Written 2011 Feb 6, by Leslie Young, SwRI

(See vt3d/vt3d_step_cn_trisol_nloc.pro)


VT3D_STEP_EXPL_NLOC

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vt3d_step_expl_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_expl_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:
  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
   Resubmitted to Icarus.
   Eq. 3.4-8a, 3.4-8b
 MODIFICATION HISTORY:
  Written 2011 Feb 6, by Leslie Young, SwRI

(See vt3d/vt3d_step_expl_nloc.pro)


VT3D_STEP_IMPL_1LOC

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vt3d_step_impl_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_1loc, 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
  beta_i:  dependance of temp[1..J] on temp[2..J+1], unitless, float[J]
  gamma_0: net flux at upper boundary, K, float
  gamma_J: net flux at lower boundary, K, float
 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.  See:
   Eq. 3.4-9b
 MODIFICATION HISTORY:
  Written 2011 Feb 6, by Leslie Young, SwRI
  2016 Mar 24 LAY. Modified for inclusion in vt3d library.

(See vt3d/vt3d_step_impl_1loc.pro)


VT3D_STEP_IMPL_NLOC

[Previous Routine] [Next Routine] [List of Routines]
 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:
  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
   Resubmitted to Icarus.
   Eq. 3.4-18
 MODIFICATION HISTORY:
  Written 2011 Feb 6, by Leslie Young, SwRI

(See vt3d/vt3d_step_impl_nloc.pro)


VT3D_STEP_IMPL_TRISOL_1LOC

[Previous Routine] [Next Routine] [List of Routines]
 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.

(See vt3d/vt3d_step_impl_trisol_1loc.pro)


VT3D_STEP_IMPL_TRISOL_NLOC

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vt3d_step_impl_trisol_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_trisol_nloc, alpha_i, beta_0, gamma_0, gamma_J, spp_tridc, spp_index, 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]
  gamma_0: net flux at upper boundary, K, float, float[n_loc]
  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[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:
  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
   Resubmitted to Icarus.
   Eq. 3.4-18
 MODIFICATION HISTORY:
  Written 2011 Feb 6, by Leslie Young, SwRI

(See vt3d/vt3d_step_impl_trisol_nloc.pro)


VT3D_TEMP_TERM0_BARE

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vt3d_temp_term0_bare
 PURPOSE: (one line)
  Return temperature in balance with mean absorbed solar flux
 DESCRIPTION:
  Return temperature in balance with mean absorbed solar flux
 CATEGORY:
  Volatile Transport
 CALLING SEQUENCE:
  temp_0 = vt3d_temp_term0_bare(sol_0, flux_int, emis)
 INPUTS:
  sol_0 : mean absorbed solar flux, erg/cm^2/s : float[n_term+1] or float[n_loc,n_term+1]
  flux_int : internal heat flux, erg/cm^2/s, float or float[n_loc]
  emis : thermal emissivity, unitless, float or float[n_loc]
 OUTPUTS:
  temp_0: temperature in balance with mean absorbed solar flux, float or float[n_loc]
 RESTRICTIONS:
 PROCEDURE:
  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
   Resubmitted to Icarus.
   Eq 3.2-8
 MODIFICATION HISTORY:
  Written 2011 Dec 31, by Leslie Young, SwRI
  2016 Mar 24 LAY. Modified for inclusion in vt3d library.

(See vt3d/vt3d_temp_term0_bare.pro)


VT3D_TEMP_TERM0_LOCAL

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vt3d_temp_term0_local
 PURPOSE: (one line)
  Return temperature in balance with mean absorbed solar flux
 DESCRIPTION:
  Return temperature in balance with mean absorbed solar flux
 CATEGORY:
  Volatile Transport
 CALLING SEQUENCE:
  temp_0 = vt3d_temp_term0_local(sol_0, flux_int, emis)
 INPUTS:
  sol_0 : mean absorbed solar flux, erg/cm^2/s : float[n_term+1] or float[n_loc,n_term+1]
  flux_int : internal heat flux, erg/cm^2/s, float or float[n_loc]
  emis : thermal emissivity, unitless, float or float[n_loc]
  latheat : latent heat of sublimation, erg/g, loat or float[n_loc]
  mflux_esc : mass flux due to escape, g/cm^2/s, float or float[n_loc]
 OUTPUTS:
  temp_0: temperature in balance with mean absorbed solar flux, float or float[n_loc]
 RESTRICTIONS:
 PROCEDURE:
  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
   Resubmitted to Icarus.
   Eq 4.2-3
 MODIFICATION HISTORY:
  Written 2011 Dec 31, by Leslie Young, SwRI
  2016 Mar 24 LAY. Modified for inclusion in vt3d library.

(See vt3d/vt3d_temp_term0_local.pro)


VT3D_TEMP_TERMS_BARE

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vt3d_temp_terms_bare
 PURPOSE: (one line)
  Return temperature in balance with mean absorbed solar flux
 DESCRIPTION:
  Return temperature in balance with mean absorbed solar flux
 CATEGORY:
  Volatile Transport
 CALLING SEQUENCE:
  temp_terms = vt3d_temp_terms_bare(sol_terms,flux_int,emis,freq,therminertia)
 INPUTS:
  sol_terms : absorbed solar flux, complex or real, erg/cm^2/s:
   complex[n_term+1] or complex[n_loc,n_term+1]
     defined so that (if _i = sqrt(-1) )
     sol = sol_terms[0] + sol_terms[1]*exp(   _i * freq * t) 
                        + sol_terms[2]*exp(2 *_i * freq * t) + ..
     
  flux_int : internal heat flux, erg/cm^2/s: float or float[n_loc]
  emis : thermal emissivity, unitless : float or float[n_loc]
  freq : frequency of forcing
  therminertia : thermal inertia, erg cm-2 s-1/2 K-1 : float or float[n_loc]
 OUTPUTS:
  temp_terms: terms of approximate temperture, K, such that
     temp = temp_terms[0] + temp_terms[1]*exp(   _i * freq * t + sqrt(    _i * zeta) ) 
                      + temp_terms[2]*exp(2 *_i * freq * t + sqrt(2 * _i * zeta) ) + ..
  temp_terms is complex[n_term+1] or complex[n_loc,n_term+1]
 RESTRICTIONS:
 PROCEDURE:
  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
   Resubmitted to Icarus.
   Eq 3.2-10
 MODIFICATION HISTORY:
  Written 2011 Dec 31, by Leslie Young, SwRI
  2016 Mar 24 LAY. Modified for inclusion in vt3d library.

(See vt3d/vt3d_temp_terms_bare.pro)


VT3D_TEMP_TERMS_BARE_ITER

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vt3d_temp_terms_bare_iter
 PURPOSE: (one line)
  Return temperature in balance with mean absorbed solar flux
 DESCRIPTION:
  Return temperature in balance with mean absorbed solar flux
 CATEGORY:
  Volatile Transport
 CALLING SEQUENCE:
  temp_terms= vt3d_temp_terms_bare_iter(sol_terms,flux_int,emis,freq,therminertia, thermcond)
 INPUTS:
  sol_terms : absorbed solar flux, complex or real, erg/cm^2/s:
   complex[n_term+1] or complex[n_loc,n_term+1]
     sol = sol_terms[0] + sol_terms[1]*exp(_i * freq * t) 
                    + sol_terms[2]*exp(2 *_i freq t) + ..
     
  flux_int : internal hear flux, erg/cm^2/s: float or float[n_loc]
  emis : thermal emissivity, unitless : float or float[n_loc]
  freq : frequency of forcing
  therminertia : thermal inertia, erg cm-2 s-1/2 K-1 : float or float[n_loc]
  thermcond : thermal conductivity, erg K-1 : float or float[n_loc]
 OUTPUTS:
  temp_t: terms of approximate temperture, K, s.t.
     temp = temp_t[0] + temp_t[1]*exp(   _i * freq * t + sqrt(    _i * zeta) ) 
                      + temp_t[2]*exp(2 *_i * freq * t + sqrt(2 * _i * zeta) ) + ..
  temp_terms is complex[n_term+1] or complex[n_loc,n_term+1]
 RESTRICTIONS:
 PROCEDURE:
  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
   Resubmitted to Icarus.
   Eq. 3.2-10
 MODIFICATION HISTORY:
  Written 2011 Dec 31, by Leslie Young, SwRI
  2016 Mar 24 LAY. Modified for inclusion in vt3d library.

(See vt3d/vt3d_temp_terms_bare_iter.pro)


VT3D_TEMP_TERMS_LOCAL

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vt3d_temp_terms_local
 PURPOSE: (one line)
  Return temperature in balance with mean absorbed solar flux
 DESCRIPTION:
  Return temperature in balance with mean absorbed solar flux
 CATEGORY:
  Volatile Transport
 CALLING SEQUENCE:
  temp_terms = vt3d_temp_terms_local(sol_terms,flux_int,emis,freq,therminertia, $
                                 is_volatile, mass_volatile, specheat_volatile, gravacc, $
                                 name_species)
 INPUTS:
  sol_terms : absorbed solar flux, complex or real, erg/cm^2/s:
   complex[n_term+1] or complex[n_loc,n_term+1]
     defined so that (if _i = sqrt(-1) )
     sol = sol_terms[0] + sol_terms[1]*exp(   _i * freq * t) 
                        + sol_terms[2]*exp(2 *_i * freq * t) + ..
     
  flux_int : internal heat flux, erg/cm^2/s: float or float[n_loc]
  emis : thermal emissivity, unitless : float or float[n_loc]
  freq : frequency of forcing
  therminertia : thermal inertia, erg cm-2 s-1/2 K-1 : float or float[n_loc]
  is_volatile         ; 1 if has a volatile, 0 if bare. flag. scalar
  mass_slab           ; mass of the volatile slab, g cm^-2, scalar
  specheat_slab       ; specific heat of the volatile slab, erg g^-1 K^-1, scalar
  gravacc             ; gravitational acceleration, cm s^-2
  name_species        ; name of the species.  String. "N2", "CO", or "CH4"
 OPTIONAL INPUTS:
  mflux_esc_terms     ; escape rate in g cm^-2 s^-1
   complex[n_term+1] or complex[n_loc,n_term+1]
     defined so that (if _i = sqrt(-1) )
     mflux_esc = mflux_esc_terms[0] + mflux_esc_terms[1]*exp(   _i * freq * t) 
                                    + mflux_esc_terms[2]*exp(2 *_i * freq * t) + ..
 OUTPUTS:
  temp_terms: terms of approximate temperture, K, such that
     temp = temp_terms[0] + temp_terms[1]*exp(   _i * freq * t + sqrt(    _i * zeta) ) 
                      + temp_terms[2]*exp(2 *_i * freq * t + sqrt(2 * _i * zeta) ) + ..
  temp_terms is complex[n_term+1] or complex[n_loc,n_term+1]
 RESTRICTIONS:
  If defined, mflux_esc_terms must be the same dimensions as sol_terms
 PROCEDURE:
  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
   Resubmitted to Icarus.
   Eq 4.2-5
 MODIFICATION HISTORY:
  Written 2011 Dec 31, by Leslie Young, SwRI
  2016 Mar 24 LAY. Modified for inclusion in vt3d library.

(See vt3d/vt3d_temp_terms_local.pro)


VT3D_TEMP_TERMS_LOCAL_ITER

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vt3d_temp_terms_local_iter
 PURPOSE: (one line)
  Return temperature in balance with mean absorbed solar flux
 DESCRIPTION:
  Return temperature in balance with mean absorbed solar flux
 CATEGORY:
  Volatile Transport
 CALLING SEQUENCE:
  temp_terms = vt3d_temp_terms_local_iter(sol_terms,flux_int,emis,freq,therminertia, thermcond,  $
                                 is_volatile, mass_volatile, specheat_volatile, gravacc, name_species)
 INPUTS:
  sol_terms : absorbed solar flux, complex or real, erg/cm^2/s:
   complex[n_term+1] or complex[n_loc,n_term+1]
     sol = sol_terms[0] + sol_terms[1]*exp(_i * freq * t) 
                    + sol_terms[2]*exp(2 *_i freq t) + ..
     
  flux_int : internal hear flux, erg/cm^2/s: float or float[n_loc]
  emis : thermal emissivity, unitless : float or float[n_loc]
  freq : frequency of forcing
  therminertia : thermal inertia, erg cm-2 s-1/2 K-1 : float or float[n_loc]
  thermcond : thermal conductivity, erg K-1 : float or float[n_loc]
  is_volatile         ; 1 if has a volatile, 0 if bare. flag. scalar
  mass_slab           ; mass of the volatile slab, g cm^-2, scalar
  specheat_slab       ; specific heat of the volatile slab, erg g^-1 K^-1, scalar
  gravacc             ; gravitational acceleration, cm s^-2
  name_species        ; name of the species.  String. "N2", "CO", or "CH4"
  mflux_esc           ; escape rate in g cm^-2 s^-1
 OPTIONAL INPUTS:
  mflux_esc_terms     ; escape rate in g cm^-2 s^-1
   complex[n_term+1] or complex[n_loc,n_term+1]
     defined so that (if _i = sqrt(-1) )
     mflux_esc = mflux_esc_terms[0] + mflux_esc_terms[1]*exp(   _i * freq * t) 
                                    + mflux_esc_terms[2]*exp(2 *_i * freq * t) + ..
 OUTPUTS:
  temp_t: terms of approximate temperture, K, s.t.
     temp = temp_t[0] + temp_t[1]*exp(   _i * freq * t + sqrt(    _i * zeta) ) 
                      + temp_t[2]*exp(2 *_i * freq * t + sqrt(2 * _i * zeta) ) + ..
  temp_terms is complex[n_term+1] or complex[n_loc,n_term+1]
 RESTRICTIONS:
  If defined, mflux_esc_terms must be the same dimensions as sol_terms
 PROCEDURE:
  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
   Resubmitted to Icarus.
   Eq. 3.2-10
 MODIFICATION HISTORY:
  Written 2011 Dec 31, by Leslie Young, SwRI
  2016 Mar 24 LAY. Modified for inclusion in vt3d library.

(See vt3d/vt3d_temp_terms_local_iter.pro)


VT3D_THERMALINERTIA

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
 vt3d_thermalinertia
 PURPOSE: (one line)
  Return thermal inertia
 DESCRIPTION:
  Return thermal inertia
 CATEGORY:
  Volatile Transport
 CALLING SEQUENCE:
  therminertia = vt3d_thermalinertia(dens, specheat, thermcond)
 INPUTS:
  dens : density of the substrate (g/cm^3)
  specheat : specific heat (erg/(g K))
  thermcond : temperature of the substrate (erg/(cm s K))
 OUTPUTS:
  therminertia : thermal inertia, cgs units (erg cm^-2 s^-1/2 K^-1)
 RESTRICTIONS:
 PROCEDURE:
  Many, many places, including:
  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
   Resubmitted to Icarus.
   Eq 3.2-5
 MODIFICATION HISTORY:
  Written 2011 Aug 10, by Leslie Young, SwRI
  2016-04-28 Moved to vt3d library

(See vt3d/vt3d_thermalinertia.pro)


VT3D_THERMPRO

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vt3d_thermpro
 PURPOSE: (one line)
  Calculate temperatures on a surface in a style
  similar to thermprojrs
 DESCRIPTION:
  vt3d_thermpro
 CATEGORY:
  VT3D
 CALLING SEQUENCE:
  vt3d_thermpro,tsurf,tod,tdeep,teq,balance,tdarr,zarr0
 INPUTS:
  rhel heliocentric distance, AU
   alb bolometric albedo
   ti  thermal inertia, erg-cgs.  Can be an array, in which case it gives the ti of each slab
   rot Rotation period, days
 KEYWORDS: 
   emvty    emissivity                   (default 1.0)
   rho      density, g cm-3.  Can be an array, giving the density of each slab. (default 1.0)
   cp       specific heat, erg g-1 K-1 (default H20)
   plat     latitude (radians)           (default 0.0)
   sunlat   Subsolar latitude (radians)  (default 0.0)
   heatflow Endogenic heat flow, erg cm-2 s-1 (default 0.0)
   zarr     Array of depths to the BASE of each slab (NOT slab thicknesses), in cm.  
            First element must be non-zero.
   zeta     array to CENTER of each slab
   nslab    Number of slabs (ignored if zarr set)        (default 30)
   ntinc    Time increments per day                      (default 5000)
   nday     Number of days per run                       (default 4)
   skdepths Model depth, in skindepths (ignored if zarr set)  (default 6)
   sol_norm_1au solar constant, erg cm-2 s-1 
   silent   If set, sends no output to the terminal (except stability warnings)
 OUTPUT
   tsurf:    diurnal surface temperature array
   tod:      time of day array (in radians)
   tdeep:    final deep temperature
   teq:      Instantaneous-equilibrium temperature array
   balance:  Thermal balance: (power out - heatflow)/(power in)
   tdarr:    temperature/depth array, with 200 diurnal increments
             (set this to a defined variable before entering thermprojrs).
   zarr0:    for convenience, returns the array of slab-bottom depths
              (same as zarr)

 MODIFICATION HISTORY:
  Written 2012 Oct 12, by Leslie Young, SwRI
  2016 Mar 24 LAY. Modified for inclusion in layoung library.

(See vt3d/vt3d_thermpro.pro)


VT3D_THERMWAVE

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vt3d_thermwave
 PURPOSE: (one line)
  Return thermal wave
 DESCRIPTION:
  Return thermal wave
 CATEGORY:
  Volatile Transport
 CALLING SEQUENCE:
  temp = vt3d_thermwave(temp_terms, phase, dtemp_dzeta, z_skin, zeta)
 INPUTS:
  temp_terms = [temp_0, temp_1, ...], K : complex[n_term+1] or complex[n_loc,n_term+1]
   s.t., at the surface, temp(phase) = temp_0 + temp_1 exp(i*phase) + temp_2 exp(2*i*phase) + ...
  phase: float or float[n_t]. freq * time. Radian
  dtemp_dzeta: float, float[n_loc]
   - internal heat flux / (sqrt(freq) * therminertia)
  z_skin: float or float[n_z] or float[n_loc,n_z]. 
   If z_skin is an array or matrix, then the skin depth varies with depth, and
   this attempts to modify the wave equation with the WKB approximation.
   This has no effect if z_skin is a scalar (constant with depth), so
   it is common to simply set z_skin = 1.
  zeta: scaled depth, float, float[n_z], float[n_loc, n_z]
 OUTPUTS:
  temp: float, float[n_t], float[n_z], or float[n_t,n_z]
        float[n_t,n_loc], 

   temp_terms      phase  substrate  temp  
   n_term+1        scalar scalar     scaler
   n_loc,n_term+1  scalar scalar     n_loc
   n_term+1        n_t    scalar     n_t
   n_loc,n_term+1  n_t    scalar     n_t,n_loc
   n_term+1        scalar n_z        n_z
   n_loc,n_term+1  scalar n_z        n_loc,n_z
   n_term+1        n_t    n_z        n_t,n_z
   n_loc,n_term+1  n_t    n_z        n_t,n_loc,n_z
   n_term+1        scalar n_loc,n_z  n_loc,n_z
   n_loc,n_term+1  scalar n_loc,n_z  n_loc,n_z
   n_term+1        n_t    n_loc,n_z  n_t,n_loc,n_z
   n_loc,n_term+1  n_t    n_loc,n_z  n_t,n_loc,n_z
   
 RESTRICTIONS:
 PROCEDURE:
  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
   Resubmitted to Icarus.
   Eq 3.2-7
  If skin depth varies with z, then the wave is calculated with a wkb approximation.
 MODIFICATION HISTORY:
  Written 2011 Dec 30, by Leslie Young, SwRI
  Modified 2012 Sep 30, LAY.  Added m inside sqrt(_i) term.
  2016 Mar 22 LAY. Modified for inclusion in vt3d library.

(See vt3d/vt3d_thermwave.pro)


VT3D_THERMWAVE_1LOC_P0_NZ

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vt3d_thermwave_1loc_p0_nz
 PURPOSE: (one line)
  Return thermal wave, one location, phase=0, n depths
 DESCRIPTION:
  Return thermal wave, one location, phase=0, n depths
 CATEGORY:
  Volatile Transport
 CALLING SEQUENCE:
  temp = vt3d_thermwave_1loc_p0_nz(temp_terms, dtemp_dzeta, zeta)
 INPUTS:
  temp_terms : complex[n_terms + 1]
  phase : float[n_t] = 0.
  dtemp_dzeta : scalar
  zeta : float[n_z]
 OUTPUTS:
  temp: float[n_z]
   
 RESTRICTIONS:
 PROCEDURE:
 MODIFICATION HISTORY:
  Written 2011 Dec 30, by Leslie Young, SwRI
  Modified 2012 Sep 30, LAY.  Added m inside sqrt(_i) term.

(See vt3d/vt3d_thermwave_1loc_p0_nz.pro)


VT3D_THERMWAVE_NLOC_P0_NZ

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vt3d_thermwave_nloc_p0_nz
 PURPOSE: (one line)
  Return thermal wave, multiple locations, phase=0, n depths
 DESCRIPTION:
  Return thermal wave, multiple locations, phase=0, n depths
 CATEGORY:
  Volatile Transport
 CALLING SEQUENCE:
  temp = vt3d_thermwave_nloc_p0_nz(temp_terms, dtemp_dzeta, zeta)
 INPUTS:
  temp_terms : complex[nloc, n_terms + 1]
  phase : float[n_t] = 0.
  dtemp_dzeta : scalar
  zeta : float[n_z]
 OUTPUTS:
  temp: float[n_z]
   
 RESTRICTIONS:
 PROCEDURE:
 MODIFICATION HISTORY:
  Written 2011 Dec 30, by Leslie Young, SwRI
  Modified 2012 Sep 30, LAY.  Added m inside sqrt(_i) term.

(See vt3d/vt3d_thermwave_nloc_p0_nz.pro)


VT3D_VTSTEP_EXPL_1LOC

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vt3d_vtstep_expl_1loc
 PURPOSE: (one line)
  step the temperature one time step
 DESCRIPTION:
  step the temperature one time step
 CATEGORY:
  Volatile Transport
 CALLING SEQUENCE:
  vt3d_step_expl_1loc, alpha_i, beta_0, beta_i, gamma_0, gamma_J, temp_0, temp_i
 INPUTS:
  alpha_i: dependance of new temp[1..J] on old temp[0..J-1], unitless, float[J]
  beta_0:  dependance of new temp[0] on old temp[1], unitless, float
  beta_i:  dependance of new temp[1..J] on old temp[2..J+1], unitless, float[J]
  gamma_0: net flux at upper boundary, K, float
  gamma_J: net flux at lower boundary, K, float
 INPUT-OUTPUTS:
  temp_0 : temperature at layer 0, K, float. 
  temp_i : temperature at layer 1..J, K, float[J]
  Old temperatures are passed, new temperatures are returned.
 COMMON BLOCKS:
  None
 SIDE EFFECTS:
  temp_0 and temp_i are modified
 RESTRICTIONS:
  nsurf = 0 or 1
 PROCEDURE:
  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
   Resubmitted to Icarus.
   Eq. 3.4-2
 MODIFICATION HISTORY:
  Written 2011 Feb 6, by Leslie Young, SwRI
  2016 Mar 24 LAY. Modified for inclusion in vt3d library.

(See vt3d/vt3d_step_expl_1loc.pro)


VT3D_ZDELTA

[Previous Routine] [Next Routine] [List of Routines]
 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.

(See vt3d/vt3d_zdelta.pro)


VT3D_ZGRID

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vt3d_zgrid
 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:
  z = vt3d_zgrid(skindepth,z_delta,n_z,nsurf=nsurf,fsurf=fsurf,ntop=ntop, $
    dtop=dtop,nbot=nbot,dbot=dbot,fbot=fbot,hp96=hp96,sm92=sm92)
 INPUTS:
  skindepth - cm. A scalar or 1-D array of length n_loc
              skin depth = sqrt(k/(rho c omega)), where
              k = thermal conductivity
              rho = density
              c = specific heat
              omega = 2 pi/period = frequency
 OPTIONAL INPUT PARAMETERS:
  none
 KEYWORD INPUT PARAMETERS:
  nsurf = 1 or 0.  Number of "surface" layers, with the
          grid point defined at the top of the layer (z=0).
          DEFAULT 1
  fsurf = if positive, there is a surface layer with width
          fsurf * dtop * skindepth.  DEFAULT 0.5
  ntop =  number of layers in the top region DEFAULT 4
  dtop =  width of layers in the top region, normalized
          by the skin depth (so that width = dtop * skindepth)
          DEFAULT 0.25
  nbot =  number of layers in the bottom region DEFAULT 20
  dbot =  width of the first layer in the bottom region, normalized.
          DEFAULT fbot * dtop
  fbot =  fractional growth of the width of successive layers
          DEFAULT 1.13
  hp96 -  set all parameters to mimic Hansen and Paige 1996
          nsurf=0, fsurf=0.0, ntop=3, dtop=0.25, nbot=57, fbot=1.13
  sm92 -  set all parameters to mimic Spencer and Moore 1992
          nsurf=1, fsurf=0.5, ntop=60, dtop=0.1, nbot=0, fbot=1.00
 KEYWORD OUTPUT PARAMETERS:
  none
 OUTPUTS:
  n_z - number of layers (nsurf + ntop + nbot)
  z - cm. depth of grid points (z <= 0)
        1-D, length n_z, if skindepgth is a scalar, 
        2-D, [n_loc, n_z] if skindepth is an array
  z_delta - cm. width of layers, same dim. as z
 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

(See vt3d/vt3d_zgrid.pro)


VTY16_FIG3_1

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vty16_fig3_1
 PURPOSE: (one line)
  Plot Young 2016, Fig3.1
 DESCRIPTION:
  Plot Young 2016, Fig3.1
  Young, L. A. 2016, Volatile transport on inhomogeneous surfaces: II. Numerical calculations (VT3D)
   Resubmitted to Icarus.
 CATEGORY:
  VT3D
 CALLING SEQUENCE:
  vty16_fig3_1, phase, flux_sol, flux_sol_1, flux_sol_7
 INPUTS:
  none
 OPTIONAL OUTPUTS:
  phase: phase of the period running from 0 to 2 pi radian
  flux_sol : absorbed solar flux, erg/cm^2/s
  flux_sol_1 : approximate absorbed solar flux, M=1, erg/cm^2/s
  flux_sol_7 : approximate absorbed solar flux, M=7, erg/cm^2/s
 SIDE EFFECTS:
  creates files 
  vty16_fig3_1.png - as in Young 2016
  vty16_fig3_1.txt 
 MODIFICATION HISTORY:
  Written 2012 Sep 29, by Leslie Young, SwRI
  2016 Mar 22 LAY. Modified for inclusion in vty16 library.

(See vty16/vty16_fig3_1.pro)


VTY16_FIG3_14

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vty16_fig3_14
 PURPOSE: (one line)
  Plot Young 2016, Fig3-14
 DESCRIPTION:
  Plot Young 2016, Fig3-14
 CATEGORY:
  VT3D
 CALLING SEQUENCE:
  vty16_fig3_14
 INPUTS:
  none
 OPTIONAL OUTPUTS:
  temp_0_map float[n_lon,n_lat,n_t]
   surface temperature as a function of longitude, latitude and 
   time (aka sub-solar longitude)
   latitudes at cell center are [-88, -84, ... 84, 88]
   longitudes at cell center are [2, 6, ... 354, 358]
   sub-solar longitudes are [0, -1, -2, ... -1079, -1080]
     (mod 360, this is [0, -1, -2, ... -359, 0]
   
 SIDE EFFECTS:
  creates the plot in the Young 2016.
   vty16_fig3_14_043.png, vty16_fig3_14_087.png, vty16_fig3_14_0167.png, vty16_fig3_14_cb.png
    vty16_fig3_14.fits - output of vt3s_thermpro
 MODIFICATION HISTORY:
  Written 2013 Jun 5, by Leslie Young, SwRI
  Based on ms_plutoseason/figs/__working_copy/_mimas.pro
  2016 Mar 24 LAY. Modified for inclusion in vty16 library.

(See vty16/vty16_fig3_14.pro)


VTY16_FIG3_2

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vty16_fig3_2
 PURPOSE: (one line)
  Plot Young 2016, Fig3.2
 DESCRIPTION:
  Plot Young 2016, Fig3.2
 CATEGORY:
  VT3D
 CALLING SEQUENCE:
  vty16_fig3_2
 INPUTS:
  none
 OPTIONAL OUTPUTS:
  zeta : scaled depth (zeta < 0), unitless, float[n_z]
  temp_phase0 : temperature as a function of depth for zero phase,
    when surface temperature equals its mean. K, float[n_z]
  temp_phase90 : temperature as a function of depth for 90 deg phase,
    when surface temperature equals its minimum. K, float[n_z]
  temp_envelope_max : envelope of temperature maxima, K, float[n_z]
  temp_envelope_min : envelope of temperature minima, K, float[n_z]
 SIDE EFFECTS:
  creates vty16_fig3_2.png and vty16_fig3_2.txt
 MODIFICATION HISTORY:
  Written 2012 Sep 30, by Leslie Young, SwRI
  2016 Mar 22 LAY. Modified for inclusion in vt3d library.

(See vty16/vty16_fig3_2.pro)


VTY16_FIG3_3

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vty16_fig3_3
 PURPOSE: (one line)
  Plot Young 2016, Fig3.3
 DESCRIPTION:
  Plot Young 2016, Fig3.3
 CATEGORY:
  VT3D
 CALLING SEQUENCE:
  vty16_fig3_3
 INPUTS:
  none
 OPTIONAL OUTPUTS:
  phase_num - phase for the numeric calculation, radian
  temp_num - surface temperature from the numeric calculation, K
  phase - phase for the approximations, radian
  temp_1_zeta0 - approximate surface temperature, M=1, K
  temp_7_zeta0 - approximate surface temperature, M=7, K
  temp_7_iter_zeta0 - approximate surface temperature, 
                      M=7 with iterated T0, K
 SIDE EFFECTS:
  creates vty16_fig3_3.png, the plot in the Young 2016.
  vty16_fig3_3_table1.txt - output of vt3s_thermpro
  vty16_fig3_3_table2.txt - approximations
  vty16_fig3_3_table3.txt - temperature terms
 MODIFICATION HISTORY:
  Written 2012 Sep 29, by Leslie Young, SwRI
  2016 Mar 24 LAY. Modified for inclusion in vty16 library.

(See vty16/vty16_fig3_3.pro)


VTY16_FIG3_7A

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vty16_fig3_7a
 PURPOSE: (one line)
  Plot Young 2016, Fig3.7a
 DESCRIPTION:
  Plot Young 2016, Fig3.7a
 CATEGORY:
  VT3D
 CALLING SEQUENCE:
  vty16_fig3_7a
 INPUTS:
  none
 OPTIONAL OUTPUTS:
 SIDE EFFECTS:
  creates vty16_fig3_7a.png
 MODIFICATION HISTORY:
  Written 2012 Sep 29, by Leslie Young, SwRI
  2016 Mar 24 LAY. Modified for inclusion in vty16 library.

(See vty16/vty16_fig3_7a.pro)


VTY16_FIG3_7B

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vty16_fig3_7b
 PURPOSE: (one line)
  Plot Young 2016, Fig3.7b
 DESCRIPTION:
  Plot Young 2016, Fig3.7b
 CATEGORY:
  VT3D
 CALLING SEQUENCE:
  vty16_fig3_7b
 INPUTS:
  none
 OPTIONAL OUTPUTS:
 SIDE EFFECTS:
  creates vty16_fig3_7b.png
 MODIFICATION HISTORY:
  Written 2012 Sep 29, by Leslie Young, SwRI
  2016 Mar 24 LAY. Modified for inclusion in vty16 library.

(See vty16/vty16_fig3_7b.pro)


VTY16_FIG3_7C

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vty16_fig3_7c
 PURPOSE: (one line)
  Plot Young 2016, Fig3.7c
 DESCRIPTION:
  Plot Young 2016, Fig3.7c
 CATEGORY:
  VT3D
 CALLING SEQUENCE:
  vty16_fig3_7c
 INPUTS:
  none
 OPTIONAL OUTPUTS:
 SIDE EFFECTS:
  creates vty16_fig3_7c.png
 MODIFICATION HISTORY:
  Written 2012 Sep 29, by Leslie Young, SwRI
  2016 Mar 24 LAY. Modified for inclusion in vty16 library.

(See vty16/vty16_fig3_7c.pro)


VTY16_FIG4_2

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vty16_fig4_2
 PURPOSE: (one line)
  Plot Young 2016, Fig4-2
 DESCRIPTION:
  Plot Young 2016, Fig4-2
 CATEGORY:
  VT3D
 CALLING SEQUENCE:
  vty16_fig4_2
 INPUTS:
  none
 OPTIONAL OUTPUTS:
  
 SIDE EFFECTS:
  creates vty16_fig_14.png, the plot in the Young 2016.
  vty16_fig3_14_temp_surf.fits - output of vt3s_thermpro
 MODIFICATION HISTORY:
  Written 2013 Jun 5, by Leslie Young, SwRI
  Based on ms_plutoseason/figs/__working_copy/_mimas.pro
  2016 Mar 24 LAY. Modified for inclusion in vty16 library.

(See vty16/vty16_fig4_2.pro)


VTY16_FIG4_2_PLOT

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vty16_fig4_2_plot
 PURPOSE: (one line)
  Plot Young 2016, Fig4-2
 DESCRIPTION:
  Plot Young 2016, Fig4-2
 CATEGORY:
  VT3D
 CALLING SEQUENCE:
  vty16_fig4_2_plot, dist_sol_au, tod, theta_sva, temp_term_0, temp, p, $
                       name_therminertia, n_per_period
 INPUTS:
  none
 OPTIONAL OUTPUTS:
  
 SIDE EFFECTS:
  creates vty16_fig_14.png, the plot in the Young 2016.
  vty16_fig3_14_temp_surf.fits - output of vt3s_thermpro
 MODIFICATION HISTORY:
  Written 2013 Jun 5, by Leslie Young, SwRI
  Based on ms_plutoseason/figs/__working_copy/_mimas.pro
  2016 Mar 24 LAY. Modified for inclusion in vty16 library.

(See vty16/vty16_fig4_2_plot.pro)


VTY16_FIG5_7

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vty16_fig5_7
 PURPOSE: (one line)
  Plot Young 2016, Fig5.7
 DESCRIPTION:
  Plot Young 2016, Fig5.7
 CATEGORY:
  VT3D
 CALLING SEQUENCE:
  vty16_fig5_7
 INPUTS:
  none
 OPTIONAL OUTPUTS:
 SIDE EFFECTS:
  creates vty16_fig5_7.png and vty16_fig5_7.txt
 MODIFICATION HISTORY:
  Written 2012 Sep 29, by Leslie Young, SwRI
  2016 Mar 22 LAY. Modified for inclusion in vty16 library.

(See vty16/vty16_fig5_7.pro)


VTY16_FIG5_7_FUNC

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vty16_fig5_7_func
 PURPOSE: (one line)
  Make the data for a Pluto movie for Young 2013 ApJL resubmission
 DESCRIPTION:
  Make the data for a Pluto movie for Young 2013 ApJL resubmission
 CATEGORY:
  VT
 CALLING SEQUENCE:
  res = vty16_fig5_7_func(run, av, ev, as, es, ti, mvtot, n_off, res_all)
 INPUTS:
  run - name of this run.  If "", then construct from parameters
  av - hemispheric albedo of the volatile
  ev - emissivity of the volatile
  as - hemispheric albedo of the substrate
  es - emissivity of the substrate
  ti - thermal inertia of the substrate, cgs units 
  mvtot - total inventory of volatiles
  n_off - number of time periods to offset start by.
      If n_off = 0, start at aphelion.
      If n_off > 0, start n_off timesteps before aphelion
 OPTIONAL INPUT PARAMETERS:
  none
 KEYWORD INPUT PARAMETERS:
  none
 KEYWORD OUTPUT PARAMETERS:
  none
 OUTPUTS:
  res: structure with the results of the run, suitable for plotting.
    ALBEDO_COMP     FLOAT     Array[3]  albedo of composition types, unitless
    EMIS_COMP       FLOAT     Array[3]  albedo of composition types, unitless
       
    AV              FLOAT     albedo of volatile, unitless
    EV              FLOAT     emissivity of volatile, unitless
    AS              FLOAT     albedo of substrate, unitless
    ES              FLOAT     emissivity of substrate, unitless
    THERMINERTIA    DOUBLE    Array[60, 19] 
                    thermal inertia, [loc, depth], cgs
    MVTOT           INT             32 total volatile, g/cm^2
    LAT             DOUBLE    Array[60] latitude of locations, radian
    N_LAT_EDGE      LONG                59 number of latitude edges
    N_TIME_PER_PERIOD
                    LONG               240 timesteps per period
    N_TIME          LONG               720 number of timesteps
    A               FLOAT     Array[240] geometric albedo, [time] 
    TV              FLOAT     Array[240] temperature of volatile, K, [time]
    TS              FLOAT     Array[60, 240] temperature of layer 0, K, [loc,time]
    MS              FLOAT     Array[60, 240] 
                       mass of slab, g/cm^2, [loc, time]
    V               FLOAT     Array[59, 240] 
                       velocity across lat boundary, cm/s, [lat_edge]
    VS              FLOAT     Array[59, 240]
                       sound speed at lat boundary, cm/s, [lat_edge]
    P               FLOAT     Array[240] pressure, microbar, [time]
    YR              DOUBLE    Array[240] time, year, [time]
    LAT0            FLOAT     Array[240] subsolar lat, radian, [time]
    R               FLOAT     Array[240] helioc. dist, AU, [time]
    ANG             FLOAT     Array[240] true anomoly, radian, [time]
    ISV             INT       Array[60, 240] 
                       is volatile-covered, boolean, [loc, time]
  res_all: structure with the results of the run, suitable for full output.
    All the tags in res, PLUS
 COMMON BLOCKS:
  None
 SIDE EFFECTS:
 RESTRICTIONS:
  None
 PROCEDURE:
 MODIFICATION HISTORY:
  Written 2012 October 13, by Leslie Young, SwRI (pluto_dps12_2.pro)
  Oct 21 2012.  Make function (pluto_mssearch_func.pro)
  Dec 13 2012.  Add start phase as options
  2016 Apr 06 LAY. Modified for inclusion in vty16 library.

(See vty16/vty16_fig5_7_func.pro)


VTY16_FIG5_7_MAT

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vty16_fig5_7_mat
 PURPOSE: (one line)
  Calculate matrix elements for Young 2013 ApJL resubmission
 DESCRIPTION:
  Calculate matrix elements that change with time for Young 2013 ApJL resubmission
 CATEGORY:
  VT
 CALLING SEQUENCE:
  vty16_fig5_7_mat, flag_frostslab, time_delta, n_loc, n_z, $
   emis, temp_surf, eflux_sol, mass_slab, specheat_frost, $   
   z_delta, z_delta_bot, dens, specheat, thermcond, eflux_int, $
   beta,alpha_bot, alpha_mid,    
 INPUTS:
   flag_frostslab  ; byte. 1 to include frost slab in matrix calc
   time_delta      ; float. duration of timestep (s)
   n_loc           ; long. number of locations
   n_z             ; long, number of layers
   emis            ; float[n_loc], emissivity (uniteless)
   temp_surf       ; float[n_loc], temp at start of timestep (K)
   eflux_sol       ; float[n_loc], solar absorption, erg/cm^2/s
   mass_slab       ; float[n_loc], mass of volatile slab, g/cm^2
   specheat_frost  ; float[n_loc], specific heat of volatile slab erg/g/K
   z_delta         ; float[n_loc, n_z]. Width of layer. cm
   z_delta_bot     ; float[n_loc, n_z], Distance from this layer to lower layer. cm.
   dens            ; float[n_loc, n_z]. Density, g/cm^3
   specheat        ; float[n_loc, n_z]. specific heat of substrate erg/g/K
   thermcond       ; float[n_loc, n_z]. thermal conductivity of substrate erg/cm/s/K
   eflux_int       ; float[n_loc]. internal heat flux, erg/cm^2/s
 INPUT/OUTPUTS:
   beta            ; float[n_loc, n_z]. constant arrays in T_new = alpha * T + beta
     This is written as "gamma" in Young 2016.  Just to keep you
     alert. 
     Only the top layer (beta[*,0]) and bottom layer (beta[*,n_z-1]) are altered
   alpha_bot       ; float[n_loc, n_z]. upper tridiagonal of alpha in
     T_new = alpha * T + beta
     called "bot" because it relates the current temperature layer
     and the one below it.
     This is written as "beta" in Young 2016, just to keep you alert.
     Only the top layer (alpha_bot[*,0]) is altered
   alpha_mid       ; float[n_loc, n_z]. Related to the diagonal of
     T_new = alpha * T + beta
     alpha_mid is eta in Young 2016, just to keep you alert.
     Only the top layer (alpha_mid[*,0]) is altered
 OUTPUTS:
   denom           
   eflux_net
 COMMON BLOCKS:
  None
 SIDE EFFECTS:
 RESTRICTIONS:
  None
 PROCEDURE:
 MODIFICATION HISTORY:
  Written 2012 October 13, by Leslie Young, SwRI (pluto_dps12_2.pro)
  Oct 21 2012.  Make function (pluto_mssearch_func.pro)
  Dec 13 2012.  Add start phase as options
  2016 Apr 06 LAY. Modified for inclusion in vty16 library.

(See vty16/vty16_fig5_7_mat.pro)


VTY16_FIG5_7_WRITE

[Previous Routine] [Next Routine] [List of Routines]
 NAME:
  vty16_fig5_7_write
 PURPOSE: (one line)
  Write the data for a Pluto movie for Young 2013 ApJL resubmission
 DESCRIPTION:
  Write the data for a Pluto movie for Young 2013 ApJL resubmission
 CATEGORY:
  VT
 CALLING SEQUENCE:
  vty16_fig5_7_write, run, res_all
 INPUTS:
  run : name of this run.
  res: structure with the results of the run, suitable for plotting.
    See vty16_fig5_7_func
 SIDE EFFECTS:
    write files _.png
    where  is tag in the structure res_all
 RESTRICTIONS:
  None
 PROCEDURE:
 MODIFICATION HISTORY:
  Written 2012 October 13, by Leslie Young, SwRI (pluto_dps12_2.pro)
  Oct 21 2012.  Make function (pluto_mssearch_func.pro)
  Dec 13 2012.  Add start phase as options
  2016 Apr 06 LAY. Modified for inclusion in vty16 library.

(See vty16/vty16_fig5_7_write.pro)


VTY16_PLUTO_MSSEARCH_RESUB_FUNC

[Previous Routine] [List of Routines]
 NAME:
  vty16_pluto_mssearch_resub_func
 PURPOSE: (one line)
  Make the data for a Pluto movie for Young 2013 ApJL resubmission
 DESCRIPTION:
  Make the data for a Pluto movie for Young 2013 ApJL resubmission
 CATEGORY:
  VT
 CALLING SEQUENCE:
  vty16_fig5_7_still, run, res, yr_still
 INPUTS:
  run : name of this run.
  res: structure with the results of the run, suitable for plotting.
    See vty16_pluto_mssearch_resub_func
  yr_still : year for which to make the plot
 SIDE EFFECTS:
    write a file fig//mov/_.png
    where  is the time index relative to the start of the last period
 RESTRICTIONS:
  None
 PROCEDURE:
 MODIFICATION HISTORY:
  Written 2012 October 13, by Leslie Young, SwRI (pluto_dps12_2.pro)
  Oct 21 2012.  Make function (pluto_mssearch_func.pro)
  Dec 13 2012.  Add start phase as options
  2016 Apr 06 LAY. Modified for inclusion in vty16 library.

(See vty16/vty16_fig5_7_still.pro)