diff --git a/CHANGELOG.md b/CHANGELOG.md
index fcd99f784..db82e7073 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -6,7 +6,11 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [PEP 440](https://www.python.org/dev/peps/pep-0440/)
 and uses [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [0.4.3]
++ Temporal interpolation of delays if the requested datetime is more than _THRESHOLD_SECONDS away from the closest weather model available time and `interpolate_time = True` (default behavior)
 + Add assert statement to raise error if the delay cube for each SAR date in a GUNW IFG is not written 
++ Verify some constants / equations and remove the comments questioning them
++ Relocate the time resolution of wmodels to one spot
 
 ## [0.4.2]
 
diff --git a/test/test_temporal_interpolate.py b/test/test_temporal_interpolate.py
new file mode 100644
index 000000000..6bd93d4bf
--- /dev/null
+++ b/test/test_temporal_interpolate.py
@@ -0,0 +1,193 @@
+import glob
+import os
+import pytest
+import shutil
+import subprocess
+import yaml
+
+import numpy as np
+import pandas as pd
+import xarray as xr
+
+from test import TEST_DIR, WM
+
+import RAiDER
+from RAiDER.logger import logger
+
+
+def makeLatLonGrid(bbox, reg, out_dir, spacing=0.1):
+    """ Make lat lons at a specified spacing """
+    S, N, W, E = bbox
+    lat_st, lat_en = S, N
+    lon_st, lon_en = W, E
+
+    lats = np.arange(lat_st, lat_en, spacing)
+    lons = np.arange(lon_st, lon_en, spacing)
+    Lat, Lon = np.meshgrid(lats, lons)
+    da_lat = xr.DataArray(Lat.T, name='data', coords={'lon': lons, 'lat': lats}, dims='lat lon'.split())
+    da_lon = xr.DataArray(Lon.T, name='data', coords={'lon': lons, 'lat': lats}, dims='lat lon'.split())
+
+    dst_lat = os.path.join(out_dir, f'lat_{reg}.nc')
+    dst_lon = os.path.join(out_dir, f'lon_{reg}.nc')
+    da_lat.to_netcdf(dst_lat)
+    da_lon.to_netcdf(dst_lon)
+
+    return dst_lat, dst_lon
+
+
+def update_yaml(dct_cfg:dict, dst:str='temp.yaml'):
+    """ Write a new yaml file from a dictionary.
+
+    Updates parameters in the default 'raider.yaml' file.
+    Each key:value pair will in 'dct_cfg' will overwrite that in the default
+    """
+
+    template_file = os.path.join(
+                    os.path.dirname(RAiDER.__file__), 'cli', 'raider.yaml')
+
+    with open(template_file, 'r') as f:
+        try:
+            params = yaml.safe_load(f)
+        except yaml.YAMLError as exc:
+            print(exc)
+            raise ValueError(f'Something is wrong with the yaml file {template_file}')
+
+    params = {**params, **dct_cfg}
+
+    with open(dst, 'w') as fh:
+        yaml.safe_dump(params, fh,  default_flow_style=False)
+
+    return dst
+
+
+def test_cube_timemean():
+    """ Test the mean interpolation by computing cube delays at 1:30PM vs mean of 12 PM / 3PM for GMAO """
+    SCENARIO_DIR = os.path.join(TEST_DIR, "INTERP_TIME")
+    os.makedirs(SCENARIO_DIR, exist_ok=True)
+    ## make the lat lon grid
+    S, N, W, E = 34, 35, -117, -116
+    date       = 20200130
+    dct_hrs    = {'GMAO': [12, 15, '13:30:00'], 'ERA5': [13, 14, '13:30:00']}
+    hr1, hr2, ti = dct_hrs[WM]
+
+    grp = {
+            'date_group': {'date_start': date},
+            'weather_model': WM,
+            'aoi_group': {'bounding_box': [S, N, W, E]},
+            'runtime_group': {'output_directory': SCENARIO_DIR},
+        }
+
+
+    ## run raider original for two exact weather model times
+    for hr in [hr1, hr2]:
+        grp['time_group'] =  {'time': f'{hr}:00:00'}
+        ## generate the default template file and overwrite it with new parms
+        cfg  = update_yaml(grp)
+
+        ## run raider for the default date
+        cmd  = f'raider.py {cfg}'
+        proc = subprocess.run(cmd.split(), stdout=subprocess.PIPE, universal_newlines=True)
+        assert np.isclose(proc.returncode, 0)
+
+    ## run interpolation in the middle of the two
+    grp['time_group'] =  {'time': ti, 'interpolate_time': True}
+    cfg  = update_yaml(grp)
+
+    cmd  = f'raider.py {cfg}'
+    proc = subprocess.run(cmd.split(), stdout=subprocess.PIPE, universal_newlines=True)
+    assert np.isclose(proc.returncode, 0)
+
+
+    with xr.open_dataset(os.path.join(SCENARIO_DIR, f'{WM}_tropo_{date}T{hr1}0000_ztd.nc')) as ds:
+        da1_tot = ds['wet'] + ds['hydro']
+
+    with xr.open_dataset(os.path.join(SCENARIO_DIR, f'{WM}_tropo_{date}T{hr2}0000_ztd.nc')) as ds:
+        da2_tot = ds['wet'] + ds['hydro']
+
+    with xr.open_dataset(os.path.join(SCENARIO_DIR, f'{WM}_tropo_{date}T{ti.replace(":", "")}_ztd.nc')) as ds:
+        da_interp_tot = ds['wet'] + ds['hydro']
+
+    da_mu = (da1_tot + da2_tot) / 2
+    assert np.allclose(da_mu, da_interp_tot)
+
+
+    # Clean up files
+    shutil.rmtree(SCENARIO_DIR)
+    [os.remove(f) for f in glob.glob(f'{WM}*')]
+    os.remove('temp.yaml')
+
+    return
+
+
+def test_cube_weighting():
+    """ Test the weighting by comparing a small crop with numpy directly """
+    from datetime import datetime
+    SCENARIO_DIR = os.path.join(TEST_DIR, "INTERP_TIME")
+    os.makedirs(SCENARIO_DIR, exist_ok=True)
+    ## make the lat lon grid
+    S, N, W, E = 34, 35, -117, -116
+    date       = 20200130
+    dct_hrs    = {'GMAO': [12, 15, '12:05:00'], 'ERA5': [13, 14, '13:05:00']}
+    hr1, hr2, ti = dct_hrs[WM]
+
+    grp = {
+            'date_group': {'date_start': date},
+            'weather_model': WM,
+            'aoi_group': {'bounding_box': [S, N, W, E]},
+            'runtime_group': {'output_directory': SCENARIO_DIR},
+        }
+
+
+    ## run raider original for two exact weather model times
+    for hr in [hr1, hr2]:
+        grp['time_group'] =  {'time': f'{hr}:00:00'}
+        ## generate the default template file and overwrite it with new parms
+        cfg  = update_yaml(grp)
+
+        ## run raider for the default date
+        cmd  = f'raider.py {cfg}'
+        proc = subprocess.run(cmd.split(), stdout=subprocess.PIPE, universal_newlines=True)
+        assert np.isclose(proc.returncode, 0)
+
+    ## run interpolation very near the first
+    grp['time_group'] =  {'time': ti, 'interpolate_time': True}
+    cfg  = update_yaml(grp)
+
+    cmd  = f'raider.py {cfg}'
+    proc = subprocess.run(cmd.split(), stdout=subprocess.PIPE, universal_newlines=True)
+
+    ## double check on weighting
+
+    with xr.open_dataset(os.path.join(SCENARIO_DIR, f'{WM}_tropo_{date}T{hr1}0000_ztd.nc')) as ds:
+        da1_tot = ds['wet'] + ds['hydro']
+
+    with xr.open_dataset(os.path.join(SCENARIO_DIR, f'{WM}_tropo_{date}T{hr2}0000_ztd.nc')) as ds:
+        da2_tot = ds['wet'] + ds['hydro']
+
+    with xr.open_dataset(os.path.join(SCENARIO_DIR, f'{WM}_tropo_{date}T{ti.replace(":", "")}_ztd.nc')) as ds:
+        da_interp_tot = ds['wet'] + ds['hydro']
+
+    dt1 = datetime.strptime(f'{date}{hr1}', '%Y%m%d%H')
+    dt2 = datetime.strptime(f'{date}{hr2}', '%Y%m%d%H')
+    dt_ref = datetime.strptime(f'{date}{ti}', '%Y%m%d%H:%M:%S')
+
+    wgts  = np.array([(dt_ref-dt1).seconds, (dt2-dt_ref).seconds])
+    da1_crop = da1_tot.isel(z=0, y=slice(0,1), x=slice(0, 2))
+    da2_crop = da2_tot.isel(z=0, y=slice(0,1), x=slice(0, 2))
+    da_out_crop = da_interp_tot.isel(z=0, y=slice(0,1), x=slice(0,2))
+
+    dat = np.vstack([da1_crop.data, da2_crop.data])
+
+    logger.info ('Tstart: %s, Tend: %s, Tref: %s', dt1, dt2, dt_ref)
+    logger.info ('Weights: %s', wgts)
+    logger.info ('Data from two dates: %s', dat)
+    logger.info ('Weighted mean: %s', da_out_crop.data)
+    assert np.allclose(da_out_crop, np.average(dat, weights=1/wgts, axis=0))
+
+    # Clean up files
+    shutil.rmtree(SCENARIO_DIR)
+    [os.remove(f) for f in glob.glob(f'{WM}*')]
+    os.remove('temp.yaml')
+
+    return
+
diff --git a/test/test_util.py b/test/test_util.py
index 176c654eb..5539fbb70 100644
--- a/test/test_util.py
+++ b/test/test_util.py
@@ -13,7 +13,7 @@
     _least_nonzero, cosd, rio_open, sind,
     writeArrayToRaster, rio_profile,
     rio_extents, getTimeFromFile, enu2ecef, ecef2enu,
-    transform_bbox, clip_bbox
+    transform_bbox, clip_bbox, get_nearest_wmtimes,
 )
 
 
@@ -462,3 +462,27 @@ def test_clip_bbox():
     snwe_in = [34.005, 35.0006, -76.999, -76.0]
     assert clip_bbox(wesn, 0.01) == wesn
     assert clip_bbox(snwe_in, 0.01) == snwe
+
+def test_get_nearest_wmtimes():
+    t0 = datetime.datetime(2020,1,1,11,35,0)
+    test_out = get_nearest_wmtimes(t0, 3)
+    true_out = [datetime.datetime(2020, 1, 1, 9, 0), datetime.datetime(2020, 1, 1, 12, 0)]
+    assert [t == t0 for t, t0 in zip(test_out, true_out)]
+
+def test_get_nearest_wmtimes_2():
+    t0 = datetime.datetime(2020,1,1,11,3,0)
+    test_out = get_nearest_wmtimes(t0, 1)
+    true_out = [datetime.datetime(2020, 1, 1, 11, 0)]
+    assert [t == t0 for t, t0 in zip(test_out, true_out)]
+
+def test_get_nearest_wmtimes_3():
+    t0 = datetime.datetime(2020,1,1,11,57,0)
+    test_out = get_nearest_wmtimes(t0, 3)
+    true_out = [datetime.datetime(2020, 1, 1, 12, 0)]
+    assert [t == t0 for t, t0 in zip(test_out, true_out)]
+
+def test_get_nearest_wmtimes_4():
+    t0 = datetime.datetime(2020,1,1,11,25,0)
+    test_out = get_nearest_wmtimes(t0, 1)
+    true_out = [datetime.datetime(2020, 1, 1, 11, 0), datetime.datetime(2020, 1, 1, 12, 0)]
+    assert [t == t0 for t, t0 in zip(test_out, true_out)]
diff --git a/tools/RAiDER/aria/prepFromGUNW.py b/tools/RAiDER/aria/prepFromGUNW.py
index e50161fe7..ad4d8b3d7 100644
--- a/tools/RAiDER/aria/prepFromGUNW.py
+++ b/tools/RAiDER/aria/prepFromGUNW.py
@@ -244,8 +244,8 @@ def main(args):
            'cube_spacing_in_m': GUNWObj.spacing_m,
            'aoi_group' : {'bounding_box': GUNWObj.SNWE},
            'height_group' : {'height_levels': GUNWObj.heights},
-           'date_group': {'date_list': str(GUNWObj.dates)},
-           'time_group': {'time': GUNWObj.ref_time},
+           'date_group': {'date_list': GUNWObj.dates},
+           'time_group': {'time': GUNWObj.mid_time, 'interpolate_time': True},
            'los_group' : {'ray_trace': True,
                           'orbit_file': GUNWObj.OrbitFiles,
                           'wavelength': GUNWObj.wavelength,
diff --git a/tools/RAiDER/cli/__init__.py b/tools/RAiDER/cli/__init__.py
index 8a0715ae1..3b6156ef2 100644
--- a/tools/RAiDER/cli/__init__.py
+++ b/tools/RAiDER/cli/__init__.py
@@ -39,5 +39,6 @@ class AttributeDict(dict):
             output_directory='.',
             weather_model_directory=None,
             output_projection='EPSG:4326',
+            interpolate_time=True,
         )
     )
