
## Does "land_surface_air__latent_heat_flux" make sense? (2/5/13)

# Copyright (c) 2001-2014, Scott D. Peckham
#
#  Sep 2014.  Fixed sign error in update_bulk_richardson_number().
#             Ability to compute separate P_snow and P_rain.
#  Aug 2014.  New CSDMS Standard Names and clean up.
#  Nov 2013.  Converted TopoFlow to a Python package.
#
#  Jan 2013. Revised handling of input/output names.
#  Oct 2012. CSDMS Standard Names (version 0.7.9) and BMI.
#  May 2012. P is now a 1D array with one element and mutable,
#            so any comp with ref to it can see it change.
#  Jun 2010. update_net_shortwave_radiation(), etc.  
#  May 2010. Changes to initialize() and read_cfg_file().
#  Aug 2009
#  Jan 2009. Converted from IDL.
#
#-----------------------------------------------------------------------
# NOTES: This file defines a "base class" for meteorology
#        components as well as any functions used by most or
#        all meteorology methods.  The methods of this class
#        should be over-ridden as necessary for different
#        methods of modeling meteorology.
#-----------------------------------------------------------------------
# Notes: Do we ever need to distinguish between a surface
#        temperature and snow temperature (in the snow) ?
#        Recall that a separate T_soil_x variable is used
#        to compute Qc.
#
#        Cp_snow is from NCAR CSM Flux Coupler web page
#
#        rho_H2O is currently not adjustable with GUI. (still true?)
#    
#-----------------------------------------------------------------------    
#
#  class met_component     (inherits from BMI_base.py)
#
#      get_component_name()
#      get_attribute()           # (10/26/11)
#      get_input_var_names()     # (5/15/12)
#      get_output_var_names()    # (5/15/12)
#      get_var_name()            # (5/15/12)
#      get_var_units()           # (5/15/12)
#      ---------------------
#      set_constants()
#      initialize()
#      update()
#      finalize()
#      ----------------------------
#      set_computed_input_vars()
#      initialize_computed_vars()
#      ----------------------------
#      update_P_integral()
#      update_P_max()
#      update_P_rain()    # (9/14/14, new method)
#      update_P_snow()    # (9/14/14, new method)
#      ------------------------------------
#      update_bulk_richardson_number()
#      update_bulk_aero_conductance()
#      update_sensible_heat_flux()
#      update_saturation_vapor_pressure()
#      update_vapor_pressure()
#      update_dew_point()                    # (7/6/10)
#      update_precipitable_water_content()   # (7/6/10)
#      ------------------------------------
#      update_latent_heat_flux()
#      update_conduction_heat_flux()
#      update_advection_heat_flux()
#      ------------------------------------
#      update_julian_day()                   # (7/1/10)
#      update_net_shortwave_radiation()      # (7/1/10)
#      update_em_air()                       # (7/1/10)
#      update_net_longwave_radiation()       # (7/1/10)
#      update_net_energy_flux()              # ("Q_sum")
#      ------------------------------------
#      open_input_files()
#      read_input_files()
#      close_input_files()
#      ------------------------------------
#      update_outfile_names()
#      open_output_files()
#      write_output_files()
#      close_output_files()
#      save_grids()
#      save_pixel_values()
#
#  Functions:
#      compare_em_air_methods()
#
#-----------------------------------------------------------------------

import numpy as np
import os

from topoflow.components import solar_funcs as solar

from topoflow.utils import BMI_base
from topoflow.utils import model_input
from topoflow.utils import model_output
from topoflow.utils import rtg_files

