;+
; 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.
;-

function vt3d_dfluxdtemp_atm, freq, temp_v, frac_varea, gravacc, name_species, latheat, a

    physconstants

    case name_species of
       'N2' : begin
          p0 = n2vp(temp_v)
          latheat = n2latentheat(temp_v)
          mass_molec = molec_weight(22,1) * !phys.m_u
       end
       'CO' : begin
          p0 = covp(temp_v)
          latheat = colatentheat(temp_v)
          mass_molec = molec_weight(5,1) * !phys.m_u
       end
       'CH4' : begin
          p0 = ch4vp(temp_v)
          latheat = ch4latentheat(temp_v) ; erg/g
          mass_molec = molec_weight(6,1) * !phys.m_u ; g/molec
       end 
       else : begin
          p0 = 0.
          latheat = 0.
          return, 0.
       end
    endcase

    a = latheat * mass_molec / (!phys.k * temp_v)
    dpdt = a * p0 / temp_v
    
    return, freq * (1./frac_varea) * (dpdt) * (latheat / gravacc)
end