diff --git a/tools/RAiDER/cli/raider.py b/tools/RAiDER/cli/raider.py
index 3f471a867..cc4785a27 100644
--- a/tools/RAiDER/cli/raider.py
+++ b/tools/RAiDER/cli/raider.py
@@ -1,9 +1,13 @@
 import argparse
+import datetime
 import os
 import shutil
 import sys
 import yaml
 
+import numpy as np
+import xarray as xr
+
 from textwrap import dedent
 
 from RAiDER import aws
@@ -12,6 +16,7 @@
 from RAiDER.cli.parser import add_out, add_cpus, add_verbose
 from RAiDER.cli.validators import DateListAction, date_type
 from RAiDER.models.allowed import ALLOWED_MODELS
+from RAiDER.utilFcns import get_dt
 
 
 HELP_MESSAGE = """
@@ -124,7 +129,7 @@ def calcDelays(iargs=None):
     from RAiDER.delay import tropo_delay
     from RAiDER.checkArgs import checkArgs
     from RAiDER.processWM import prepareWeatherModel
-    from RAiDER.utilFcns import writeDelays
+    from RAiDER.utilFcns import writeDelays, get_nearest_wmtimes
     examples = 'Examples of use:' \
         '\n\t raider.py customTemplatefile.cfg' \
         '\n\t raider.py -g'
@@ -218,21 +223,58 @@ def calcDelays(iargs=None):
         logger.debug('Starting to run the weather model calculation')
         logger.debug(f'Date: {t.strftime("%Y%m%d")}')
         logger.debug('Beginning weather model pre-processing')
-        try:
-            weather_model_file = prepareWeatherModel(
-                model, t,
-                ll_bounds=ll_bounds, # SNWE
-                wmLoc=params['weather_model_directory'],
-                makePlots=params['verbose'],
-            )
-        except RuntimeError:
-            logger.exception("Date %s failed", t)
-            continue
+
+        # Grab the closest two times unless the user specifies 'nearest'
+        # If the model time_delta is not specified then use 6
+        # The two datetimes will be combined to a single file and processed
+        times = get_nearest_wmtimes(t, [model.dtime() if model.dtime() is not None else 6][0]) if params['interpolate_time'] else [t]
+        wfiles = []
+        for tt in times:
+            try:
+                wfiles.append(
+                    prepareWeatherModel(
+                        model, tt,
+                        ll_bounds=ll_bounds, # SNWE
+                        wmLoc=params['weather_model_directory'],
+                        makePlots=params['verbose'],
+                    )
+                )
+            except RuntimeError:
+                logger.exception("Date %s failed", t)
+                continue
 
         # dont process the delays for download only
         if dl_only:
             continue
 
+        if len(wfiles)>1:
+            ds1 = xr.open_dataset(wfiles[0])
+            ds2 = xr.open_dataset(wfiles[1])
+            
+            # calculate relative weights of each dataset
+            date1 = datetime.datetime.strptime(ds1.attrs['datetime'], '%Y_%m_%dT%H_%M_%S')
+            date2 = datetime.datetime.strptime(ds2.attrs['datetime'], '%Y_%m_%dT%H_%M_%S')
+            wgts  = [ 1 - get_dt(t, date1) / get_dt(date2, date1), 1 - get_dt(date2, t) / get_dt(date2, date1)]
+            try:
+                assert np.sum(wgts)==1
+            except AssertionError:
+                logging.error('Time interpolation weights do not sum to one, something is off with the dates')
+                raise ValueError('Time interpolation weights do not sum to one, something is off with the dates')
+
+            # combine datasets
+            ds = ds1
+            for var in ['wet', 'hydro', 'wet_total', 'hydro_total']:
+                ds[var] = (wgts[0] * ds1[var]) + (wgts[1] * ds2[var])
+            ds.attrs['Date1'] = 0
+            ds.attrs['Date2'] = 0
+            weather_model_file = os.path.join(
+                os.path.dirname(wfiles[0]),
+                os.path.basename(wfiles[0]).split('_')[0] + '_' + t.strftime('%Y_%m_%dT%H_%M_%S') + '_timeInterp_' + '_'.join(wfiles[0].split('_')[-4:]),
+            )
+            ds.to_netcdf(weather_model_file)
+        else:
+            weather_model_file = wfiles[0]
+
         # Now process the delays
         try:
             wet_delay, hydro_delay = tropo_delay(
@@ -246,10 +288,7 @@ def calcDelays(iargs=None):
             logger.exception("Date %s failed", t)
             continue
 
-        ###########################################################
-        # Write the delays to file
         # Different options depending on the inputs
-
         if los.is_Projected():
             out_filename = w.replace("_ztd", "_std")
             f = f.replace("_ztd", "_std")
diff --git a/tools/RAiDER/cli/raider.yaml b/tools/RAiDER/cli/raider.yaml
index 4c84e8724..33fc5be05 100644
--- a/tools/RAiDER/cli/raider.yaml
+++ b/tools/RAiDER/cli/raider.yaml
@@ -58,6 +58,7 @@ date_group:
 time_group:
   time:
   end_time:
+  interpolate_time: True # if True will interpolate two closest times, if False will take the nearest time
 
 
 ########## 4. Area of Interest
diff --git a/tools/RAiDER/cli/validators.py b/tools/RAiDER/cli/validators.py
index e07a778b1..b5d999148 100755
--- a/tools/RAiDER/cli/validators.py
+++ b/tools/RAiDER/cli/validators.py
@@ -114,7 +114,7 @@ def get_heights(args, out, station_file, bounding_box=None):
         out['height_levels'] = np.array([float(ll) for ll in l])
         if np.any(out['height_levels'] < 0):
             logger.warning('Weather model only extends to the surface topography; '
-            'height levels below the topography will be interpolated from the surface'
+            'height levels below the topography will be interpolated from the surface '
             'and may be inaccurate.')
 
     return out
@@ -202,6 +202,8 @@ def parse_dates(arg_dict):
         l = arg_dict['date_list']
         if isinstance(l, str):
             l = re.findall('[0-9]+', l)
+        elif isinstance(l, int):
+            l = [l]
         L = [enforce_valid_dates(d) for d in l]
 
     else:
diff --git a/tools/RAiDER/constants.py b/tools/RAiDER/constants.py
index 7d7c5cbe7..0b3308ea3 100644
--- a/tools/RAiDER/constants.py
+++ b/tools/RAiDER/constants.py
@@ -11,10 +11,14 @@
 _ZREF = np.float64(15000)  # maximum requierd height
 _STEP = np.float64(15.0)     # integration step size in meters
 
+_g0 = np.float64(9.80665) # Standard gravitational constant
+_g1 = np.float64(9.80616) # Gravitational constant @ 45° latitude used for corrections of earth's centrifugal force
+_RE = np.float64(6371008.7714)
+
+R_EARTH_MAX_WGS84 = 6378137
+R_EARTH_MIN_WGS84 = 6356752
+
 _CUBE_SPACING_IN_M = float(2000)  # Horizontal spacing of cube
+_THRESHOLD_SECONDS = 1 * 60 # Threshold delta_time in seconds
 
-_g0 = np.float64(9.80665)
-_RE = np.float64(6371008.7714)
 
-R_EARTH_MAX = 6378137
-R_EARTH_MIN = 6356752
diff --git a/tools/RAiDER/delay.py b/tools/RAiDER/delay.py
index da986ce92..18d354ad5 100755
--- a/tools/RAiDER/delay.py
+++ b/tools/RAiDER/delay.py
@@ -43,7 +43,7 @@ def tropo_delay(
         look_dir: str='right',
     ):
     """