#-----------------------------------------------------------------------
class met_component( BMI_base.BMI_component ):

    #-------------------------------------------------------------------
    _att_map = {
        'model_name':         'TopoFlow_Meteorology',
        'version':            '3.1',
        'author_name':        'Scott D. Peckham',
        'grid_type':          'uniform',
        'time_step_type':     'fixed',
        'step_method':        'explicit',
        #-------------------------------------------------------------
        'comp_name':          'Meteorology',
        'model_family':       'TopoFlow',
        'cfg_template_file':  'Meteorology.cfg.in',
        'cfg_extension':      '_meteorology.cfg',
        'cmt_var_prefix':     '/Meteorology/Input/Var/',
        'gui_xml_file':       '/home/csdms/cca/topoflow/3.1/src/share/cmt/gui/Meteorology.xml',
        'dialog_title':       'Meteorology: Method 1 Parameters',
        'time_units':         'seconds' }

    #---------------------------------------------------------
    # Note that SWE = "snow water equivalent", but it really
    # just means "liquid_equivalent".
    #---------------------------------------------------------   
    _input_var_names = [
        'snowpack__z_mean_of_mass-per-volume_density', # rho_snow
        'snowpack__depth',                             # h_snow
        'snowpack__liquid-equivalent_depth',           # h_swe
        'snowpack__melt_volume_flux' ]                 # SM   (MR used for ice?)

    #-----------------------------------------------------------
    # albedo, emissivity and transmittance are dimensionless.
    #-----------------------------------------------------------
    # "atmosphere_aerosol_dust__reduction_of_transmittance" vs.
    # This TF parameter comes from Dingman, App. E, p. 604.
    #-----------------------------------------------------------
    # There is an Optical_Air_Mass function in solar_funcs.py.
    # However, this quantity is not saved in comp state.
    #
    # "optical_path_length_ratio" vs. "optical_air_mass" OR
    # "airmass_factor" OR "relative_airmass" OR
    # "relative_optical_path_length"
    #-----------------------------------------------------------
    # Our term "liquid_equivalent_precipitation" is widely
    # used on the Internet, with 374,000 Google hits.
    #--------------------------------------------------------------
    # Note: "bulk exchange coefficient" has 2460 Google hits.
    #       It is closely related to a "transfer coefficient"
    #       for mass, momentum or heat.  There are no CF
    #       Standard Names with "bulk", "exchange" or "transfer".
    #
    # Zhang et al. (2000) use "bulk exchange coefficient" in a
    # nonstandard way, with units of velocity vs. unitless.
    #
    # Dn = bulk exchange coeff for the conditions of
    #      neutral atmospheric stability [m/s]
    # Dh = bulk exchange coeff for heat  [m/s]
    # De = bulk exchange coeff for vapor [m/s]
    #---------------------------------------------------------------
    # Now this component uses T_air to break the liquid-equivalent
    # precip rate into separate P_rain and P_snow components.
    # P_rain is used by channel_base.update_R()
    # P_snow is used by snow_base.update_depth()
    #---------------------------------------------------------------
    _output_var_names = [
        # 'atmosphere__optical_path_length_ratio',                           # M_opt [1]  (in solar_funcs.py)
        # 'atmosphere__von_karman_constant',                                 # kappa
        'atmosphere_aerosol_dust__reduction_of_transmittance',               # dust_atten  ##### (from GUI)
        'atmosphere_air-column_water-vapor__liquid-equivalent_depth',        # W_p ("precipitable depth")
        'atmosphere_bottom_air__brutsaert_emissivity_canopy_factor',         # canopy_factor
        'atmosphere_bottom_air__brutsaert_emissivity_cloud_factor',          # cloud_factor
        'atmosphere_bottom_air__bulk_latent_heat_aerodynamic_conductance',   # De [m s-1], latent
        'atmosphere_bottom_air__bulk_sensible_heat_aerodynamic_conductance', # Dh [m s-1], sensible
        'atmosphere_bottom_air__emissivity',                                 # em_air
        'atmosphere_bottom_air__mass-per-volume_density',                    # rho_air
        'atmosphere_bottom_air__mass-specific_isobaric_heat_capacity',       # Cp_air
        'atmosphere_bottom_air__neutral_bulk_aerodynamic_conductance',       # Dn [m s-1], neutral
        'atmosphere_bottom_air__pressure',                                   # p0
        'atmosphere_bottom_air__temperature',                                # T_air
        'atmosphere_bottom_air_flow__bulk_richardson_number',                # Ri [1]
        'atmosphere_bottom_air_flow__log_law_roughness_length',              # z0_air
        'atmosphere_bottom_air_flow__reference-height_speed',                # uz
        'atmosphere_bottom_air_flow__speed_reference_height',                # z
        'atmosphere_bottom_air_land_net-latent-heat__energy_flux',           # Qe [W m-2]
        'atmosphere_bottom_air_land_net-sensible-heat__energy_flux',         # Qh [W m-2]
        'atmosphere_bottom_air_water-vapor__dew_point_temperature',          # T_dew
        'atmosphere_bottom_air_water-vapor__partial_pressure',               # e_air # (insert "reference_height" ??)
        'atmosphere_bottom_air_water-vapor__relative_saturation',            # RH
        'atmosphere_bottom_air_water-vapor__saturated_partial_pressure',     # e_sat_air         
        'atmosphere_water__domain_time_integral_of_precipitation_leq-volume_flux',  # vol_P
        'atmosphere_water__domain_time_max_of_precipitation_leq-volume_flux',     # P_max
        'atmosphere_water__precipitation_leq-volume_flux',                        # P [m s-1]
        'atmosphere_water__rainfall_volume_flux',            # P_rain [m s-1] (liquid)      
        'atmosphere_water__snowfall_leq-volume_flux',        # P_snow [m s-1]
        'earth__standard_gravity_constant',                  # g   [m s-2]
        'land_surface__albedo',                              # albedo
        'land_surface__aspect_angle',                        # alpha  (from GUI)
        'land_surface__emissivity',                          # em_surf
        'land_surface__latitude',                            # lat_deg [degrees]
        'land_surface__longitude',                           # lon_deg [degrees]
        'land_surface__slope_angle',                         # beta  (from GUI)
        'land_surface__temperature',                         # T_surf   ### OR JUST "land__temperature"?
        # 'land_surface_air__temperature',                          # T_air
        'land_surface_air_water-vapor__partial_pressure',           # e_surf # (insert "reference_height" ??)
        'land_surface_air_water-vapor__saturated_partial_pressure', # e_sat_surf
        'land_surface_net-longwave-radiation__energy_flux',  # Qn_LW [W m-2]
        'land_surface_net-shortwave-radiation__energy_flux', # Qn_SW [W m-2]
        'land_surface_net-total-energy__energy_flux',        # Q_sum [W w-2]
        'model__time_step',                                  # dt
        'physics__stefan_boltzmann_constant',                # sigma  [W m-2 K-4]
        'physics__von_karman_constant',                      # kappa  [1]
        'water__mass-specific_latent_fusion_heat',           # Lf     [J kg-1]
        'water__mass-specific_latent_vaporization_heat',     # Lv     [J kg-1]
        'water-liquid__mass-per-volume_density' ]            # rho_H2O
        
    #-----------------------------------------
    # These are used only in solar_funcs.py
    # Later, create a Radiation component.
    #---------------------------------------------
    # Should we allow "day" as a base quantity ?
    # "day_length" is confusing.  Think about "date" also.
    # Maybe something like:
    #
    #    "earth__mean_solar_rotation_period"
    #    "earth__sidereal_rotation_period"
    #    "earth__stellar_rotation_period"   (relative to "fixed stars")
    #         maybe:  "earth__complete_rotation_period" ??
    #
    #    OR:
    #    "earth_mean_solar_day__duration"
    #    "earth_sidereal_day__duration"
    #    "earth_stellar_day__duration"
    #
    #    OR perhaps:
    #    "earth_mean_solar_day__rotation_period"
    #    "earth_sidereal_day__rotation_period"
    #    "earth_stellar_day__rotation_period"
    #
    #    "stellar rotation period" gives 84,500 Google hits.
    #    "solar_rotation_period" gives 41,100 Google hits.
    #    "sidereal_roation_period" gives 86,000 Google hits.
    #    "stellar day" gives 136,000 Google hits (but many unrelated).
    #
    #    NB! "stellar_rotation_period" is ambiguous since it is also
    #         used for the rotation period of a star.
    #
    #    "earth_mean_solar_day__hour_count"  ("standard_day" ?)
    #    "earth_sidereal_day__hour_count"
    #    "earth_sidereal_day__duration"
    #    "earth__rotation_period"   = "sidereal_day"
    #
    #    "earth_stellar_day__period"  ??
    #    "earth_stellar_day__duration" ??
    #
    #------------------------------------------------------------------
    # For "earth__rotation_rate", it seems this should be based on
    # the sidereal day (23.93 hours) instead of the mean solar day.
    #------------------------------------------------------------------
    # There are at least a few online sources that use both terms:
    # "equivalent latitude" and "equivalent longitude".  See:
    # "The Construction and Application of a Martian Snowpack Model".
    #------------------------------------------------------------------
    # Adopt the little-used term:  "topographic_sunrise" ?
    # Or maybe "illuminated_topography", or "local_sunrise" ??
    #------------------------------------------------------------------
    # For angle relations between the earth and the sun, should we
    # just use the adjective "solar" in the quantity name or include
    # sun in the object name?  We could also use terms like:
    #     earth_to_sun__declination_angle
    #     earth_to_sun__right_ascension_angle
    #
    #------------------------------------------------------------------
    # The adjective "local" in "earth_local_apparent_noon__time"
    # may be helpful in other contexts such as:
    # 'earth__local_longitude' and 'land_surface__local_elevation'.
    #------------------------------------------------------------------    
    # 'earth__autumnal_equinox_date',
    # 'earth__autumnal_equinox_time',
    # 'earth_axis__ecliptic_tilt_angle',  # tilt_angle
    # 'earth__julian_day_number',         ########
    # 'earth__julian_day_angle',
    # 'earth__local_apparent_noon_time'
    # 'earth__mean_radius', 
    # 'earth__mean_solar_day_duration',   # (exactly 24 hours)
    # 'earth_orbit__eccentricity',
    # 'earth_orbit__period',     # (one year)
    # 'earth__perihelion_julian_day',  ######
    # 'earth__rotation_period',        ######
    # 'earth__rotation_rate', # Omega       ###### What about Angular Velocity ?
    # 'earth__sidereal_day_duration',     # (one rotation = 23.934470 hours)
    # 'earth__solar_declination_angle',
    # 'earth__solar_hour_angle',
    # 'earth__solar_irradiation_constant',  ## (or "insolation_constant" ??)
    # 'earth__solar_right_ascension_angle',
    # 'earth__solar_vertical_angle',   (complement of zenith angle)
    # 'earth__solar_zenith_angle',
    # 'earth__stellar_day_duration',  # (relative to the "fixed stars")
    # 'earth__summer_solstice_date',
    # 'earth__summer_solstice_time',
    # 'earth__topographic_sunrise_equivalent_latitude',
    # 'earth__topographic_sunrise_equivalent_longitude',  (flat_lon + offset) 
    # 'earth__topographic_sunrise_equivalent_longitude_offset',
    # 'earth__topographic_sunrise_time',
    # 'earth__topographic_sunset_time',
    # 'earth_true_solar_noon___time',  #####
    #     'earth_clock__true_solar_noon_time'
    # 'earth__vernal_equinox_date',
    # 'earth__vernal_equinox_time',
    # 'earth__winter_solstice_date',
    # 'earth__winter_solstice_time',
    #
    # What about a "slope_corrected" or "topographic" version of K_dir ?
    #
    # 'land_surface__backscattered_shortwave_irradiation_flux', # K_bs
    # 'land_surface__diffuse_shortwave_irradiation_flux',       # K_dif
    # 'land_surface__direct_shortwave_irradiation_flux',        # K_dir
    # 'land_surface__global_shortwave_irradiation_flux',        # K_glob = K_dif + K_dir
    #------------------------------------------------------------------


    #------------------------------------------------------------------   
    # Maybe we should rename "z" to "z_ref" and "uz" to "uz_ref" ?
    #------------------------------------------------------------------   
    _var_name_map = {
        'snowpack__z_mean_of_mass-per-volume_density': 'rho_snow',
        'snowpack__depth': 'h_snow',
        'snowpack__liquid-equivalent_depth': 'h_swe',
        'snowpack__melt_volume_flux': 'SM',              # (MR is used for ice)
        #-----------------------------------------------------------------
        #'atmosphere__optical_path_length_ratio': 'M_opt',    # (in solar_funcs.py)
        # 'atmosphere__von_karman_constant': 'kappa',
        'atmosphere_aerosol_dust__reduction_of_transmittance': 'dust_atten',
        'atmosphere_air-column_water-vapor__liquid-equivalent_depth': 'W_p',   #########
        'atmosphere_bottom_air__brutsaert_emissivity_canopy_factor': 'canopy_factor',
        'atmosphere_bottom_air__brutsaert_emissivity_cloud_factor':  'cloud_factor',
        'atmosphere_bottom_air__bulk_latent_heat_aerodynamic_conductance':   'De',
        'atmosphere_bottom_air__bulk_sensible_heat_aerodynamic_conductance': 'Dh',
        'atmosphere_bottom_air__emissivity': 'em_air',               
        'atmosphere_bottom_air__mass-per-volume_density': 'rho_air',
        'atmosphere_bottom_air__mass-specific_isobaric_heat_capacity': 'Cp_air',
        'atmosphere_bottom_air__neutral_bulk_heat_aerodynamic_conductance':  'Dn',
        'atmosphere_bottom_air__pressure':                           'p0',
        'atmosphere_bottom_air__temperature':                        'T_air',
        'atmosphere_bottom_air_flow__bulk_richardson_number':        'Ri',        
        'atmosphere_bottom_air_flow__log_law_roughness_length':      'z0_air', ## (not "z0")
        'atmosphere_bottom_air_flow__reference-height_speed':        'uz',
        'atmosphere_bottom_air_flow__speed_reference_height':        'z',
        'atmosphere_bottom_air_land_net-latent-heat__energy_flux':   'Qe',
        'atmosphere_bottom_air_land_net-sensible-heat__energy_flux': 'Qh',
        'atmosphere_bottom_air_water-vapor__dew_point_temperature':  'T_dew',  
        'atmosphere_bottom_air_water-vapor__partial_pressure':       'e_air',
        'atmosphere_bottom_air_water-vapor__relative_saturation':    'RH',
        'atmosphere_bottom_air_water-vapor__saturated_partial_pressure': 'e_sat_air',
        'atmosphere_water__domain_time_integral_of_precipitation_leq-volume_flux': 'vol_P', 
        'atmosphere_water__domain_time_max_of_precipitation_leq-volume_flux': 'P_max',       
        'atmosphere_water__precipitation_leq-volume_flux': 'P',
        'atmosphere_water__rainfall_volume_flux':          'P_rain',    
        'atmosphere_water__snowfall_leq-volume_flux':      'P_snow',
        'earth__standard_gravity_constant':                'g',
        'land_surface__albedo':                            'albedo',
        'land_surface__aspect_angle':                      'alpha',
        'land_surface__emissivity':                        'em_surf',
        'land_surface__latitude':                          'lat_deg',
        'land_surface__longitude':                         'lon_deg',
        'land_surface__slope_angle':                       'beta',
        'land_surface__temperature':                       'T_surf',
         # 'land_surface_air__temperature': 'T_surf',
        'land_surface_air_water-vapor__partial_pressure':           'e_surf',
        'land_surface_air_water-vapor__saturated_partial_pressure': 'e_sat_surf',        
        'land_surface_net-longwave-radiation__energy_flux':         'Qn_LW',
        'land_surface_net-shortwave-radiation__energy_flux':        'Qn_SW',
        'land_surface_net-total-energy__energy_flux':               'Q_sum',
        'model__time_step':                                         'dt',
        'physics__stefan_boltzmann_constant':                       'sigma',
        'physics__von_karman_constant':                             'kappa',
        'water__mass-specific_latent_fusion_heat':                  'Lf',
        'water__mass-specific_latent_vaporization_heat':            'Lv',
        'water-liquid__mass-per-volume_density': 'rho_H2O' }

    #-----------------------------------------------------------------
    # Note: The "update()" function calls several functions with the
    #       MBAR keyword set to get units of "mbar" vs. "kPa".
    #-----------------------------------------------------------------
    # Note: We need to be careful with whether units are C or K,
    #       for all "thermal" quantities (e.g. thermal_capacity).
    #-----------------------------------------------------------------
    # Note: ARHYTHM had 3 "bulk exchange coefficients" that are all
    #       equal and therefore have the same units of [m s-1].
    #       Double-check that this is what is intended.   ##########
    #-----------------------------------------------------------------
    # Note: "atmosphere_column_water__liquid_equivalent_depth" has
    #       units of "cm", as in Dingman's book.  Make sure it gets
    #       used correctly in equations.
    #-----------------------------------------------------------------
    # Note: slope_angle and aspect_angle have units of RADIANS.
    #       aspect_angle is measured CW from north.
    #       RT files ending in "_mf-angle.rtg" and "fd-aspect.rtg"
    #       contain aspect values.  The former are in [0, 2 Pi]
    #       while the latter are in [-Pi, Pi] and both measure
    #       CCW from due east.  They are converted for use here.
    #-----------------------------------------------------------------    
    _var_units_map = {
        'snowpack__z_mean_of_mass-per-volume_density':     'kg m-3',
        'snowpack__depth':                   'm',
        'snowpack__liquid-equivalent_depth': 'm',
        'snowpack__melt_volume_flux':        'm s-1',
        #-------------------------------------------------------------
        # 'atmosphere__optical_path_length_ratio': '1',
        # 'atmosphere__von_karman_constant': '1',
        'atmosphere_aerosol_dust__reduction_of_transmittance':        '1',
        'atmosphere_air-column_water-vapor__liquid-equivalent_depth': 'cm',           # (see Notes above)
        'atmosphere_bottom_air__brutsaert_emissivity_canopy_factor':  '1',
        'atmosphere_bottom_air__brutsaert_emissivity_cloud_factor':   '1',
        'atmosphere_bottom_air__bulk_latent_heat_aerodynamic_conductance': 'm s-1',   # (see Notes above)
        'atmosphere_bottom_air__bulk_sensible_heat_aerodynamic_conductance': 'm s-1', # (see Notes above)  
        'atmosphere_bottom_air__emissivity': '1',                            
        'atmosphere_bottom_air__mass-per-volume_density': 'kg m-3',
        'atmosphere_bottom_air__mass-specific_isobaric_heat_capacity': 'J kg-1 K-1',   # (see Notes above)
        'atmosphere_bottom_air__neutral_bulk_heat_aerodynamic_conductance': 'm s-1',   # (see Notes above)
        'atmosphere_bottom_air__pressure':                               'mbar',
        'atmosphere_bottom_air__temperature':                            'deg_C',      # (see Notes above)
        'atmosphere_bottom_air_flow__bulk_richardson_number':            '1',
        'atmosphere_bottom_air_flow__log_law_roughness_length':          'm',
        'atmosphere_bottom_air_flow__reference-height_speed':            'm s-1',
        'atmosphere_bottom_air_flow__speed_reference_height':            'm',
        'atmosphere_bottom_air_land_net-latent-heat__energy_flux':       'W m-2',
        'atmosphere_bottom_air_land_net-sensible-heat__energy_flux':     'W m-2',
        'atmosphere_bottom_air_water-vapor__dew_point_temperature':      'deg_C',
        'atmosphere_bottom_air_water-vapor__partial_pressure':           'mbar', # (see Notes above)
        'atmosphere_bottom_air_water-vapor__relative_saturation':        '1',
        'atmosphere_bottom_air_water-vapor__saturated_partial_pressure': 'mbar',     # (see Notes above)
        'atmosphere_water__domain_time_integral_of_precipitation_leq-volume_flux': 'm3',
        'atmosphere_water__domain_time_max_of_precipitation_leq-volume_flux': 'm s-1',        
        'atmosphere_water__precipitation_leq-volume_flux': 'm s-1',
        'atmosphere_water__rainfall_volume_flux': 'm s-1',      # (see Notes above)  
        'atmosphere_water__snowfall_leq-volume_flux': 'm s-1',  # (see Notes above)
        'earth__standard_gravity_constant': 'm s-2',
        'land_surface__albedo': '1',
        'land_surface__aspect_angle': 'radians',                      # (see Notes above)
        'land_surface__emissivity': '1',
        'land_surface__latitude': 'degrees',
        'land_surface__longitude': 'degrees',
        'land_surface__slope_angle': 'radians',
        'land_surface__temperature': 'deg_C',
        # 'land_surface_air__temperature': 'deg_C',
        'land_surface_air_water-vapor__partial_pressure':           'mbar',
        'land_surface_air_water-vapor__saturated_partial_pressure': 'mbar',
        'land_surface_net-longwave-radiation__energy_flux': 'W m-2',
        'land_surface_net-shortwave-radiation__energy_flux': 'W m-2',
        'land_surface_net-total-energy__energy_flux': 'W m-2',
        'model__time_step': 's',
        'physics__stefan_boltzmann_constant': 'W m-2 K-4',
        'physics__von_karman_constant': '1',
        'water__mass-specific_latent_fusion_heat': 'J kg-1',
        'water__mass-specific_latent_vaporization_heat': 'J kg-1',
        'water-liquid__mass-per-volume_density': 'kg m-3' }

    #------------------------------------------------    
    # Return NumPy string arrays vs. Python lists ?
    #------------------------------------------------
    ## _input_var_names  = np.array( _input_var_names )
    ## _output_var_names = np.array( _output_var_names )

    #-------------------------------------------------------------------
    def get_component_name(self):
  
        return 'TopoFlow_Meteorology'

    #   get_component_name()         
    #-------------------------------------------------------------------
    def get_attribute(self, att_name):

        try:
            return self._att_map[ att_name.lower() ]
        except:
            print '###################################################'
            print ' ERROR: Could not find attribute: ' + att_name
            print '###################################################'
            print ' '

    #   get_attribute()
    #-------------------------------------------------------------------
    def get_input_var_names(self):

        #--------------------------------------------------------
        # Note: These are currently variables needed from other
        #       components vs. those read from files or GUI.
        #--------------------------------------------------------   
        return self._input_var_names
    
    #   get_input_var_names()
    #-------------------------------------------------------------------
    def get_output_var_names(self):
 
        return self._output_var_names
    
    #   get_output_var_names()
    #-------------------------------------------------------------------
    def get_var_name(self, long_var_name):
            
        return self._var_name_map[ long_var_name ]

    #   get_var_name()
    #-------------------------------------------------------------------
    def get_var_units(self, long_var_name):

        return self._var_units_map[ long_var_name ]
   
    #   get_var_units()
    #-------------------------------------------------------------------
##    def get_var_type(self, long_var_name):
##
##        #---------------------------------------
##        # So far, all vars have type "double",
##        # but use the one in BMI_base instead.
##        #---------------------------------------
##        return 'float64'
##    
##    #   get_var_type()
    #-------------------------------------------------------------------
    def set_constants(self):

        #---------------------------------
        # Define some physical constants
        #---------------------------------
        self.g        = np.float64(9.81)    # [m s-2, gravity]
        self.kappa    = np.float64(0.408)   # [1]  (von Karman)
        self.rho_H2O  = np.float64(1000)    # [kg m-3]
        self.rho_air  = np.float64(1.2614)  # [kg m-3]
        self.Cp_air   = np.float64(1005.7)  # [J kg-1 K-1]
        self.Lv       = np.float64(2500000) # [J kg-1] Latent heat of vaporiz.
        self.Lf       = np.float64(334000)  # [J kg-1 = W s kg-1], Latent heat of fusion
        self.sigma    = np.float64(5.67E-8) # [W m-2 K-4]  (Stefan-Boltzman constant)
        self.C_to_K   = np.float64(273.15)  # (add to convert deg C to K)

        self.twopi         = np.float64(2) * np.pi
        self.one_seventh   = np.float64(1) / 7
        self.hours_per_day = np.float64(24)
        self.secs_per_day  = np.float64(3600) * self.hours_per_day

        #---------------------------
        # See update_latent_heat()
        #-----------------------------------------------------------        
        # According to Dingman (2002, p. 273), constant should
        # be 0.622 instead of 0.662 (Zhang et al., 2000, p. 1002).
        # Is this constant actually the dimensionless ratio of
        # the molecular weight of water to that of dry air ?
        #-----------------------------------------------------------
        ## self.latent_heat_constant = np.float64(0.622)
        self.latent_heat_constant = np.float64(0.662)
        
        #----------------------------------------
        # Constants related to precip (9/24/09)
        #----------------------------------------
        self.mmph_to_mps = (np.float64(1) / np.float64(3600000))
        self.mps_to_mmph = np.float64(3600000)
        self.forever     = np.float64(999999999)  # [minutes]
        
        #------------------------------------------------
        # Only needed for method 1, where all rates and
        # durations are read as 1D arrays from GUI.
        # Method 1 may be removed in a future version.
        #------------------------------------------------
##        self.method1_rates     = None
##        self.method1_durations = None
##        self.method1_n_rates   = 0
        
    #   set_constants()             
    #-------------------------------------------------------------------
    def initialize(self, cfg_file=None, mode="nondriver",
                   SILENT=False):

        if not(SILENT):
            print ' '
            print 'Meteorology component: Initializing...'
            
        self.status   = 'initializing'  # (OpenMI 2.0 convention)
        self.mode     = mode
        self.cfg_file = cfg_file
                
        #-----------------------------------------------
        # Load component parameters from a config file
        #-----------------------------------------------
        self.set_constants()
        self.initialize_config_vars()     
        ## print '    Calling read_grid_info()...'
        self.read_grid_info()
        ## print '    Calling initialize_basin_vars()...'
        self.initialize_basin_vars()  # (5/14/10)
        #----------------------------------------------------
        # NB! This read_input_files() uses self.time_index.
        #     Also needs to be before "Disabled" test.
        #----------------------------------------------------
        ## print '    Calling initialize_time_vars()...'
        self.initialize_time_vars()

        #-------------------------------------------------
        # (5/19/12) This makes P "mutable", which allows
        # its updated values to be seen by any component
        # that has a reference to it.
        # Note that P will typically be read from a file.
        #-------------------------------------------------
        self.P      = self.initialize_var( self.P_type )
        self.P_rain = self.initialize_var( self.P_type )
        self.P_snow = self.initialize_var( self.P_type )
                    
        #------------------------------------------------------------
        # If "Enabled", will call initialize_computed_vars() below.
        #------------------------------------------------------------
        if (self.comp_status == 'Disabled'):
            if not(SILENT):
                print 'Meteorology component: Disabled in CFG file.'
            self.e_air    = self.initialize_scalar(0, dtype='float64')
            self.e_surf   = self.initialize_scalar(0, dtype='float64')
            self.em_air   = self.initialize_scalar(0, dtype='float64')
            self.Qn_SW    = self.initialize_scalar(0, dtype='float64')
            self.Qn_LW    = self.initialize_scalar(0, dtype='float64')
            self.Q_sum    = self.initialize_scalar(0, dtype='float64')
            self.Qc       = self.initialize_scalar(0, dtype='float64')
            self.Qa       = self.initialize_scalar(0, dtype='float64')
            self.DONE     = True
            self.status   = 'initialized'
            return

        #-----------------------------------------------
        # Read from files as needed to initialize vars 
        #-----------------------------------------------
        self.open_input_files()
        self.read_input_files()  # (initializes P)
        
        ## self.check_input_types()  # (not needed so far)
        
        #-----------------------
        # Initialize variables
        #-----------------------
        ## print '    Calling initialize_computed_vars()...'
        self.initialize_computed_vars() # (after read_input_files)
        
        if not(self.PRECIP_ONLY):
            self.open_output_files() 
        self.status = 'initialized'  # (OpenMI 2.0 convention)
        
    #   initialize()
    #-------------------------------------------------------------------
    ## def update(self, dt=-1.0, time_seconds=None):
    def update(self, dt=-1.0):

        #----------------------------------------------------------
        # Note: The read_input_files() method is first called by
        #       the initialize() method.  Then, the update()
        #       method is called one or more times, and it calls
        #       other update_*() methods to compute additional
        #       variables using input data that was last read.
        #       Based on this pattern, read_input_files() should
        #       be called at end of update() method as done here.
        #       If the input files don't contain any additional
        #       data, the last data read persists by default.
        #----------------------------------------------------------
        if (self.comp_status == 'Disabled'): return
        self.status = 'updating'  # (OpenMI 2.0 convention)

        #-------------------------------------------
        # Update computed values related to precip
        #-------------------------------------------
        self.update_P_integral()
        self.update_P_max()
        self.update_P_rain()
        self.update_P_snow()

        #-------------------------
        # Update computed values
        #-------------------------
        if not(self.PRECIP_ONLY):
            self.update_bulk_richardson_number()
            self.update_bulk_aero_conductance()
            self.update_sensible_heat_flux()
            self.update_saturation_vapor_pressure(MBAR=True)
            self.update_saturation_vapor_pressure(MBAR=True, SURFACE=True)  ########
            self.update_vapor_pressure(MBAR=True)
            self.update_dew_point() ###
            self.update_precipitable_water_content()  ###
            self.update_vapor_pressure(MBAR=True, SURFACE=True)   ########
            self.update_latent_heat_flux()      # (uses e_air and e_surf)
            self.update_conduction_heat_flux()
            self.update_advection_heat_flux()
            self.update_julian_day()
            self.update_net_shortwave_radiation()
            self.update_em_air()
            self.update_net_longwave_radiation()
            self.update_net_energy_flux()  # (at the end)
            
        #----------------------------------------
        # Read next met vars from input files ?
        #-------------------------------------------
        # Note that read_input_files() is called
        # by initialize() and these values must be
        # used for "update" calls before reading
        # new ones.
        #-------------------------------------------
        if (self.time_index > 0):
            self.read_input_files()
      
        #----------------------------------------------
        # Write user-specified data to output files ?
        #----------------------------------------------
        # Components use own self.time_sec by default.
        #-----------------------------------------------
        if not(self.PRECIP_ONLY):
            self.write_output_files()
            ## self.write_output_files( time_seconds )

        #-----------------------------
        # Update internal clock
        # after write_output_files()
        #-----------------------------
        self.update_time( dt )
        self.status = 'updated'  # (OpenMI)
        
    #   update()
    #-------------------------------------------------------------------
    def finalize(self):

        self.status = 'finalizing'  # (OpenMI)
        if (self.comp_status == 'Enabled'):
            self.close_input_files()   ##  TopoFlow input "data streams"
            if not(self.PRECIP_ONLY):
                self.close_output_files()
        self.status = 'finalized'  # (OpenMI)

        self.print_final_report(comp_name='Meteorology component')
     
    #   finalize()
    #-------------------------------------------------------------------
    def set_computed_input_vars(self):

        #-----------------------------------------------
        # Convert precip rate units from mm/h to m/s ?
        #-----------------------------------------------
        # NB! read_input_files() does this for files.
        #-----------------------------------------------
        if (self.P_type == 'Scalar'):