-    Calculate integrated delays on query points. Options are: 
+    Calculate integrated delays on query points. Options are:
     1. Zenith delays (ZTD)
     2. Zenith delays projected to the line-of-sight (STD-projected)
     3. Slant delays integrated along the raypath (STD-raytracing)
@@ -70,7 +70,7 @@ def tropo_delay(
        wm_proj = CRS.from_epsg(4326)
     else:
         wm_proj = CRS.from_wkt(wm_proj.to_wkt())
-    
+
     # get heights
     if height_levels is None:
         if aoi.type() == 'Geocube':
@@ -128,7 +128,7 @@ def _get_delays_on_cube(dt, weather_model_file, wm_proj, ll_bounds, heights, los
     snwe = transform_bbox(ll_bounds, src_crs=4326, dest_crs=crs)
     out_spacing = get_output_spacing(cube_spacing_m, weather_model_file, wm_proj, crs)
     out_snwe = clip_bbox(snwe, out_spacing)
-    
+
     logger.debug(f"Output SNWE: {out_snwe}")
     logger.debug(f"Output cube spacing: {out_spacing}")
 
diff --git a/tools/RAiDER/models/ecmwf.py b/tools/RAiDER/models/ecmwf.py
index 532abdd93..6f134893b 100755
--- a/tools/RAiDER/models/ecmwf.py
+++ b/tools/RAiDER/models/ecmwf.py
@@ -14,8 +14,8 @@
     A_137_HRES,
     B_137_HRES,
 )
-from RAiDER.models.weatherModel import WeatherModel
 
+from RAiDER.models.weatherModel import WeatherModel, TIME_RES
 
 class ECMWF(WeatherModel):
     '''
@@ -31,7 +31,7 @@ def __init__(self):
         self._k2 = 0.233   # [K/Pa]
         self._k3 = 3.75e3  # [K^2/Pa]
 
-        self._time_res = 1
+        self._time_res = TIME_RES['ECMWF']
 
         self._lon_res = 0.2
         self._lat_res = 0.2
@@ -164,6 +164,8 @@ def _get_from_ecmwf(self, lat_min, lat_max, lat_step, lon_min, lon_max,
         server = ecmwfapi.ECMWFDataServer()
 
         corrected_date = util.round_date(time, datetime.timedelta(hours=6))
+        if not corrected_date == time:
+            logger.warning('Rounded given datetime from  %s to %s', time, corrected_date)
 
         server.retrieve({
             "class": self._classname,  # ERA-Interim
@@ -204,6 +206,7 @@ def _get_from_cds(
         acqTime,
         outname
     ):
+        """ Used for ERA5 """
         import cdsapi
         c = cdsapi.Client(verify=0)
 
@@ -247,7 +250,9 @@ def _get_from_cds(
             logger.exception(e)
             raise Exception
 
+
     def _download_ecmwf(self, lat_min, lat_max, lat_step, lon_min, lon_max, lon_step, time, out):
+        """ Used for HRES """
         from ecmwfapi import ECMWFService
 
         server = ECMWFService("mars")
@@ -255,7 +260,7 @@ def _download_ecmwf(self, lat_min, lat_max, lat_step, lon_min, lon_max, lon_step
         # round to the closest legal time
         corrected_date = util.round_date(time, datetime.timedelta(hours=self._time_res))
         if not corrected_date == time:
-            logger.warning('Rounded given datetime from  %s to %s', time, corrected_date)
+            logger.warning('Rounded given hour from  %d to %d', time.hour, corrected_date.hour)
 
         if self._model_level_type == 'ml':
             param = "129/130/133/152"
diff --git a/tools/RAiDER/models/gmao.py b/tools/RAiDER/models/gmao.py
index ad24655ed..e0c9d3fd0 100755
--- a/tools/RAiDER/models/gmao.py
+++ b/tools/RAiDER/models/gmao.py
@@ -7,7 +7,7 @@
 import pydap.client
 from pyproj import CRS
 
-from RAiDER.models.weatherModel import WeatherModel
+from RAiDER.models.weatherModel import WeatherModel, TIME_RES
 from RAiDER.logger import logger
 from RAiDER.utilFcns import writeWeatherVars2NETCDF4, round_time, requests_retry_session
 from RAiDER.models.model_levels import (
@@ -37,7 +37,8 @@ def __init__(self):
         self._k2 = 0.233  # [K/Pa]
         self._k3 = 3.75e3  # [K^2/Pa]
 
-        self._time_res = 1
+        self._time_res = TIME_RES[self._dataset.upper()]
+
 
         # horizontal grid spacing
         self._lat_res = 0.25
@@ -69,7 +70,7 @@ def _fetch(self, out):
         T0 = dt.datetime(2017, 12, 1, 0, 0, 0)
         # round time to nearest third hour
         time1 = time
-        time = round_time(time, 3 * 60 * 60)
+        time = round_time(time, self._time_res * 60 * 60)
         if not time1 == time:
             logger.warning('Rounded given hour from  %d to %d. May be incorrect near beginning/end of day.', time1.hour, time.hour)
 
diff --git a/tools/RAiDER/models/hres.py b/tools/RAiDER/models/hres.py
index 4dfc45dd8..97be42c87 100755
--- a/tools/RAiDER/models/hres.py
+++ b/tools/RAiDER/models/hres.py
@@ -5,7 +5,7 @@
 from pyproj import CRS
 
 from RAiDER.models.ecmwf import ECMWF
-from RAiDER.models.weatherModel import WeatherModel
+from RAiDER.models.weatherModel import WeatherModel, TIME_RES
 from RAiDER.models.model_levels import (
     LEVELS_91_HEIGHTS,
     LEVELS_25_HEIGHTS,
@@ -28,7 +28,6 @@ def __init__(self, level_type='ml'):
         self._k2 = 0.233   # [K/Pa]
         self._k3 = 3.75e3  # [K^2/Pa]
 
-        self._time_res = 6
 
         # 9 km horizontal grid spacing. This is only used for extending the download-buffer, i.e. not in subsequent processing.
         self._lon_res = 9. / 111  # 0.08108115
@@ -44,6 +43,7 @@ def __init__(self, level_type='ml'):
         self._Name = 'HRES'
         self._proj = CRS.from_epsg(4326)
 
+        self._time_res = TIME_RES[self._dataset.upper()]
         # Tuple of min/max years where data is available.
         self._valid_range = (datetime.datetime(1983, 4, 20), "Present")
         # Availability lag time in days
diff --git a/tools/RAiDER/models/hrrr.py b/tools/RAiDER/models/hrrr.py
index 88599bf93..5bcb7497e 100644
--- a/tools/RAiDER/models/hrrr.py
+++ b/tools/RAiDER/models/hrrr.py
@@ -12,8 +12,10 @@
 from pyproj import CRS, Transformer
 
 from RAiDER.logger import logger
-from RAiDER.utilFcns import rio_profile, rio_extents
-from RAiDER.models.weatherModel import WeatherModel, transform_coords
+from RAiDER.utilFcns import rio_profile, rio_extents, round_date
+from RAiDER.models.weatherModel import (
+    WeatherModel, transform_coords, TIME_RES
+)
 from RAiDER.models.model_levels import (
     LEVELS_137_HEIGHTS,
 )
@@ -30,7 +32,7 @@ def __init__(self):
         self._classname = 'hrrr'
         self._dataset = 'hrrr'
 
-        self._time_res = 1
+        self._time_res = TIME_RES[self._dataset.upper()]
 
         # Tuple of min/max years where data is available.
         self._valid_range = (datetime.datetime(2016, 7, 15), "Present")
@@ -262,8 +264,14 @@ def _download_hrrr_file(self, DATE, out, model='hrrr', product='prs', fxx=0, ver
         '''
         Download a HRRR model
         '''
+        ## TODO: Check how Herbie does temporal interpolation
+        # corrected_DT = round_date(DATE, datetime.timedelta(hours=self._time_res))
+        # if not corrected_DT == DATE:
+        #     logger.warning('Rounded given datetime from  %s to %s', DATE, corrected_DT)
+        corrected_DT = DATE
+
         H = Herbie(
-            DATE.strftime('%Y-%m-%d %H:%M'),
+            corrected_DT.strftime('%Y-%m-%d %H:%M'),
             model=model,
             product=product,
             overwrite=False,
diff --git a/tools/RAiDER/models/ncmr.py b/tools/RAiDER/models/ncmr.py
index 84fa62d35..86e3b418a 100755
--- a/tools/RAiDER/models/ncmr.py
+++ b/tools/RAiDER/models/ncmr.py
@@ -10,7 +10,7 @@
 
 from pyproj import CRS
 
-from RAiDER.models.weatherModel import WeatherModel
+from RAiDER.models.weatherModel import WeatherModel, TIME_RES
 from RAiDER.logger import logger
 from RAiDER.utilFcns import (
     writeWeatherVars2NETCDF4,
@@ -36,12 +36,12 @@ def __init__(self):
         self._classname = 'ncmr'                     # name of the custom weather model
         self._dataset = 'ncmr'                       # same name as above
         self._Name = 'NCMR'                          # name of the new weather model (in Capital)
+        self._time_res = TIME_RES[self._dataset.upper()]
 
         # Tuple of min/max years where data is available.
         self._valid_range = (datetime.datetime(2015, 12, 1), "Present")
         # Availability lag time in days/hours
         self._lag_time = datetime.timedelta(hours=6)
-        self._time_res = 1
 
         # model constants
         self._k1 = 0.776   # [K/Pa]
@@ -67,7 +67,7 @@ def _fetch(self, out):
         Fetch weather model data from NCMR: note we only extract the lat/lon bounds for this weather model;
         fetching data is not needed here as we don't actually download data , data exist in same system
         '''
-        time = self._time 
+        time = self._time
 
         # Auxillary function:
         '''
diff --git a/tools/RAiDER/models/weatherModel.py b/tools/RAiDER/models/weatherModel.py
index 53abee761..f0c408612 100755
--- a/tools/RAiDER/models/weatherModel.py
+++ b/tools/RAiDER/models/weatherModel.py
@@ -22,6 +22,15 @@
     robmax, robmin, write2NETCDF4core, calcgeoh, transform_coords
 )
 
+TIME_RES = {'GMAO': 3,
+            'ECMWF': 1,
+            'HRES': 6,
+            'HRRR': 1,
+            'WRF': 1,
+            'NCMR': 1
+            }
+
+
 
 class WeatherModel(ABC):
     '''
@@ -122,6 +131,9 @@ def __str__(self):
     def Model(self):
         return self._Name
 
+    def dtime(self):
+        return self._time_res
+
     def fetch(self, out, ll_bounds, time):
         '''
         Checks the input datetime against the valid date range for the model and then
@@ -272,10 +284,10 @@ def _convertmb2Pa(self, pres):
 
     def _get_heights(self, lats, geo_hgt, geo_ht_fill=np.nan):
         '''
-        Transform geo heights to actual heights
+        Transform geo heights to WGS84 ellipsoidal heights
         '''
         geo_ht_fix = np.where(geo_hgt != geo_ht_fill, geo_hgt, np.nan)
-        self._zs = util._geo_to_ht(lats, geo_ht_fix)
+        self._zs = util.geo_to_ht(lats, geo_ht_fix)
 
     def _find_e(self):
         """Check the type of e-calculation needed"""
diff --git a/tools/RAiDER/models/wrf.py b/tools/RAiDER/models/wrf.py
index d35f795bd..827495a54 100644
--- a/tools/RAiDER/models/wrf.py
+++ b/tools/RAiDER/models/wrf.py
@@ -2,7 +2,7 @@
 import scipy.io.netcdf as netcdf
 from pyproj import CRS, Transformer
 
-from RAiDER.models.weatherModel import WeatherModel
+from RAiDER.models.weatherModel import WeatherModel, TIME_RES
 
 
 # Need to incorporate this snippet into this part of the code.
@@ -27,9 +27,11 @@ def __init__(self):
         self._k2 = 0.233  # K/Pa
         self._k3 = 3.75e3  # K^2/Pa
 
+
         # Currently WRF is using RH instead of Q to get E
         self._humidityType = 'rh'
         self._Name = 'WRF'
+        self._time_res = TIME_RES[self._Name]
 
     def _fetch(self):
         pass
diff --git a/tools/RAiDER/utilFcns.py b/tools/RAiDER/utilFcns.py
index 987a9cbfc..e36b22a3e 100755
--- a/tools/RAiDER/utilFcns.py
+++ b/tools/RAiDER/utilFcns.py
@@ -18,8 +18,10 @@
 
 from RAiDER.constants import (
     _g0 as g0,
-    R_EARTH_MAX as Rmax,
-    R_EARTH_MIN as Rmin,
+    _g1 as G1,
+    R_EARTH_MAX_WGS84 as Rmax,
+    R_EARTH_MIN_WGS84 as Rmin,
+    _THRESHOLD_SECONDS,
 )
 from RAiDER.logger import logger
 
@@ -291,7 +293,7 @@ def writeResultsToXarray(dt, xpts, ypts, zpts, crs, wetDelay, hydroDelay, weathe
         ds.x.attrs["standard_name"] = "projection_x_coordinate"
         ds.x.attrs["long_name"] = "x-coordinate in projected coordinate system"
         ds.x.attrs["units"] = "m"
-    
+
     return ds
 
 
@@ -396,30 +398,59 @@ def _get_g_ll(lats):
     '''
     Compute the variation in gravity constant with latitude
     '''
-    # TODO: verify these constants. In particular why is the reference g different from self._g0?
-    return 9.80616 * (1 - 0.002637 * cosd(2 * lats) + 0.0000059 * (cosd(2 * lats))**2)
+    return G1 * (1 - 0.002637 * cosd(2 * lats) + 0.0000059 * (cosd(2 * lats))**2)
 
 
-def _get_Re(lats):
+def get_Re(lats):
     '''
-    Returns: the ellipsoid as a fcn of latitude
+    Returns earth radius as a function of latitude for WGS84
+    
+    Args: 
+        lats    - ndarray of geodetic latitudes in degrees
+    
+    Returns:
+        ndarray of earth radius at each latitude
+    
+    Example: 
+    >>> import numpy as np
+    >>> from RAiDER.utilFcns import get_Re
+    >>> output = get_Re(np.array([0, 30, 45, 60, 90]))
+    >>> output
+     array([6378137., 6372770.5219805, 6367417.56705189, 6362078.07851428, 6356752.])
+    >>> assert output[0] == 6378137 # (Rmax)
+    >>> assert output[-1] == 6356752 # (Rmin)
     '''
-    # TODO: verify constants, add to base class constants?
     return np.sqrt(1 / (((cosd(lats)**2) / Rmax**2) + ((sind(lats)**2) / Rmin**2)))
 
 
-def _geo_to_ht(lats, hts):
-    """Convert geopotential height to altitude."""
-    # Convert geopotential to geometric height. This comes straight from
-    # TRAIN
-    # Map of g with latitude (I'm skeptical of this equation - Ray)
-    g_ll = _get_g_ll(lats)
-    Re = _get_Re(lats)
+def geo_to_ht(lats, hts):
+    """
+    Convert geopotential height to ellipsoidal heights referenced to WGS84.
+    
+    Note that this formula technically computes height above geoid (geometric height)
+    but the geoid is actually a perfect sphere;
+    Thus returned heights are above a reference ellipsoid, which most assume to be 
+    a sphere (e.g., ECMWF - see https://confluence.ecmwf.int/display/CKB/ERA5%3A+compute+pressure+and+geopotential+on+model+levels%2C+geopotential+height+and+geometric+height#ERA5:computepressureandgeopotentialonmodellevels,geopotentialheightandgeometricheight-Geopotentialheight 
+    - "Geometric Height" and also https://confluence.ecmwf.int/display/CKB/ERA5%3A+data+documentation#ERA5:datadocumentation-Earthmodel).
+    However, by calculating the ellipsoid here we directly reference to WGS84. 
+
+    Compare to MetPy: 
+    (https://unidata.github.io/MetPy/latest/api/generated/metpy.calc.geopotential_to_height.html)
+    # h = (geopotential * Re) / (g0 * Re - geopotential)
+    # Assumes a sphere instead of an ellipsoid
+
+    Args: 
+        lats    - latitude of points of interest
+        hts     - geopotential height at points of interest
+    
+    Returns: 
+        ndarray: geometric heights. These are approximate ellipsoidal heights referenced to WGS84
+    """
+    g_ll = _get_g_ll(lats) # gravity function of latitude
+    Re = get_Re(lats) # Earth radius function of latitude
 
     # Calculate Geometric Height, h
     h = (hts * Re) / (g_ll / g0 * Re - hts)
-    # from metpy
-    # return (geopotential * Re) / (g0 * Re - geopotential)
 
     return h
 
@@ -1135,3 +1166,59 @@ def transform_coords(proj1, proj2, x, y):
     """
     transformer = Transformer.from_crs(proj1, proj2, always_xy=True)
     return transformer.transform(x, y)
+
+
+def get_nearest_wmtimes(t0, time_delta):
+    """"
+    Get the nearest two available times to the requested time given a time step
+
+    Args:
+        t0         - user-requested Python datetime
+        time_delta  - time interval of weather model
+
+    Returns:
+        tuple: list of datetimes representing the one or two closest
+        available times to the requested time
+
+    Example:
+    >>> import datetime
+    >>> from RAiDER.utilFcns import get_nearest_wmtimes
+    >>> t0 = datetime.datetime(2020,1,1,11,35,0)
+    >>> get_nearest_wmtimes(t0, 3)
+     (datetime.datetime(2020, 1, 1, 9, 0), datetime.datetime(2020, 1, 1, 12, 0))
+    """
+    # get the closest time available
+    tclose = round_time(t0, roundTo = time_delta * 60 *60)
+
+    # Just calculate both options and take the closest
+    t2_1 = tclose + timedelta(hours=time_delta)
+    t2_2 = tclose - timedelta(hours=time_delta)
+    t2 = [t2_1 if get_dt(t2_1, t0) < get_dt(t2_2, t0) else t2_2][0]
+
+    # If you're within 5 minutes just take the closest time
+    if get_dt(tclose, t0) < _THRESHOLD_SECONDS:
+        return [tclose]
+    else:
+        if t2 > tclose:
+            return [tclose, t2]
+        else:
+            return [t2, tclose]
+
+def get_dt(t1,t2):
+    '''
+    Helper function for getting the absolute difference in seconds between
+    two python datetimes
+
+    Args:
+        t1, t2  - Python datetimes
+
+    Returns:
+        Absolute difference in seconds between the two inputs
+
+    Examples:
+    >>> import datetime
+    >>> from RAiDER.utilFcns import get_dt
+    >>> get_dt(datetime.datetime(2020,1,1,5,0,0), datetime.datetime(2020,1,1,0,0,0))
+     18000.0
+    '''
+    return np.abs((t1 - t2).total_seconds())