##            print '######## self.P_type  =', self.P_type
##            print '######## type(self.P) =', type(self.P)
##            print '######## self.P       =', self.P
##            print '######## Converting scalar P from MMPH to MPS.'
            #-----------------------------------------------------
            # (2/7/13) Must use "*=" here to preserve reference.
            #-----------------------------------------------------
            self.P *= self.mmph_to_mps
            ## self.P = self.P * self.mmph_to_mps
            print 'Scalar rainrate set to:', self.P, ' [mmph]'
            
        #---------------------------------
        # Process the PRECIP_ONLY toggle
        #---------------------------------
        if not(hasattr(self, 'PRECIP_ONLY')):
            self.PRECIP_ONLY = False
        elif (self.PRECIP_ONLY.lower() == 'yes'):
            self.PRECIP_ONLY = True
        else:
            self.PRECIP_ONLY = False

        #---------------------------------------
        # Print info message about PRECIP_ONLY
        #---------------------------------------
        if (self.PRECIP_ONLY):
            print '-----------------------------------------'
            print ' NOTE: Since PRECIP_ONLY = True, output'
            print '       variables will not be computed'
            print '       or saved to files.'
            print '-----------------------------------------'
            print' '
            
        #----------------------------------------------------
        # Toggle to use SATTERLUND or BRUTSAERT methods
        # for computing e_air and em_air. (Not in GUI yet.)
        #----------------------------------------------------
        if not(hasattr(self, 'SATTERLUND')):
            self.SATTERLUND = False

        #---------------------------------------------
        # Convert GMT_offset from string to int
        # because GUI can't use ints in droplist yet
        #---------------------------------------------
        self.GMT_offset = np.int16( self.GMT_offset )

        #------------------------------------------------
        # Convert start_month from string to integer
        # January should be 1.  See solar.Julian_Day().
        #------------------------------------------------
        month_list = ['January', 'February', 'March', 'April',
                      'May', 'June', 'July', 'August', 'September',
                      'October', 'November', 'December']
        self.start_month = month_list.index( self.start_month ) + 1
                 
        #-------------------------------
        # Initialize some more toggles
        #-------------------------------                   
        if not(hasattr(self, 'SAVE_QSW_GRIDS')):
            self.SAVE_QSW_GRIDS = False
        if not(hasattr(self, 'SAVE_QLW_GRIDS')):
            self.SAVE_QLW_GRIDS = False
        #-------------------------------------------
        if not(hasattr(self, 'SAVE_QSW_PIXELS')):
            self.SAVE_QSW_PIXELS = False
        if not(hasattr(self, 'SAVE_QLW_PIXELS')):
            self.SAVE_QLW_PIXELS = False
   
        #---------------------------------------------------------
        # Make sure that all "save_dts" are larger or equal to
        # the specified process dt.  There is no point in saving
        # results more often than they change.
        # Issue a message to this effect if any are smaller ??
        #---------------------------------------------------------
        self.save_grid_dt   = np.maximum(self.save_grid_dt,    self.dt)
        self.save_pixels_dt = np.maximum(self.save_pixels_dt,  self.dt)
        
    #   set_computed_input_vars()
    #-------------------------------------------------------------------
    def initialize_computed_vars(self):

        #------------------------------------------------------
        # Note: Some of these require "self.rti", which is
        #       only stored by read_grid_info() after the
        #       set_computed_input_vars() function is called.
        #       So these parts can't go there.
        #------------------------------------------------------
        
        #---------------------------------------
        # Add self.in_directory to:
        #   slope_grid_file & aspect_grid_file
        #---------------------------------------
        self.slope_grid_file  = (self.in_directory + self.slope_grid_file) 
        self.aspect_grid_file = (self.in_directory + self.aspect_grid_file)
        
        #-------------------------------------------------
        # Read slope grid & convert to slope angle, beta
        # NB!  RT slope grids have NaNs on edges.
        #-------------------------------------------------
        slopes = rtg_files.read_grid( self.slope_grid_file, self.rti,
                                      RTG_type='FLOAT' )
        beta   = np.arctan( slopes )
        beta   = (self.twopi + beta) % self.twopi
        #---------------------------------------------
        w_nan = np.where( np.logical_not(np.isfinite(beta)) )
        n_nan = np.size(w_nan[0])
        if (n_nan != 0):    
            beta[ w_nan ] = np.float64(0)
        #------------------------------------------------------------------
        w_bad = np.where( np.logical_or( (beta < 0), (beta > np.pi / 2) ) )
        n_bad = np.size(w_bad[0])
        if (n_bad != 0):    
            msg = array(['ERROR:  Some slope angles are out of range.', ' '])
            for line in msg:
                print line
            ## result = GUI_Message(msg, INFO=True, TITLE='ERROR MESSAGE')
            return
        self.beta = beta  ######
        
        #------------------------------------------------------
        # Read aspect grid.  Alpha must be CW from north.
        # NB!  RT aspect grids have NaNs on edges.
        #---------------------------------------------------------
        # RT files ending in "_mf-angle.rtg" and "fd-aspect.rtg"
        # contain aspect values.  The former are in [0, 2 Pi]
        # while the latter are in [-Pi, Pi] and both measure
        # CCW from due east.
        #---------------------------------------------------------
        aspects = rtg_files.read_grid( self.aspect_grid_file, self.rti,
                                       RTG_type='FLOAT' )
        alpha   = (np.pi / 2) - aspects
        alpha   = (self.twopi + alpha) % self.twopi
        #-----------------------------------------------
        w_nan = np.where( np.logical_not( np.isfinite(alpha) ) )
        n_nan = np.size( w_nan[0] )
        if (n_nan != 0):    
            alpha[ w_nan ] = np.float64(0)
        self.alpha = alpha  ######

        #---------------------------        
        # Create lon and lat grids
        #---------------------------
        if (self.rti.pixel_geom == 0):
            self.lon_deg = solar.Longitude_Grid( self.rti )
            self.lat_deg = solar.Latitude_Grid( self.rti )

##            print 'Lon grid ='
##            print self.lon_deg
##            print 'Lat grid ='
##            print self.lat_deg
            
            #-----------------------------
            # Write grids to RTG files ?
            #-----------------------------
##            lon_file = (self.out_directory + self.site_prefix + '_lons.bin')
##            rtg_files.write_grid( self.lon_deg, lon_file, self.rti )
##            lat_file = (self.out_directory + self.site_prefix + '_lats.bin')
##            rtg_files.write_grid( self.lat_deg, lat_file, self.rti )
        else:    
            print 'SORRY: Cannot yet create lon and lat grids for'
            print '       this DEM because it uses UTM coordinates.'
            print '       Will use lat/lon for Denver, Colorado.'
            print ' '
            #--------------------------------------------
            # For now, use scalar values for Denver, CO
            #--------------------------------------------
            self.lon_deg = np.float64( -104.9841667 )
            self.lat_deg = np.float64( 39.7391667 )
            ## return

        #-------------------------------------------------
        # Initialize max precip rate with the first rate
        #------------------------------------------------
        # Note: Need this here because rate may be
        #       zero at the end of update_precip_rate()
        #------------------------------------------------
        # vol_P is used for mass balance check.
        #------------------------------------------------
        P_max = self.P.max()       # (after read_input_files)
        ## self.P_max = self.P.max()
        self.P_max = self.initialize_scalar( P_max, dtype='float64')
        self.vol_P = self.initialize_scalar( 0, dtype='float64')


        #----------------------------------------------------------
        # For using new framework which embeds references from
        # meteorology to snow, etc., these need to be defined
        # in the initialize step.  However, they will most likely
        # change from scalar to grid during update, so we need to
        # check that the reference isn't broken when the dtype
        # changes. (5/17/12)
        #---------------------------------------------------------- 
        # These depend on grids alpha and beta, so will be grids.
        #----------------------------------------------------------                
        self.Qn_SW  = np.zeros([self.ny, self.nx], dtype='float64')
        self.Qn_LW  = np.zeros([self.ny, self.nx], dtype='float64')
        self.Qn_tot = np.zeros([self.ny, self.nx], dtype='float64')
        self.Q_sum  = np.zeros([self.ny, self.nx], dtype='float64')
        #----------------------------------------------------------  
#         self.Qn_SW  = self.initialize_scalar( 0, dtype='float64')
#         self.Qn_LW  = self.initialize_scalar( 0, dtype='float64')
#         self.Qn_tot = self.initialize_scalar( 0, dtype='float64')
#         self.Q_sum  = self.initialize_scalar( 0, dtype='float64')
        #---------------------------------------------------------- 
        # These may be scalars or grids.
        #---------------------------------         
        self.Qe     = self.initialize_scalar( 0, dtype='float64')
        self.e_air  = self.initialize_scalar( 0, dtype='float64')
        self.e_surf = self.initialize_scalar( 0, dtype='float64')
        self.em_air = self.initialize_scalar( 0, dtype='float64')
        self.Qc     = self.initialize_scalar( 0, dtype='float64')
        self.Qa     = self.initialize_scalar( 0, dtype='float64')
         
        #------------------------------------
        # Initialize the decimal Julian day
        #------------------------------------
        self.julian_day = solar.Julian_Day( self.start_month,
                                            self.start_day,
                                            self.start_hour )
        ## print '    julian_day =', self.julian_day
        
    #   initialize_computed_vars()
    #-------------------------------------------------------------------
    def update_P_integral(self):

        #---------------------------------------------------
        # Notes: This can be used for mass balance checks,
        #        such as now done by update_mass_totals()
        #        in topoflow.py.  The "dt" here should be
        #        TopoFlow's "main dt" vs. the process dt.
        
        #        dV[i] = P[i] * da[i] * dt, dV = sum(dV[i])
        #---------------------------------------------------
        if (self.DEBUG):
            print 'Calling update_P_integral()...'
            
        #------------------------------------------------
        # Update mass total for P, sum over all pixels
        #------------------------------------------------   
        volume = np.double(self.P * self.da * self.dt)  # [m^3]
        if (np.size(volume) == 1):
            self.vol_P += (volume * self.rti.n_pixels)
        else:
            self.vol_P += np.sum(volume)
        
    #   update_P_integral()
    #-------------------------------------------------------------------
    def update_P_max(self):

        if (self.DEBUG):
            print 'Calling update_P_max()...'
            
        #-----------------------------------------
        # Save the maximum precip. rate in [m/s]
        #-------------------------------------------
        # Must use "fill()" to preserve reference.
        #-------------------------------------------
        self.P_max.fill( np.maximum(self.P_max, self.P.max()) )
        ## self.P_max = np.maximum(self.P_max, self.P.max())

        #---------------
        # For testing
        #--------------
        ## print '##### P =', self.P * 1000 * 3600  # (mmph)
        ## print '##### P_max =', self.P_max * 1000 * 3600  # (mmph)

    #   update_P_max()  
    #-------------------------------------------------------------------
    def update_P_rain(self):

        #-----------------------------------------------------------
        # Note:  This routine is written so that it doesn't matter
        #        whether P and T_air are grids or scalars.
        #        For scalars: 1.5 * True = 1.5, 1.5 * False = 0.
        #        Here are the possible combinations for checking.
        #-----------------------------------------------------------
        # P       T_air     P_rain
        #----------------------------
        # scalar  scalar    scalar    
        # scalar  grid      grid
        # grid    scalar    grid
        # grid    grid      grid
        #----------------------------
        if (self.DEBUG):
            print 'Calling update_P_rain()...'

        #-------------------------------------------------
        # P_rain is the precip that falls as liquid that
        # can contribute to runoff production.
        #-------------------------------------------------
        # P_rain is used by channel_base.update_R.
        #-------------------------------------------------
        P_rain = self.P * (self.T_air > 0)
        self.update_var( 'P_rain', P_rain )  ## (2/14/17)
#         if (np.ndim( self.P_rain ) == 0):
#             self.P_rain.fill( P_rain )   #### (mutable scalar)
#         else:
#             self.P_rain[:] = P_rain
  
        if (self.DEBUG):
            if (self.P_rain.max() > 0):
                print '   >> Rain is falling...'

        #--------------
        # For testing
        #--------------
        ## print 'shape(P)      =', shape(self.P)
        ## print 'shape(T_air)  =', shape(self.T_air)
        ## print 'shape(P_rain) =', shape(self.P_rain)
        ## print 'T_air         =', self.T_air

        #########################################
        #### Old note, to remember for later.
        #--------------------------------------------------
        # (2/7/13) We must use "*=" to preserve reference
        # if P is a "mutable scalar".
        #--------------------------------------------------
                             
    #   update_P_rain()
    #-------------------------------------------------------------------
    def update_P_snow(self):

        #----------------------------------------------------
        # Notes:  Rain and snow may fall simultaneously at
        #         different grid cells in the model domain.
        #----------------------------------------------------
        if (self.DEBUG):
            print 'Calling update_P_snow()...'
            
        #-------------------------------------------------
        # P_snow is the precip that falls as snow or ice
        # that contributes to the snow depth.  This snow
        # may melt to contribute to runoff later on.
        #-------------------------------------------------
        # P_snow is used by snow_base.update_depth.
        #-------------------------------------------------
        P_snow = self.P * (self.T_air <= 0)
        self.update_var( 'P_snow', P_snow )   ## (2/14/17)
#         if (np.ndim( self.P_snow ) == 0):
#             self.P_snow.fill( P_snow )   #### (mutable scalar)
#         else:
#             self.P_snow[:] = P_snow

        if (self.DEBUG):
            if (self.P_snow.max() > 0):
                print '   >> Snow is falling...'
                     
    #   update_P_snow()
    #-------------------------------------------------------------------
    def update_bulk_richardson_number(self):

        if (self.DEBUG):
            print 'Calling update_bulk_richardson_number()...'

        #---------------------------------------------------------------
        # (9/6/14)  Found a typo in the Zhang et al. (2000) paper,
        # in the definition of Ri.  Also see Price and Dunne (1976).
        # We should have (Ri > 0) and (T_surf > T_air) when STABLE.
        # This also removes problems/singularities in the corrections
        # for the stable and unstable cases in the next function.      
        #---------------------------------------------------------------
        # Notes: Other definitions are possible, such as the one given
        #        by Dingman (2002, p. 599).  However, this one is the
        #        one given by Zhang et al. (2000) and is meant for use
        #        with the stability criterion also given there.
        #---------------------------------------------------------------
        #### top     = self.g * self.z * (self.T_air - self.T_surf)  # BUG.
        top     = self.g * self.z * (self.T_surf - self.T_air)
        bot     = (self.uz)**2.0 * (self.T_air + np.float64(273.15))
        self.Ri = (top / bot)

    #   update_bulk_richardson_number()
    #-------------------------------------------------------------------        
    def update_bulk_aero_conductance(self):

        if (self.DEBUG):
            print 'Calling update_bulk_aero_conductance()...'

        #----------------------------------------------------------------
        # Notes: Dn       = bulk exchange coeff for the conditions of
        #                   neutral atmospheric stability [m/s]
        #        Dh       = bulk exchange coeff for heat  [m/s]
        #        De       = bulk exchange coeff for vapor [m/s]
        #        h_snow   = snow depth [m]
        #        z0_air   = surface roughness length scale [m]
        #                   (includes vegetation not covered by snow)
        #        z        = height that has wind speed uz [m]
        #        uz       = wind speed at height z [m/s]
        #        kappa    = 0.408 = von Karman's constant [unitless]
        #        RI       = Richardson's number (see function)
        #----------------------------------------------------------------
        h_snow = self.h_snow  # (ref from new framework)
        
        #---------------------------------------------------
        # Compute bulk exchange coeffs (neutral stability)
        # using the logarithm "law of the wall".
        #-----------------------------------------------------
        # Note that "arg" = the drag coefficient (unitless).
        #-----------------------------------------------------        
        arg = self.kappa / np.log((self.z - h_snow) / self.z0_air)
        Dn  = self.uz * (arg)**2.0
        
        #-----------------------------------------------
        # NB! Dn could be a scalar or a grid, so this
        #     must be written to handle both cases.
        #     Note that WHERE can be used on a scalar:
        
        #     IDL> a = 1
        #     IDL> print, size(a)
        #     IDL> w = where(a ge 1, nw)
        #     IDL> print, nw
        #     IDL> a[w] = 2
        #     IDL> print, a
        #     IDL> print, size(a)
        #-----------------------------------------------

        ###########################################################
        #  NB!  If T_air and T_surf are both scalars, then next
        #       few lines won't work because we can't index the
        #       resulting empty "w" (even if T_air == T_surf).
        ###########################################################
##        w  = np.where(self.T_air != self.T_surf)
##        nw = np.size(w[0])
##        ## nw = np.size(w,0)  # (doesn't work if 2 equal scalars)
        #----------------------------------------------------------
        T_AIR_SCALAR  = (np.ndim( self.T_air )  == 0)
        T_SURF_SCALAR = (np.ndim( self.T_surf ) == 0)
        if (T_AIR_SCALAR and T_SURF_SCALAR):
            if (self.T_air == self.T_surf):  nw=1
            else: nw=0      
        else:
            w  = np.where(self.T_air != self.T_surf)
            nw = np.size(w[0])
        
        if (nw == 0):
            #--------------------------------------------
            # All pixels are neutral. Set Dh = De = Dn.
            #--------------------------------------------
            self.Dn = Dn
            self.Dh = Dn
            self.De = Dn
            return
        
        #-------------------------------------
        # One or more pixels are not neutral
        # so make a correction using RI
        #---------------------------------------------
        # NB!  RI could be a grid when Dn is a
        # scalar, and this will change Dn to a grid.
        #---------------------------------------------
        # Ri = Richardson_Number(z, uz, T_air, T_surf)
        #--------------------------------------------
        # Before 12/21/07.  Has bug if RI is a grid
        #--------------------------------------------
        # w_stable = where(*T_air gt *T_surf, n_stable)
        # if (n_stable ne 0) then begin
        #     Dn[w_stable] = Dn[w_stable]/(1d + (10d * RI))
        # endif
        # w_unstable = where(*T_air lt *T_surf, n_unstable)
        # if (n_unstable ne 0) then begin
        #----------------------------------------------
        # Multiplication and substraction vs. opposites
        # for the stable case.  Zhang et al. (2000)
        # Hopefully not just a typo.
        #----------------------------------------------
        #    Dn[w_unstable] = Dn[w_unstable]*(1d - (10d * self.Ri))
        # endif
        
        #-----------------
        # After 12/21/07
        #------------------------------------------------------------
        # If T_air, T_surf or uz is a grid, then Ri will be a grid.
        # This version makes only one call to WHERE, so its faster.
        #------------------------------------------------------------
        # Multiplication and substraction vs. opposites for the
        # stable case (Zhang et al., 2000); hopefully not a typo.
        # It plots as a smooth curve through Ri=0.
        #------------------------------------------------------------
        # (9/7/14)  Modified so that Dn is saved, but Dh = De.
        #------------------------------------------------------------        
        Dh = Dn.copy()   ### (9/7/14.  Save Dn also.)
        nD = np.size( Dh )
        nR = np.size( self.Ri )
        if (nR > 1):    
            #--------------------------
            # Case where RI is a grid
            #--------------------------
            ws = np.where( self.Ri > 0 )
            ns = np.size( ws[0] )
            wu = np.where( np.invert(self.Ri > 0) )
            nu = np.size( wu[0] )
            if (nD == 1):    
                #******************************************
                # Convert Dn to a grid here or somewhere
                # Should stop with an error message
                #******************************************
                dum = np.int16(0)
            if (ns != 0):
                #----------------------------------------------------------
                # If (Ri > 0), or (T_surf > T_air), then STABLE. (9/6/14)
                #----------------------------------------------------------   
                Dh[ws] = Dh[ws] / (np.float64(1) + (np.float64(10) * self.Ri[ws]))
            if (nu != 0):    
                Dh[wu] = Dh[wu] * (np.float64(1) - (np.float64(10) * self.Ri[wu]))
        else:    
            #----------------------------
            # Case where Ri is a scalar
            #--------------------------------
            # Works if Dh is grid or scalar
            #--------------------------------
            if (self.Ri > 0):    
                Dh = Dh / (np.float64(1) + (np.float64(10) * self.Ri))
            else:    
                Dh = Dh * (np.float64(1) - (np.float64(10) * self.Ri))

        #----------------------------------------------------
        # NB! We currently assume that these are all equal.
        #----------------------------------------------------
        self.Dn = Dn
        self.Dh = Dh
        self.De = Dh   ## (assumed equal)
        
    #   update_bulk_aero_conductance()
    #-------------------------------------------------------------------
    def update_sensible_heat_flux(self):

        #--------------------------------------------------------
        # Notes: All the Q's have units of W/m^2 = J/(m^2 s).
        #        Dh is returned by Bulk_Exchange_Coeff function
        #        and is not a pointer.
        #--------------------------------------------------------
        if (self.DEBUG):
            print 'Callilng update_sensible_heat_flux()...'
            
        #---------------------
        # Physical constants
        #---------------------
        # rho_air = 1.225d   ;[kg m-3, at sea-level]
        # Cp_air  = 1005.7   ;[J kg-1 K-1]
        
        #-----------------------------
        # Compute sensible heat flux
        #-----------------------------
        delta_T = (self.T_air - self.T_surf)
        self.Qh = (self.rho_air * self.Cp_air) * self.Dh * delta_T

    #   update_sensible_heat_flux()
    #-------------------------------------------------------------------
    def update_saturation_vapor_pressure(self, MBAR=False,
                                         SURFACE=False):

        if (self.DEBUG):
            print 'Calling update_saturation_vapor_pressure()...'

        #----------------------------------------------------------------
        #Notes:  Saturation vapor pressure is a function of temperature.
        #        T is temperature in Celsius.  By default, the method
        #        of Brutsaert (1975) is used.  However, the SATTERLUND
        #        keyword is set then the method of Satterlund (1979) is
        #        used.  When plotted, they look almost identical.  See
        #        the Compare_em_air_Method routine in Qnet_file.pro.
        #        Dingman (2002) uses the Brutsaert method.
        #        Liston (1995, EnBal) uses the Satterlund method.

        #        By default, the result is returned with units of kPa.
        #        Set the MBAR keyword for units of millibars.
        #        100 kPa = 1 bar = 1000 mbars
        #                => 1 kPa = 10 mbars
        #----------------------------------------------------------------
        #NB!     Here, 237.3 is correct, and not a misprint of 273.2.
        #        See footnote on p. 586 in Dingman (Appendix D).
        #----------------------------------------------------------------
        if (SURFACE):
##            if (self.T_surf_type in ['Scalar', 'Grid']):
##                return
            T = self.T_surf
        else:
##            if (self.T_air_type in ['Scalar', 'Grid']):
##                return
            T = self.T_air
        
        if not(self.SATTERLUND):    
            #------------------------------
            # Use Brutsaert (1975) method
            #------------------------------
            term1 = (np.float64(17.3) * T) / (T + np.float64(237.3))
            e_sat = np.float64(0.611) * np.exp(term1)        # [kPa]
        else:    
            #-------------------------------
            # Use Satterlund (1979) method     ############ DOUBLE CHECK THIS (7/26/13)
            #-------------------------------
            term1 = np.float64(2353) / (T + np.float64(273.15))
            e_sat = np.float64(10) ** (np.float64(11.4) - term1)   # [Pa]
            e_sat = (e_sat / np.float64(1000))   # [kPa]

        #-----------------------------------
        # Convert units from kPa to mbars?
        #-----------------------------------
        if (MBAR):    
            e_sat = (e_sat * np.float64(10))   # [mbar]

        if (SURFACE):
            self.e_sat_surf = e_sat
        else:
            self.e_sat_air  = e_sat

    #   update_saturation_vapor_pressure()
    #------------------------------------------------------------------- 
    def update_vapor_pressure(self, MBAR=False,
                              SURFACE=False):

        if (self.DEBUG):
            print 'Calling update_vapor_pressure()...'

        #---------------------------------------------------
        # Notes: T is temperature in Celsius
        #        RH = relative humidity, in [0,1]
        #             by definition, it equals (e / e_sat)
        #        e has units of kPa.
        #---------------------------------------------------
        if (SURFACE):
##            if (self.T_surf_type in ['Scalar', 'Grid']) and \
##               (self.RH_type in ['Scalar', 'Grid']):
##                return
            e_sat = self.e_sat_surf
        else:
##            if (self.T_air_type in ['Scalar', 'Grid']) and \
##               (self.RH_type in ['Scalar', 'Grid']):
##                return
            e_sat = self.e_sat_air
            
        e = (self.RH * e_sat)
        
        #-----------------------------------
        # Convert units from kPa to mbars?
        #-----------------------------------
        if (MBAR):    
            e = (e * np.float64(10))   # [mbar]

        if (SURFACE):
            self.e_surf = e
        else:
            self.e_air  = e
 
    #   update_vapor_pressure()
    #-------------------------------------------------------------------
    def update_dew_point(self):

        if (self.DEBUG):
            print 'Calling update_dew_point()...'

        #-----------------------------------------------------------
        # Notes:  The dew point is a temperature in degrees C and
        #         is a function of the vapor pressure, e_air.
        #         Vapor pressure is a function of air temperature,
        #         T_air, and relative humidity, RH.
        #         The formula used here needs e_air in kPa units.
        #         See Dingman (2002, Appendix D, p. 587).
        #-----------------------------------------------------------
        e_air_kPa = self.e_air / np.float64(10)   # [kPa]
        log_vp    = np.log( e_air_kPa )

        top = log_vp + np.float64(0.4926)
        bot = np.float64(0.0708) - (np.float64(0.00421) * log_vp)

        self.T_dew = (top / bot)    # [degrees C]
    
    #   update_dew_point()
    #-------------------------------------------------------------------
    def update_precipitable_water_content(self):

        if (self.DEBUG):
            print 'Calling update_precipitable_water_content()...'

        #------------------------------------------------------------
        # Notes:  W_p is precipitable water content in centimeters,
        #         which depends on air temp and relative humidity.
        #------------------------------------------------------------
        arg      = np.float64( 0.0614 * self.T_dew )
        self.W_p = np.float64(1.12) * np.exp( arg )  # [cm]

    #   update_precipitable_water_content()
    #-------------------------------------------------------------------
    def update_latent_heat_flux(self):

        if (self.DEBUG):
            print 'Calling update_latent_heat_flux()...'

        #--------------------------------------------------------
        # Notes:  Pressure units cancel out because e_air and
        #         e_surf (in numer) have same units (mbar) as
        #         p0 (in denom).
        #--------------------------------------------------------        
        # According to Dingman (2002, p. 273), constant should
        # be 0.622 instead of 0.662 (Zhang et al., 2000).
        #--------------------------------------------------------
        const   = self.latent_heat_constant
        factor  = (self.rho_air * self.Lv * self.De)
        delta_e = (self.e_air - self.e_surf)
        self.Qe = factor * delta_e * (const / self.p0)

    #   update_latent_heat_flux()
    #-------------------------------------------------------------------
    def update_conduction_heat_flux(self):

        if (self.DEBUG):
            print 'Calling update_conduction_heat_flux()...'

        #-----------------------------------------------------------------
        # Notes: The conduction heat flux from snow to soil for computing
        #        snowmelt energy, Qm, is close to zero.

        #        However, the conduction heat flux from surface and sub-
        #        surface for computing Qet is given by Fourier's Law,
        #        namely Qc = Ks(Tx - Ts)/x.

        #        All the Q's have units of W/m^2 = J/(m^2 s).
        #-----------------------------------------------------------------
        pass  # (initialized at start)

    #   update_conduction_heat_flux()
    #-------------------------------------------------------------------
    def update_advection_heat_flux(self):

        if (self.DEBUG):
            print 'Calling update_advection_heat_flux()...'

        #------------------------------------------------------
        # Notes: All the Q's have units of W/m^2 = J/(m^2 s).
        #------------------------------------------------------
        pass  # (initialized at start)
        
    #   update_advection_heat_flux()
    #-------------------------------------------------------------------
    def update_julian_day(self):

        if (self.DEBUG):
            print 'Calling update_julian_day()...'

        #----------------------------------
        # Update the *decimal* Julian day
        #----------------------------------
        self.julian_day += (self.dt / self.secs_per_day) # [days]
  
        #------------------------------------------
        # Compute the offset from True Solar Noon
        # clock_hour is in 24-hour military time
        # but it can have a decimal part.
        #------------------------------------------
        dec_part   = self.julian_day - np.int16(self.julian_day)
        clock_hour = dec_part * self.hours_per_day
        ## print '    Computing solar_noon...'
        solar_noon = solar.True_Solar_Noon( self.julian_day,
                                            self.lon_deg,
                                            self.GMT_offset )
        ## print '    Computing TSN_offset...'
        self.TSN_offset = (clock_hour - solar_noon)    # [hours]
    
    #   update_julian_day()
    #-------------------------------------------------------------------
    def update_net_shortwave_radiation(self):

        #---------------------------------------------------------
        # Notes:  If time is before local sunrise or after local
        #         sunset then Qn_SW should be zero.
        #---------------------------------------------------------
        if (self.DEBUG):
            print 'Calling update_net_shortwave_radiation()...'

        #---------------------------------------
        # Compute Qn_SW for this time  [W m-2]
        #---------------------------------------
        Qn_SW = solar.Clear_Sky_Radiation( self.lat_deg,
                                           self.julian_day,
                                           self.W_p,
                                           self.TSN_offset,
                                           self.alpha,
                                           self.beta,
                                           self.albedo,
                                           self.dust_atten )

        self.update_var( 'Qn_SW', Qn_SW )  ## (2/14/17)
#         if (np.ndim( self.Qn_SW ) == 0):
#             self.Qn_SW.fill( Qn_SW )   #### (mutable scalar)
#         else:
#             self.Qn_SW[:] = Qn_SW  # [W m-2]
        
    #   update_net_shortwave_radiation()
    #-------------------------------------------------------------------
    def update_em_air(self):

        if (self.DEBUG):
            print 'Calling update_em_air()...'

        #---------------------------------------------------------
        # NB!  The Brutsaert and Satterlund formulas for air
        #      emissivity as a function of air temperature are in
        #      close agreement; see compare_em_air_methods().
        #      However, we must pay close attention to whether
        #      equations require units of kPa, Pa, or mbar.
        #
        #             100 kPa = 1 bar = 1000 mbars
        #                => 1 kPa = 10 mbars
        #---------------------------------------------------------
        # NB!  Temperatures are assumed to be given with units
        #      of degrees Celsius and are converted to Kelvin
        #      wherever necessary by adding C_to_K = 273.15.
        #
        #      RH = relative humidity [unitless]
        #---------------------------------------------------------
        # NB!  I'm not sure about how F is added at end because
        #      of how the equation is printed in Dingman (2002).
        #      But it reduces to other formulas as it should.
        #---------------------------------------------------------
        T_air_K = self.T_air + self.C_to_K
        
        if not(self.SATTERLUND):
            #-----------------------------------------------------
            # Brutsaert (1975) method for computing emissivity
            # of the air, em_air.  This formula uses e_air with
            # units of kPa. (From Dingman (2002, p. 196).)
            # See notes for update_vapor_pressure().
            #-----------------------------------------------------
            e_air_kPa = self.e_air / np.float64(10)  # [kPa]
            F       = self.canopy_factor
            C       = self.cloud_factor
            term1   = (1.0 - F) * 1.72 * (e_air_kPa / T_air_K) ** self.one_seventh
            term2   = (1.0 + (0.22 * C ** 2.0))
            self.em_air  = (term1 * term2) + F
        else:
            #--------------------------------------------------------
            # Satterlund (1979) method for computing the emissivity
            # of the air, em_air, that is intended to "correct
            # apparent deficiencies in this formulation at air
            # temperatures below 0 degrees C" (see G. Liston)
            # Liston cites Aase and Idso(1978), Satterlund (1979)
            #--------------------------------------------------------
            e_air_mbar = self.e_air
            eterm  = np.exp(-1 * (e_air_mbar)**(T_air_K / 2016) )
            self.em_air = 1.08 * (1.0 - eterm)
 
        #--------------------------------------------------------------  
        # Can't do this yet.  em_air is always initialized scalar now
        # but may change to grid on assignment. (9/23/14)
        #--------------------------------------------------------------
#         if (np.ndim( self.em_air ) == 0):
#             self.em_air.fill( em_air )   #### (mutable scalar)
#         else:
#             self.em_air[:] = em_air
           
    #   update_em_air()    
    #-------------------------------------------------------------------
    def update_net_longwave_radiation(self):

        #----------------------------------------------------------------
        # Notes: Net longwave radiation is computed using the
        #        Stefan-Boltzman law.  All four data types
        #        should be allowed (scalar, time series, grid or
        #        grid stack).
        #
        #        Qn_LW = (LW_in - LW_out)
        #        LW_in   = em_air  * sigma * (T_air  + 273.15)^4
        #        LW_out  = em_surf * sigma * (T_surf + 273.15)^4
        #
        #        Temperatures in [deg_C] must be converted to
        #        [K].  Recall that absolute zero occurs at
        #        0 [deg_K] or -273.15 [deg_C].
        #
        #----------------------------------------------------------------
        # First, e_air is computed as:
        #   e_air = RH * 0.611 * exp[(17.3 * T_air) / (T_air + 237.3)]
        # Then, em_air is computed as:
        #   em_air = (1 - F) * 1.72 * [e_air / (T_air + 273.15)]^(1/7) *
        #             (1 + 0.22 * C^2) + F
        #----------------------------------------------------------------
        if (self.DEBUG):
            print 'Calling update_net_longwave_radiation()...'

        #--------------------------------
        # Compute Qn_LW for this time
        #--------------------------------
        T_air_K  = self.T_air  + self.C_to_K
        T_surf_K = self.T_surf + self.C_to_K
        LW_in    = self.em_air  * self.sigma * (T_air_K)** 4.0
        LW_out   = self.em_surf * self.sigma * (T_surf_K)** 4.0
        LW_out   = LW_out + ((1.0 - self.em_surf) * LW_in)
               
        self.Qn_LW = (LW_in - LW_out)   # [W m-2]

        #--------------------------------------------------------------  
        # Can't do this yet.  Qn_LW is always initialized grid now
        # but will often be created above as a scalar. (9/23/14)
        #--------------------------------------------------------------
#         if (np.ndim( self.Qn_LW ) == 0):
#             self.Qn_LW.fill( Qn_LW )   #### (mutable scalar)
#         else:
#             self.Qn_LW[:] = Qn_LW  # [W m-2]
        
    #   update_net_longwave_radiation()
    #-------------------------------------------------------------------
    def update_net_total_radiation(self):

        #-----------------------------------------------
        # Notes: Added this on 9/11/14.  Not used yet.
        #------------------------------------------------------------
        #        Qn_SW = net shortwave radiation flux (solar)
        #        Qn_LW = net longwave radiation flux (air, surface)
        #------------------------------------------------------------       
        if (self.DEBUG):
            print 'Calling update_net_total_radiation()...'
            
        Qn_tot = self.Qn_SW + self.Qn_LW   # [W m-2]

        self.update_var( 'Qn_tot', Qn_tot )   ## (2/14/17)
#         if (np.ndim( self.Qn_tot ) == 0):
#             self.Qn_tot.fill( Qn_tot )   #### (mutable scalar)
#         else:
#             self.Qn_tot[:] = Qn_tot  # [W m-2]
                       
    #   update_net_total_radiation()
    #-------------------------------------------------------------------
    def update_net_energy_flux(self):

        if (self.DEBUG):
            print 'Calling update_net_energy_flux()...'

        #------------------------------------------------------
        # Notes: Q_sum is used by "snow_energy_balance.py".
        #------------------------------------------------------
        #        Qm    = energy used to melt snowpack (if > 0)
        #        Qn_SW = net shortwave radiation flux (solar)
        #        Qn_LW = net longwave radiation flux (air, surface)
        #        Qh    = sensible heat flux from turbulent convection
        #                between snow surface and air
        #        Qe    = latent heat flux from evaporation, sublimation,
        #                and condensation
        #        Qa    = energy advected by moving water (i.e. rainfall)
        #                (ARHYTHM assumes this to be negligible; Qa=0.)
        #        Qc    = energy flux via conduction from snow to soil
        #                (ARHYTHM assumes this to be negligible; Qc=0.)
        #        Ecc   = cold content of snowpack = amount of energy
        #                needed before snow can begin to melt [J m-2]

        #        All Q's here have units of [W m-2].
        #        Are they all treated as positive quantities ?

        #        rho_air  = density of air [kg m-3]
        #        rho_snow = density of snow [kg m-3]
        #        Cp_air   = specific heat of air [J kg-1 K-1]
        #        Cp_snow  = heat capacity of snow [J kg-1 K-1]
        #                 = ???????? = specific heat of snow
        #        Kh       = eddy diffusivity for heat [m2 s-1]
        #        Ke       = eddy diffusivity for water vapor [m2 s-1]
        #        Lv       = latent heat of vaporization [J kg-1]
        #        Lf       = latent heat of fusion [J kg-1]
        #        ------------------------------------------------------
        #        Dn       = bulk exchange coeff for the conditions of
        #                   neutral atmospheric stability [m/s]
        #        Dh       = bulk exchange coeff for heat
        #        De       = bulk exchange coeff for vapor
        #        ------------------------------------------------------
        #        T_air    = air temperature [deg_C]
        #        T_surf   = surface temperature [deg_C]
        #        T_snow   = average snow temperature [deg_C]
        #        RH       = relative humidity [unitless] (in [0,1])
        #        e_air    = air vapor pressure at height z [mbar]
        #        e_surf   = surface vapor pressure [mbar]
        #        ------------------------------------------------------
        #        h_snow   = snow depth [m]
        #        z        = height where wind speed is uz [m]
        #        uz       = wind speed at height z [m/s]
        #        p0       = atmospheric pressure [mbar]
        #        T0       = snow temperature when isothermal [deg_C]
        #                   (This is usually 0.)
        #        z0_air   = surface roughness length scale [m]
        #                   (includes vegetation not covered by snow)
        #                   (Values from page 1033: 0.0013, 0.02 [m])
        #        kappa    = von Karman's constant [unitless] = 0.41
        #        dt       = snowmelt timestep [seconds]
        #----------------------------------------------------------------
        Q_sum = self.Qn_SW + self.Qn_LW + self.Qh + \
                self.Qe + self.Qa + self.Qc    # [W m-2]

        self.update_var( 'Q_sum', Q_sum )   ## (2/14/17)
#         if (np.ndim( self.Q_sum) == 0):
#             self.Q_sum.fill( Q_sum )   #### (mutable scalar)
#         else:
#             self.Q_sum[:] = Q_sum  # [W m-2]
            
    #   update_net_energy_flux()   
    #-------------------------------------------------------------------  
    def open_input_files(self):

        if (self.DEBUG):
            print 'Calling open_input_files()...'
       
        self.P_file      = self.in_directory + self.P_file
        self.T_air_file  = self.in_directory + self.T_air_file
        self.T_surf_file = self.in_directory + self.T_surf_file
        self.RH_file     = self.in_directory + self.RH_file
        self.p0_file     = self.in_directory + self.p0_file
        self.uz_file     = self.in_directory + self.uz_file
        self.z_file      = self.in_directory + self.z_file
        self.z0_air_file = self.in_directory + self.z0_air_file

        self.albedo_file        = self.in_directory + self.albedo_file
        self.em_surf_file       = self.in_directory + self.em_surf_file
        self.dust_atten_file    = self.in_directory + self.dust_atten_file
        self.cloud_factor_file  = self.in_directory + self.cloud_factor_file
        self.canopy_factor_file = self.in_directory + self.canopy_factor_file

        self.P_unit      = model_input.open_file(self.P_type,      self.P_file)
        self.T_air_unit  = model_input.open_file(self.T_air_type,  self.T_air_file)
        self.T_surf_unit = model_input.open_file(self.T_surf_type, self.T_surf_file)
        self.RH_unit     = model_input.open_file(self.RH_type,     self.RH_file)
        self.p0_unit     = model_input.open_file(self.p0_type,     self.p0_file)
        self.uz_unit     = model_input.open_file(self.uz_type,     self.uz_file)
        self.z_unit      = model_input.open_file(self.z_type,      self.z_file)
        self.z0_air_unit = model_input.open_file(self.z0_air_type, self.z0_air_file)
               
        #-----------------------------------------------
        # These are needed to compute Qn_SW and Qn_LW.
        #-----------------------------------------------
        self.albedo_unit        = model_input.open_file(self.albedo_type,
                                                        self.albedo_file)
        self.em_surf_unit       = model_input.open_file(self.em_surf_type,
                                                        self.em_surf_file)
        self.dust_atten_unit    = model_input.open_file(self.dust_atten_type,
                                                        self.dust_atten_file)
        self.cloud_factor_unit  = model_input.open_file(self.cloud_factor_type,
                                                        self.cloud_factor_file)
        self.canopy_factor_unit = model_input.open_file(self.canopy_factor_type,
                                                        self.canopy_factor_file)
        #----------------------------------------------------------------------------
        # Note: GMT_offset plus slope and aspect grids will be read separately.
        #----------------------------------------------------------------------------

##        self.Qn_SW_unit  = model_input.open_file(self.Qn_SW_type,  self.Qn_SW_file)
##        self.Qn_LW_unit  = model_input.open_file(self.Qn_LW_type,  self.Qn_LW_file)
        
    #   open_input_files()
    #-------------------------------------------------------------------  
    def read_input_files(self):

        if (self.DEBUG):
            print 'Calling read_input_files()...'

        rti = self.rti

        #--------------------------------------------------------
        # All grids are assumed to have a data type of Float32.
        #--------------------------------------------------------
        # NB! read_next() returns None if TYPE arg is "Scalar".
        #--------------------------------------------------------
        P = model_input.read_next(self.P_unit, self.P_type, rti,
                                  factor=self.mmph_to_mps)

        ## print '######### self.P_type = ' + self.P_type
        ## print '######### np.ndim( P ) = ' + str(np.ndim(P))

        if (P is not None):
            ## print 'MET: (time,P) =', self.time, P
            
            self.update_var( 'P', P )  ### 11/15/16

            ## if (self.P_type.lower() != 'scalar'):
#             if (np.ndim( self.P ) == 0):
#                 self.P.fill( P )  #### (2/7/13, mutable scalar)
#             else:
#                 self.P = P

            if (self.DEBUG or (self.time_index == 0)):
                print 'In read_input_files():'
                print '   min(P) =', P.min() * self.mps_to_mmph, ' [mmph]'
                print '   max(P) =', P.max() * self.mps_to_mmph, ' [mmph]'
                print ' '
        else:
            #-----------------------------------------------
            # Either self.P_type is "Scalar" or we've read
            # all of the data in the rain_rates file.
            #-----------------------------------------------
            if (self.P_type.lower() != 'scalar'):
                #------------------------------------
                # Precip is unique in this respect.
                #--------------------------------------------------
                # 2/7/13. Note that we don't change P from grid
                # to scalar since that could cause trouble for
                # other comps that use P, so we just zero it out.
                #--------------------------------------------------
                self.P.fill( 0 )
                if (self.DEBUG):
                    print 'Reached end of file:', self.P_file
                    print '  P set to 0 by read_input_files().'
            elif (self.time_sec >= self.dt):
                self.P.fill( 0 )
                if (self.DEBUG):
                    print 'Reached end of scalar rainfall duration.'
                    print '  P set to 0 by read_input_files().'
                    ## print 'time_sec =', self.time_sec
                    ## print 'met dt   =', self.dt
                    
##        print '######### In met_base.read_input_files() #######'
##        print 'self.P_type =', self.P_type
##        print 'self.P      =', self.P

        #------------------------------------------------------------
        # Read variables from files into scalars or grids while
        # making sure to preserve references (in-place). (11/15/16)
        #------------------------------------------------------------
        model_input.read_next2(self, 'T_air',   rti)
        model_input.read_next2(self, 'T_surf',  rti)
        model_input.read_next2(self, 'RH',      rti)
        model_input.read_next2(self, 'p0',      rti)
        model_input.read_next2(self, 'uz',      rti)
        model_input.read_next2(self, 'z',       rti)
        model_input.read_next2(self, 'z0_air',  rti)
        #----------------------------------------------------
        model_input.read_next2(self, 'albedo',  rti)
        model_input.read_next2(self, 'em_surf', rti)
        model_input.read_next2(self, 'dust_atten',    rti)
        model_input.read_next2(self, 'cloud_factor',  rti)
        model_input.read_next2(self, 'canopy_factor', rti)

        ###############################################################
        # If any of these are scalars (read from a time series file)
        # then we'll need to use "fill()" method to prevent breaking
        # the reference to the "mutable scalar". (2/7/13)
        ###############################################################
        T_air = model_input.read_next(self.T_air_unit, self.T_air_type, rti)
        self.update_var( 'T_air', T_air )
        # if (T_air is not None): self.T_air = T_air

        T_surf = model_input.read_next(self.T_surf_unit, self.T_surf_type, rti)
        self.update_var( 'T_surf', T_surf )
        # if (T_surf is not None): self.T_surf = T_surf

        RH = model_input.read_next(self.RH_unit, self.RH_type, rti)
        self.update_var( 'RH', RH )
        # if (RH is not None): self.RH = RH

        p0 = model_input.read_next(self.p0_unit, self.p0_type, rti)
        self.update_var( 'p0', p0 )
        # if (p0 is not None): self.p0 = p0

        uz = model_input.read_next(self.uz_unit, self.uz_type, rti)
        self.update_var( 'uz', uz )
        # if (uz is not None): self.uz = uz

        z = model_input.read_next(self.z_unit, self.z_type, rti)
        self.update_var( 'z', z )
        # if (z is not None): self.z = z

        z0_air = model_input.read_next(self.z0_air_unit, self.z0_air_type, rti)
        self.update_var( 'z0_air', z0_air )
        # if (z0_air is not None): self.z0_air = z0_air

        #----------------------------------------------------------------------------
        # These are needed to compute Qn_SW and Qn_LW.
        #----------------------------------------------------------------------------
        # Note: We could later write a version of read_next() that takes "self"
        #       and "var_name" as args and that uses "exec()".
        #----------------------------------------------------------------------------
        albedo = model_input.read_next(self.albedo_unit, self.albedo_type, rti)
        if (albedo is not None): self.albedo = albedo

        em_surf = model_input.read_next(self.em_surf_unit, self.em_surf_type, rti)
        if (em_surf is not None): self.em_surf = em_surf

        dust_atten = model_input.read_next(self.dust_atten_unit, self.dust_atten_type, rti)
        if (dust_atten is not None): self.dust_atten = dust_atten

        cloud_factor = model_input.read_next(self.cloud_factor_unit, self.cloud_factor_type, rti)
        if (cloud_factor is not None): self.cloud_factor = cloud_factor

        canopy_factor = model_input.read_next(self.canopy_factor_unit, self.canopy_factor_type, rti)
        if (canopy_factor is not None): self.canopy_factor = canopy_factor

        #-------------------------------------------------------------
        # Compute Qsw_prefactor from cloud_factor and canopy factor.
        #-------------------------------------------------------------
        ## self.Qsw_prefactor = 
        
        #-------------------------------------------------------------
        # These are currently treated as input data, but are usually
        # generated by functions in Qnet_file.py.  Later on, we'll
        # provide the option to compute them "on the fly" with new
        # functions called "update_net_shortwave_radiation()" and
        # "update_net_longwave_radiation()", called from update().
        #-------------------------------------------------------------        
##        Qn_SW = model_input.read_next(self.Qn_SW_unit, self.Qn_SW_type, rti)
##        if (Qn_SW is not None): self.Qn_SW = Qn_SW
##
##        Qn_LW = model_input.read_next(self.Qn_LW_unit, self.Qn_LW_type, rti)
##        if (Qn_LW is not None): self.Qn_LW = Qn_LW
         
    #   read_input_files()
    #-------------------------------------------------------------------  
    def close_input_files(self):

        if (self.DEBUG):
            print 'Calling close_input_files()...'

        if (self.P_type      != 'Scalar'): self.P_unit.close()
        if (self.T_air_type  != 'Scalar'): self.T_air_unit.close()
        if (self.T_surf_type != 'Scalar'): self.T_surf_unit.close()
        if (self.RH_type     != 'Scalar'): self.RH_unit.close()
        if (self.p0_type     != 'Scalar'): self.p0_unit.close()
        if (self.uz_type     != 'Scalar'): self.uz_unit.close()
        if (self.z_type      != 'Scalar'): self.z_unit.close()
        if (self.z0_air_type != 'Scalar'): self.z0_air_unit.close()

        #---------------------------------------------------
        # These are needed to compute Qn_SW and Qn_LW.
        #---------------------------------------------------
        if (self.albedo_type        != 'Scalar'): self.albedo_unit.close()
        if (self.em_surf_type       != 'Scalar'): self.em_surf_unit.close()
        if (self.dust_atten_type    != 'Scalar'): self.dust_atten_unit.close()
        if (self.cloud_factor_type  != 'Scalar'): self.cloud_factor_unit.close()
        if (self.canopy_factor_type != 'Scalar'): self.canopy_factor_unit.close()
        
##        if (self.Qn_SW_type  != 'Scalar'): self.Qn_SW_unit.close()        
##        if (self.Qn_LW_type  != 'Scalar'): self.Qn_LW_unit.close()
        
##        if (self.P_file      != ''): self.P_unit.close()
##        if (self.T_air_file  != ''): self.T_air_unit.close()
##        if (self.T_surf_file != ''): self.T_surf_unit.close()
##        if (self.RH_file     != ''): self.RH_unit.close()
##        if (self.p0_file     != ''): self.p0_unit.close()
##        if (self.uz_file     != ''): self.uz_unit.close()
##        if (self.z_file      != ''): self.z_unit.close()
##        if (self.z0_air_file != ''): self.z0_air_unit.close()
##        #--------------------------------------------------------
##        if (self.Qn_SW_file  != ''): self.Qn_SW_unit.close()        
##        if (self.Qn_LW_file  != ''): self.Qn_LW_unit.close()
        
    #   close_input_files()
    #-------------------------------------------------------------------  
    def update_outfile_names(self):

        if (self.DEBUG):
            print 'Calling update_outfile_names()...'
        
        #-------------------------------------------------
        # Notes:  Append out_directory to outfile names.
        #-------------------------------------------------
        self.ea_gs_file  = (self.out_directory + self.ea_gs_file  )
        self.es_gs_file  = (self.out_directory + self.es_gs_file  )
        self.Qsw_gs_file = (self.out_directory + self.Qsw_gs_file )
        self.Qlw_gs_file = (self.out_directory + self.Qlw_gs_file )
        self.ema_gs_file = (self.out_directory + self.ema_gs_file )
        #------------------------------------------------------------
        self.ea_ts_file  = (self.out_directory + self.ea_ts_file  )
        self.es_ts_file  = (self.out_directory + self.es_ts_file  )
        self.Qsw_ts_file = (self.out_directory + self.Qsw_ts_file )
        self.Qlw_ts_file = (self.out_directory + self.Qlw_ts_file )
        self.ema_ts_file = (self.out_directory + self.ema_ts_file )
        
##        self.ea_gs_file = (self.case_prefix + '_2D-ea.rts')
##        self.es_gs_file = (self.case_prefix + '_2D-es.rts')
##        #-----------------------------------------------------
##        self.ea_ts_file = (self.case_prefix + '_0D-ea.txt')
##        self.es_ts_file = (self.case_prefix + '_0D-es.txt')

    #   update_outfile_names()   
    #-------------------------------------------------------------------  
    def open_output_files(self):

        if (self.DEBUG):
            print 'Calling open_output_files()...'
        model_output.check_netcdf()
        self.update_outfile_names()
        
        #--------------------------------------
        # Open new files to write grid stacks
        #--------------------------------------
        if (self.SAVE_EA_GRIDS):
            model_output.open_new_gs_file( self, self.ea_gs_file, self.rti,
                                           ## var_name='e_air',
                                           var_name='ea',
                                           long_name='vapor_pressure_in_air',
                                           units_name='mbar')
            
        if (self.SAVE_ES_GRIDS):
            model_output.open_new_gs_file( self, self.es_gs_file, self.rti,
                                           ## var_name='e_surf',
                                           var_name='es',
                                           long_name='vapor_pressure_at_surface',
                                           units_name='mbar')
        if (self.SAVE_QSW_GRIDS):
            model_output.open_new_gs_file( self, self.Qsw_gs_file, self.rti,
                                           var_name='Qsw',
                                           long_name='net_shortwave_radiation',
                                           units_name='W/m^2')
            
        if (self.SAVE_QLW_GRIDS):
            model_output.open_new_gs_file( self, self.Qlw_gs_file, self.rti,
                                           var_name='Qlw',
                                           long_name='net_longwave_radiation',
                                           units_name='W/m^2')
        if (self.SAVE_EMA_GRIDS):
            model_output.open_new_gs_file( self, self.ema_gs_file, self.rti,
                                           var_name='ema',
                                           long_name='air_emissivity',
                                           units_name='none')
            
        #--------------------------------------
        # Open new files to write time series
        #--------------------------------------
        IDs = self.outlet_IDs 
        if (self.SAVE_EA_PIXELS):
            model_output.open_new_ts_file( self, self.ea_ts_file, IDs,
                                           ## var_name='e_air',
                                           var_name='ea',
                                           long_name='vapor_pressure_in_air',
                                           units_name='mbar')

        if (self.SAVE_ES_PIXELS):
            model_output.open_new_ts_file( self, self.es_ts_file, IDs,
                                           ## var_name='e_surf',
                                           var_name='es',
                                           long_name='vapor_pressure_at_surface',
                                           units_name='mbar')

        if (self.SAVE_QSW_PIXELS):
            model_output.open_new_ts_file( self, self.Qsw_ts_file, IDs,
                                           var_name='Qsw',
                                           long_name='net_shortwave_radiation',
                                           units_name='W/m^2')

        if (self.SAVE_QLW_PIXELS):
            model_output.open_new_ts_file( self, self.Qlw_ts_file, IDs,
                                           var_name='Qlw',
                                           long_name='net_longwave_radiation',
                                           units_name='W/m^2')
            
        if (self.SAVE_EMA_PIXELS):
            model_output.open_new_ts_file( self, self.ema_ts_file, IDs,
                                           var_name='ema',
                                           long_name='air_emissivity',
                                           units_name='none')
            
    #   open_output_files()
    #-------------------------------------------------------------------
    def write_output_files(self, time_seconds=None):

        if (self.DEBUG):
            print 'Calling write_output_files()...'

        #-----------------------------------------
        # Allows time to be passed from a caller
        #-----------------------------------------
        if (time_seconds is None):
            time_seconds = self.time_sec
        model_time = int(time_seconds)
        
        #----------------------------------------
        # Save computed values at sampled times
        #----------------------------------------
        if (model_time % int(self.save_grid_dt) == 0):
            self.save_grids()
        if (model_time % int(self.save_pixels_dt) == 0):
            self.save_pixel_values()

    #  write_output_files()
    #-------------------------------------------------------------------
    def close_output_files(self):
    
        if (self.SAVE_EA_GRIDS):   model_output.close_gs_file( self, 'ea')
        if (self.SAVE_ES_GRIDS):   model_output.close_gs_file( self, 'es')
        if (self.SAVE_QSW_GRIDS):  model_output.close_gs_file( self, 'Qsw')
        if (self.SAVE_QLW_GRIDS):  model_output.close_gs_file( self, 'Qlw')
        if (self.SAVE_EMA_GRIDS):  model_output.close_gs_file( self, 'ema')
        #-------------------------------------------------------------------
        if (self.SAVE_EA_PIXELS):  model_output.close_ts_file( self, 'ea') 
        if (self.SAVE_ES_PIXELS):  model_output.close_ts_file( self, 'es') 
        if (self.SAVE_QSW_PIXELS): model_output.close_ts_file( self, 'Qsw') 
        if (self.SAVE_QLW_PIXELS): model_output.close_ts_file( self, 'Qlw')
        if (self.SAVE_EMA_PIXELS): model_output.close_ts_file( self, 'ema')
        
    #   close_output_files()        
    #-------------------------------------------------------------------  
    def save_grids(self):
       
        if (self.SAVE_EA_GRIDS):
            model_output.add_grid( self, self.e_air,  'ea', self.time_min )
            
        if (self.SAVE_ES_GRIDS):
            model_output.add_grid( self, self.e_surf, 'es', self.time_min )

        if (self.SAVE_QSW_GRIDS):
            model_output.add_grid( self, self.Qn_SW, 'Qsw', self.time_min )
            
        if (self.SAVE_QLW_GRIDS):
            model_output.add_grid( self, self.Qn_LW, 'Qlw', self.time_min )

        if (self.SAVE_EMA_GRIDS):
            model_output.add_grid( self, self.em_air, 'ema', self.time_min )
            
    #   save_grids()            
    #-------------------------------------------------------------------  
    def save_pixel_values(self):

        IDs  = self.outlet_IDs
        time = self.time_min   ######
        
        if (self.SAVE_EA_PIXELS):
            model_output.add_values_at_IDs( self, time, self.e_air,  'ea', IDs )
            
        if (self.SAVE_ES_PIXELS):
            model_output.add_values_at_IDs( self, time, self.e_surf, 'es', IDs )

        if (self.SAVE_QSW_PIXELS):
            model_output.add_values_at_IDs( self, time, self.Qn_SW, 'Qsw', IDs )
            
        if (self.SAVE_QLW_PIXELS):
            model_output.add_values_at_IDs( self, time, self.Qn_LW, 'Qlw', IDs )

        if (self.SAVE_EMA_PIXELS):
            model_output.add_values_at_IDs( self, time, self.em_air, 'ema', IDs )
            
    #   save_pixel_values()
    #-------------------------------------------------------------------
#---------------------------------------------------------------------------------
def compare_em_air_methods():

    #--------------------------------------------------------------
    # Notes:  There are two different methods that are commonly
    #         used to compute the vapor pressure of air, e_air,
    #         and then the emissivity of air, em_air, for use in
    #         longwave radiation calculations.  This routine
    #         compares them graphically.
    #
    # NB!     This hasn't been tested since conversion from IDL.
    #-------------------------------------------------------------
    import matplotlib.pyplot
    
    T_air = np.arange(80, dtype='Float32') - np.float64(40)   #[Celsius]  (-40 to 40)
    RH  = np.float64(1.0)
    C2K = np.float64(273.15)
    
    #--------------------------
    # Brutsaert (1975) method
    #--------------------------
    term1   = (np.float64(17.3) * T_air) / (T_air + np.float64(237.3))   ######### DOUBLE CHECK THIS (7/26/13)
    e_air1  = RH * np.float64(0.611) * np.exp( term1 )  # [kPa]
    em_air1 = np.float64(1.72) * (e_air1 / (T_air + C2K)) ** (np.float64(1) / 7)
    
    #---------------------------
    # Satterlund (1979) method
    #----------------------------
    # NB! e_air has units of Pa
    #----------------------------
    term2   = np.float64(2353) / (T_air + C2K)
    e_air2  = RH * np.float64(10) ** (np.float64(11.40) - term2)   # [Pa]
    eterm   = np.exp(-np.float64(1) * (e_air2 / np.float64(100)) ** ((T_air + C2K) / np.float64(2016)))
    em_air2 = np.float64(1.08) * (np.float64(1) - eterm)
    
    #----------------------------
    # Plot the two e_air curves
    #--------------------------------
    # These two agree quite closely
    #--------------------------------
    matplotlib.pyplot.figure(figsize=(8, 6), dpi=80)
    matplotlib.pyplot.show()
    matplotlib.pyplot.plot(T_air, e_air1)
    matplotlib.pyplot.show()
    ## oplot(T_air, (e_air2 / np.float64(1000)), psym=-3)   # [Pa -> kPa]
    
    #-----------------------------
    # Plot the two em_air curves
    #--------------------------------------------------
    # These two don't agree very well for some reason
    #--------------------------------------------------
    matplotlib.pyplot.figure(figsize=(8, 6), dpi=80)
    matplotlib.pyplot.show()
    matplotlib.pyplot.plot(T_air, em_air1)
    matplotlib.pyplot.show()
    ## oplot(T_air, em_air2, psym=-3)
    
#   compare_em_air_Methods
#---------------------------------------------------------------------------------
    
        
