Repository: UW-SASWE/RAT (language: python)

Issue:
Landsat 9 script for surface area calculation is missing from TMSOS. 
Please add landsat9 also in run_sarea_for_res() in run_sarea.py 
```
    # Sentinel-2
    log.debug(f"Reservoir: {reservoir_name}; Downloading Sentinel-2 data from {start_date} to {end_date}")
    sarea_s2(reservoir_name, reservoir_polygon, start_date, end_date, os.path.join(datadir, 's2'))
    s2_dfpath = os.path.join(datadir, 's2', reservoir_name+'.csv')

    # Landsat-8
    log.debug(f"Reservoir: {reservoir_name}; Downloading Landsat-8 data from {start_date} to {end_date}")
    sarea_l8(reservoir_name, reservoir_polygon, start_date, end_date, os.path.join(datadir, 'l8'))
    l8_dfpath = os.path.join(datadir, 'l8', reservoir_name+'.csv')

    # Sentinel-1
    log.debug(f"Reservoir: {reservoir_name}; Downloading Sentinel-1 data from {start_date} to {end_date}")
    s1_dfpath = sarea_s1(reservoir_name, reservoir_polygon, start_date, end_date, os.path.join(datadir, 'sar'))
    s1_dfpath = os.path.join(datadir, 'sar', reservoir_name+'_12d_sar.csv')

    tmsos = TMS(reservoir_name, reservoir_area)
    result = tmsos.tms_os(l8_dfpath=l8_dfpath, s2_dfpath=s2_dfpath, s1_dfpath=s1_dfpath)
```

Also add sarea_cli_l9.py script which got missed in some previous commit.


Discussion / hints:

Fix patch:
diff --git a/.gitignore b/.gitignore
index a7feba65..56dee962 100644
--- a/.gitignore
+++ b/.gitignore
@@ -207,6 +207,7 @@ backend/scripts/data_processing/_download_data.py
 
 
 #SM
+testing/*
 data/*
 global_data/*
 trial_notebooks/*
diff --git a/docs/development_guide.md b/docs/development_guide.md
new file mode 100644
index 00000000..c703cae1
--- /dev/null
+++ b/docs/development_guide.md
@@ -0,0 +1,15 @@
+# Development guide
+
+## Getting RAT from github
+
+1. Create a conda environment where RAT will be installed - `conda create -n rat_dev`
+
+2. Clone the `RAT` project into a directory and navigate into it - `git clone https://github.com/UW-SASWE/RAT.git; cd RAT`. This will clone the RAT project in the directory.
+
+3. Install `RAT` requirements into the `rat_dev` environment, optionally using mamba (recommended). `conda activate rat_dev; conda install mamba; mamba env update --file environment.yml`. This will install the requirements of `RAT` in the environment.
+
+4. Now "install" `RAT` in development mode - `conda develop src/`. You should see a message that states that the src was added to path. 
+
+Congratulations, `RAT` is now installed! You should be able to import it using `import rat`.
+
+5. `RAT` uses a Command Line Interface (CLI) to perform various tasks such as initialization and run. The CLI functionality gets automatically activated when it is installed via a conda command (`conda install rat`), but it doesn't automatically get registered when installing in the "development" mode. To use the CLI functionality, you can directly call the python script that handles the CLI and pass along any of the arguments. For instance, to initialize rat (which would be `rat init` in a normal install), you would write `python src/cli/rat_cli.py init`.
\ No newline at end of file
diff --git a/environment.yml b/environment.yml
index 61660015..81041bc6 100644
--- a/environment.yml
+++ b/environment.yml
@@ -212,7 +212,6 @@ dependencies:
   - ripgrep=13.0.0
   - rsa=4.9
   - rtree=1.0.1
-  - ruamel_yaml=0.15.80
   - scikit-learn=1.1.2
   - scipy=1.9.1
   - setuptools=65.5.0
@@ -262,3 +261,4 @@ dependencies:
   - zipp=3.9.0
   - zlib=1.2.13
   - zstd=1.5.2
+  - ruamel.yaml
\ No newline at end of file
diff --git a/src/rat/cli/rat_cli.py b/src/rat/cli/rat_cli.py
index 742d679f..f84281a5 100644
--- a/src/rat/cli/rat_cli.py
+++ b/src/rat/cli/rat_cli.py
@@ -5,7 +5,7 @@
 import shutil
 import requests, zipfile, io
 import yaml
-import ruamel.yaml as ryaml
+import ruamel_yaml as ryaml
 
 
 def init_func(args):
@@ -30,6 +30,10 @@ def init_func(args):
         global_data = args.global_data
     global_data_dir = project_dir.joinpath("global_data")
 
+    if args.secrets is not None:
+        secrets_fp = Path(args.secrets).resolve()
+        assert secrets_fp.exists(), f"Secrets file {secrets_fp} does not exist"
+
     # create additional directories
     data_dir = project_dir.joinpath('data')
     data_dir.mkdir(exist_ok=True)
@@ -37,6 +41,7 @@ def init_func(args):
     models_dir.mkdir(exist_ok=True)
     params_dir = project_dir.joinpath('params')
     params_dir.mkdir(exist_ok=True)
+    rat_config_fp = params_dir.joinpath('rat_config.yaml')
 
     #### Model installation
     # install metsim
@@ -52,7 +57,7 @@ def init_func(args):
     subprocess.run(cmd)
     
     # import download links
-    from rat_init_config import DOWNLOAD_LINKS
+    from rat.cli.rat_init_config import DOWNLOAD_LINKS
 
     # install route
     print(f"Downloading source code of routing model")
@@ -110,19 +115,83 @@ def init_func(args):
     #### update params
     update_param_file(
         project_dir,
-        global_data == 'Y',
-        n_cores
+        config_path=rat_config_fp,
+        global_data_downloaded=(global_data == 'Y'),
+        n_cores=n_cores,
+        secrets=secrets_fp
+    )
+
+def test_func(args):
+    from rat.cli.rat_test_config import DOWNLOAD_LINK, PATHS, PARAMS
+
+    test_basin_options = PARAMS.keys()
+    test_basin = args.test_basin
+    assert test_basin in test_basin_options, f"Please specify the correct test basin. Acceptable values are {test_basin_options})"
+
+    project_dir = Path(args.project_dir).resolve()
+    assert project_dir.exists(), f"{project_dir} does not exist - please pass a valid rat project directory after initialization using the `rat init` command."
+
+    test_param_fp = project_dir / 'params' / 'test_config.yml'
+    
+    secrets_fp = Path(args.secrets).resolve()
+    assert secrets_fp.exists(), f"{secrets_fp} does not exist - rat requires secrets.ini file to be passed. Please refer to documentation for more details."
+
+    # Download test data
+    data_dir = project_dir / 'data'
+    assert data_dir.exists(), f"{data_dir} does not exist - rat has not been initialized properly."
+
+    test_data_link = DOWNLOAD_LINK['test_data']
+    r = requests.get(test_data_link)
+    z = zipfile.ZipFile(io.BytesIO(r.content))
+    z.extractall(data_dir)
+
+    n_cores = 4
+    try:
+        from multiprocessing import cpu_count
+        n_cores = cpu_count()
+    except Exception as e:
+        print(f"Failed to determine number of cores: {e}")
+
+    update_param_file(
+        project_dir=project_dir, 
+        config_path=test_param_fp,
+        global_data_downloaded=False, 
+        n_cores=n_cores,
+        secrets=secrets_fp,
+        any_other_suffixes=PATHS[test_basin],
+        any_other_args=PARAMS[test_basin]
     )
 
+    run_args = argparse.Namespace()
+    run_args.param = test_param_fp
+    run_args.operational_latency = None
+    run_func(run_args)
+
 
 def update_param_file(
         project_dir: Path,
-        global_data_downloaded: bool = False, # if global data was downloaded, then we can use the global_data dir
-        n_cores: int = None
+        config_path: Path = None,
+        global_data_downloaded: bool = False,   # if global data was downloaded, then we can use the global_data dir
+        n_cores: int = None,
+        secrets: Path = None,
+        any_other_suffixes: dict = None,        # any other suffixes that need to be added to the config file
+        any_other_args: dict = None             # any other args that need to be added to the config file
     ):
+    """Creates RAT config file with project specific information.
+
+    Args:
+        project_dir (Path): Directory of the rat project.
+        config_path (Path, optional): Where to save the configuration path? Defaults to `project_dir/params/rat_config.yml`.
+        global_data_downloaded (bool, optional): Do global datasets need to be downloaded to initialize rat? . Defaults to False.
+        n_cores (int, optional): Number of cores to use. Defaults to None.
+        secrets (Path, optional): Path of the `secrets.ini` file. Please refer to documentation for more details. Defaults to None.
+        any_other_suffixes (dict, optional): _description_. Defaults to None.
+    """
 
     config_template_path = (project_dir / 'params' / 'rat_config_template.yml')
-    config_path = (project_dir / 'params' / 'rat_config.yml')
+
+    if config_path is None:
+        config_path = (project_dir / 'params' / 'rat_config.yml')
 
     ryaml_client = ryaml.YAML()
     config_template = ryaml_client.load(config_template_path.read_text())
@@ -140,16 +209,29 @@ def update_param_file(
             for k2, v2 in v1.items():
                 config_template[k1][k2] = str(project_dir.joinpath(v2))
 
+    if any_other_args is not None:
+        for k1, v1 in any_other_args.items():
+            for k2, v2 in v1.items():
+                config_template[k1][k2] = v2
+    
+    if any_other_suffixes is not None:
+        for k1, v1 in any_other_suffixes.items():
+            for k2, v2 in v1.items():
+                config_template[k1][k2] = str(project_dir.joinpath(v2))
+
     # update the number of cores
     config_template['GLOBAL']['multiprocessing'] = n_cores
 
+    # if secrets were provided, update the config file
+    if secrets is not None:
+        config_template['CONFIDENTIAL']['secrets'] = str(secrets)
+
     ryaml_client.dump(config_template, config_path.open('w'))
 
 
 def run_func(args):
     from rat.run_rat import run_rat
-
-    run_rat(args.param)
+    run_rat(args.param, args.operational_latency)
 
 def main():
     ## CLI interface
@@ -175,7 +257,14 @@ def main():
         action='store_true',
         dest='global_data',
         required=False,
-        default=None
+        default=False
+    )
+    init_parser.add_argument(
+        '-s', '--secrets', 
+        help='Specify the path of secrets.ini file', 
+        action='store',
+        dest='secrets',
+        required=False
     )
 
     init_parser.set_defaults(func=init_func)
@@ -190,8 +279,44 @@ def main():
         dest='param',
         required=True
     )
+    run_parser.add_argument(
+        '-o', '--operational',
+        help='RAT Operational Latency in days',
+        action='store',
+        dest='operational_latency',
+        required=False
+    )
     
     run_parser.set_defaults(func=run_func)
+    
+    # Test command
+    test_parser = command_parsers.add_parser('test', help='Test RAT')
+
+    test_parser.add_argument(
+        '-b', '--basin', 
+        help='Specify name of test basin', 
+        action='store',
+        dest='test_basin',
+        required=True
+    )
+
+    test_parser.add_argument(
+        '-d', '--dir', 
+        help='Specify RAT project directory (where rat init was run)', 
+        action='store',
+        dest='project_dir',
+        required=True
+    )
+    
+    test_parser.add_argument(
+        '-s', '--secrets', 
+        help='Specify the path of secrets.ini file', 
+        action='store',
+        dest='secrets',
+        required=True
+    )
+    
+    test_parser.set_defaults(func=test_func)
 
     args = p.parse_args()
     args.func(args)
diff --git a/src/rat/cli/rat_init_config.py b/src/rat/cli/rat_init_config.py
index 3aa6c02e..7e2f89ee 100644
--- a/src/rat/cli/rat_init_config.py
+++ b/src/rat/cli/rat_init_config.py
@@ -1,8 +1,8 @@
 DOWNLOAD_LINKS = {
     'route_model': "https://www.dropbox.com/s/9jwep2g5pyj8sni/routing.zip?dl=1",
-    'params': "https://www.dropbox.com/s/9n6b1pn73fgl0y9/params.zip?dl=1",
+    'params': "https://www.dropbox.com/s/obb5bd24eb8u12h/params.zip?dl=1",
     'global_data': "https://www.dropbox.com/s/u8vc3oxujmaak97/global_data.zip?dl=1",
-    'global_vic_params': "https://www.dropbox.com/s/jsg2wu62qi2ltwz/global_vic_params.zip?dl=1"
+    'global_vic_params': "https://www.dropbox.com/s/jsg2wu62qi2ltwz/global_vic_params.zip?dl=1",
 }
 
 SUFFIXES_GLOBAL = {
diff --git a/src/rat/core/run_postprocessing.py b/src/rat/core/run_postprocessing.py
index 45b7a121..58d3a2bb 100644
--- a/src/rat/core/run_postprocessing.py
+++ b/src/rat/core/run_postprocessing.py
@@ -4,6 +4,7 @@
 import xarray as xr
 import geopandas as gpd
 import warnings
+import datetime
 
 warnings.filterwarnings("ignore")
 
@@ -125,11 +126,20 @@ def calc_E(res_data, start_date, end_date, forcings_path, vic_res_path, sarea_pa
 
 
 def calc_outflow(inflowpath, dspath, epath, area, savepath):
-    inflow = pd.read_csv(inflowpath, parse_dates=["date"])
-    E = pd.read_csv(epath, parse_dates=['time'])
+    if os.path.isfile(inflowpath):
+        inflow = pd.read_csv(inflowpath, parse_dates=["date"])
+    else:
+        print('Inflow file does not exist and Outflow cannot be calculated.')
+    if os.path.isfile(epath):
+        E = pd.read_csv(epath, parse_dates=['time'])
+    else:
+        print('Evaporation file does not exist and Outflow cannot be calculated.')
 
     if isinstance(dspath, str):
-        df = pd.read_csv(dspath, parse_dates=['date'])
+        if os.path.isfile(dspath):
+            df = pd.read_csv(dspath, parse_dates=['date'])
+        else:
+            print('Storage Change file does not exist and Outflow cannot be calculated.')
     else:
         df = dspath
     
@@ -159,12 +169,16 @@ def calc_outflow(inflowpath, dspath, epath, area, savepath):
     df.to_csv(savepath, index=False)
 
 
-def run_postprocessing(basin_name, basin_data_dir, reservoir_shpfile, reservoir_shpfile_column_dict, aec_dir_path, start_date, end_date,
+def run_postprocessing(basin_name, basin_data_dir, reservoir_shpfile, reservoir_shpfile_column_dict, aec_dir_path, start_date, end_date, rout_init_state_save_file, use_rout_state,
                             evap_datadir, dels_savedir, outflow_savedir, vic_status, routing_status, gee_status):
     # read file defining mapped resrvoirs
     # reservoirs_fn = os.path.join(project_dir, 'backend/data/ancillary/RAT-Reservoirs.geojson')
     reservoirs = gpd.read_file(reservoir_shpfile)
     start_date_str = start_date.strftime("%Y-%m-%d")
+    if(use_rout_state):
+        start_date_str_evap = (start_date-datetime.timedelta(days=100)).strftime("%Y-%m-%d")
+    else:
+        start_date_str_evap = start_date_str
     end_date_str = end_date.strftime("%Y-%m-%d")
     # Defining flags
     EVAP_STATUS = 0
@@ -180,17 +194,20 @@ def run_postprocessing(basin_name, basin_data_dir, reservoir_shpfile, reservoir_
         aec_dir = aec_dir_path
 
         for reservoir_no,reservoir in reservoirs.iterrows():
-            # Reading reservoir information
-            reservoir_name = str(reservoir[reservoir_shpfile_column_dict['unique_identifier']])
-            sarea_path = os.path.join(sarea_raw_dir, reservoir_name + ".csv")
-            savepath = os.path.join(dels_savedir, reservoir_name + ".csv")
-            aecpath = os.path.join(aec_dir, reservoir_name + ".csv")
-
-            if os.path.isfile(sarea_path):
-                log.debug(f"Calculating ∆S for {reservoir_name}, saving at: {savepath}")
-                calc_dels(aecpath, sarea_path, savepath)
-            else:
-                log.debug(f"{sarea_path} not found; skipping ∆S calculation")
+            try:
+                # Reading reservoir information
+                reservoir_name = str(reservoir[reservoir_shpfile_column_dict['unique_identifier']])
+                sarea_path = os.path.join(sarea_raw_dir, reservoir_name + ".csv")
+                savepath = os.path.join(dels_savedir, reservoir_name + ".csv")
+                aecpath = os.path.join(aec_dir, reservoir_name + ".csv")
+
+                if os.path.isfile(sarea_path):
+                    log.debug(f"Calculating ∆S for {reservoir_name}, saving at: {savepath}")
+                    calc_dels(aecpath, sarea_path, savepath)
+                else:
+                    log.debug(f"{sarea_path} not found; skipping ∆S calculation")
+            except:
+                log.exception(f"∆S for {reservoir_name} could not be calculated.")
         DELS_STATUS=1
     else:
         log.debug("Cannot Calculate ∆S because GEE Run Failed.")
@@ -199,21 +216,27 @@ def run_postprocessing(basin_name, basin_data_dir, reservoir_shpfile, reservoir_
     # Evaporation
     if(vic_status and gee_status):
         log.debug("Retrieving Evaporation")
-        vic_results_path = os.path.join(basin_data_dir,'vic', "vic_outputs/nc_fluxes."+start_date_str+".nc")
+        if(use_rout_state):
+            vic_results_path = rout_init_state_save_file
+        else:
+            vic_results_path = os.path.join(basin_data_dir,'vic', "vic_outputs/nc_fluxes."+start_date_str+".nc")
         forcings_path = os.path.join(basin_data_dir,'vic', "vic_inputs/*.nc")
 
         for reservoir_no,reservoir in reservoirs.iterrows():
-            # Reading reservoir information
-            reservoir_name = str(reservoir[reservoir_shpfile_column_dict['unique_identifier']])
-            sarea_path = os.path.join(sarea_raw_dir, reservoir_name + ".csv")
-            e_path = os.path.join(evap_datadir, reservoir_name + ".csv")
-            
-            if os.path.isfile(sarea_path):
-                log.debug(f"Calculating Evaporation for {reservoir_name}")
-                # calc_E(e_path, respath, vic_results_path)
-                calc_E(reservoir, start_date_str, end_date_str, forcings_path, vic_results_path, sarea_path, e_path)
-            else:
-                log.debug(f"{sarea_path} not found; skipping evaporation calculation")
+            try:
+                # Reading reservoir information
+                reservoir_name = str(reservoir[reservoir_shpfile_column_dict['unique_identifier']])
+                sarea_path = os.path.join(sarea_raw_dir, reservoir_name + ".csv")
+                e_path = os.path.join(evap_datadir, reservoir_name + ".csv")
+                
+                if os.path.isfile(sarea_path):
+                    log.debug(f"Calculating Evaporation for {reservoir_name}")
+                    # calc_E(e_path, respath, vic_results_path)
+                    calc_E(reservoir, start_date_str_evap, end_date_str, forcings_path, vic_results_path, sarea_path, e_path)
+                else:
+                    log.debug(f"{sarea_path} not found; skipping evaporation calculation")
+            except:
+                log.exception(f"Evaporation for {reservoir_name} could not be calculated.")
         EVAP_STATUS = 1
     elif((not vic_status) and (not gee_status)):
         log.debug("Cannot Retrieve Evaporation because both VIC and GEE Run Failed.")
@@ -228,16 +251,19 @@ def run_postprocessing(basin_name, basin_data_dir, reservoir_shpfile, reservoir_
         inflow_dir = os.path.join(basin_data_dir,'ro', "rout_inflow")
 
         for reservoir_no,reservoir in reservoirs.iterrows():
-            # Reading reservoir information
-            reservoir_name = str(reservoir[reservoir_shpfile_column_dict['unique_identifier']])
-            deltaS = os.path.join(dels_savedir, reservoir_name + ".csv")
-            inflowpath = os.path.join(inflow_dir, reservoir_name[:5] + ".csv")  ## Routing model keeps only first 5 letters of the reservoir name as file name
-            epath = os.path.join(evap_datadir, reservoir_name + ".csv")
-            a = float(reservoir[reservoir_shpfile_column_dict['area_column']])
-
-            savepath = os.path.join(outflow_savedir, reservoir_name + ".csv")
-            log.debug(f"Calculating Outflow for {reservoir_name} saving at: {savepath}")
-            calc_outflow(inflowpath, deltaS, epath, a, savepath)
+            try:
+                # Reading reservoir information
+                reservoir_name = str(reservoir[reservoir_shpfile_column_dict['unique_identifier']])
+                deltaS = os.path.join(dels_savedir, reservoir_name + ".csv")
+                inflowpath = os.path.join(inflow_dir, reservoir_name[:5] + ".csv")  ## Routing model keeps only first 5 letters of the reservoir name as file name
+                epath = os.path.join(evap_datadir, reservoir_name + ".csv")
+                a = float(reservoir[reservoir_shpfile_column_dict['area_column']])
+
+                savepath = os.path.join(outflow_savedir, reservoir_name + ".csv")
+                log.debug(f"Calculating Outflow for {reservoir_name} saving at: {savepath}")
+                calc_outflow(inflowpath, deltaS, epath, a, savepath)
+            except:
+                log.exception(f"Outflow for {reservoir_name} could not be calculated")
         OUTFLOW_STATUS = 1
     else:
         log.debug("Cannot Calculate Outflow because either evaporation, ∆S or Inflow is missing.")
diff --git a/src/rat/core/run_sarea.py b/src/rat/core/run_sarea.py
index 1b4f7b4e..ab5d357a 100644
--- a/src/rat/core/run_sarea.py
+++ b/src/rat/core/run_sarea.py
@@ -6,6 +6,7 @@
 
 from rat.core.sarea.sarea_cli_s2 import sarea_s2
 from rat.core.sarea.sarea_cli_l8 import sarea_l8
+from rat.core.sarea.sarea_cli_l9 import sarea_l9
 from rat.core.sarea.sarea_cli_sar import sarea_s1
 from rat.core.sarea.TMS import TMS
 
@@ -35,6 +36,11 @@ def run_sarea_for_res(reservoir_name, reservoir_area, reservoir_polygon, start_d
     log.debug(f"Reservoir: {reservoir_name}; Downloading Landsat-8 data from {start_date} to {end_date}")
     sarea_l8(reservoir_name, reservoir_polygon, start_date, end_date, os.path.join(datadir, 'l8'))
     l8_dfpath = os.path.join(datadir, 'l8', reservoir_name+'.csv')
+    
+    # Landsat-9
+    log.debug(f"Reservoir: {reservoir_name}; Downloading Landsat-9 data from {start_date} to {end_date}")
+    sarea_l9(reservoir_name, reservoir_polygon, start_date, end_date, os.path.join(datadir, 'l9'))
+    l9_dfpath = os.path.join(datadir, 'l9', reservoir_name+'.csv')
 
     # Sentinel-1
     log.debug(f"Reservoir: {reservoir_name}; Downloading Sentinel-1 data from {start_date} to {end_date}")
@@ -42,7 +48,7 @@ def run_sarea_for_res(reservoir_name, reservoir_area, reservoir_polygon, start_d
     s1_dfpath = os.path.join(datadir, 'sar', reservoir_name+'_12d_sar.csv')
 
     tmsos = TMS(reservoir_name, reservoir_area)
-    result = tmsos.tms_os(l8_dfpath=l8_dfpath, s2_dfpath=s2_dfpath, s1_dfpath=s1_dfpath)
+    result = tmsos.tms_os(l9_dfpath=l9_dfpath, l8_dfpath=l8_dfpath, s2_dfpath=s2_dfpath, s1_dfpath=s1_dfpath)
 
     tmsos_savepath = os.path.join(datadir, reservoir_name+'.csv')
     log.debug(f"Saving surface area of {reservoir_name} at {tmsos_savepath}")
diff --git a/src/rat/core/run_vic.py b/src/rat/core/run_vic.py
index abf471e6..0f0dc3d0 100644
--- a/src/rat/core/run_vic.py
+++ b/src/rat/core/run_vic.py
@@ -9,6 +9,8 @@
 import time
 import math
 import datetime
+from pathlib import Path
+from functools import partial
 
 from logging import getLogger
 from rat.utils.logging import LOG_NAME, NOTIFICATION
@@ -38,29 +40,37 @@ def run_vic(self, np=16, cd=None):
 
         ret_code = run_command(arg, shell=True, executable="/bin/bash")
 
-    def generate_routing_input_state(self, ndays, rout_input_state_file, save_path):
-        if os.path.isfile(rout_input_state_file):
+    def generate_routing_input_state(self, ndays, rout_input_state_file, save_path, use_rout_state):
+        if os.path.isfile(rout_input_state_file) and (use_rout_state):
             print('Routing input state fle exists at '+rout_input_state_file)
-            prev_vic_output = xr.open_mfdataset(rout_input_state_file).load()
-            last_existing_time = prev_vic_output.time[-1]
+            new_vic_output = xr.open_mfdataset(self.vic_result).load()
+            first_existing_time = new_vic_output.time[0]
+            new_vic_output.close()
+
+            #Preprocessing function for merging netcdf files
+            def _remove_coinciding_days(ds, cutoff_time, ndays):
+                file_name = ds.encoding["source"]
+                file_stem = Path(file_name).stem
+                if('ro_init' in file_stem):
+                    return ds.sel(time=slice(cutoff_time - np.timedelta64(ndays,'D') , cutoff_time - np.timedelta64(1,'D')))
+                else:
+                    return ds
+            remove_coinciding_days_func = partial(_remove_coinciding_days, cutoff_time=first_existing_time, ndays=ndays)
             
+            # Merging previous and new vic outputs
             try:
-                prev_new_combine_vic_output = xr.open_mfdataset([rout_input_state_file,self.vic_result],{'time':365})
-                save_vic_output = prev_new_combine_vic_output.sel(time=slice(last_existing_time - np.timedelta64(ndays,'D') , None))
+                save_vic_output = xr.open_mfdataset([rout_input_state_file,self.vic_result],{'time':365}, preprocess=remove_coinciding_days_func)
+                save_vic_output.to_netcdf(save_path)
+                print('Latest Routing input state fle saved at '+save_path)
             except:
                 ## In case routing state file has dates matching with the vic output file. 
                 print('Rout input state has same dates as vic_output_dates')
-                starting_date_routing_file = prev_vic_output.time[0].values.astype('datetime64[us]').astype(datetime.datetime)
-                prev_vic_output.close()
-                return starting_date_routing_file
-            prev_vic_output.close()
         else:
             new_vic_output = xr.open_mfdataset(self.vic_result,{'time':365})
             save_vic_output = new_vic_output[dict(time=slice(-ndays,None))]
             new_vic_output.close()
-        vic_output_start_date = save_vic_output.time[0].values.astype('datetime64[us]').astype(datetime.datetime)
-        save_vic_output.to_netcdf(save_path)
-        return vic_output_start_date
+            save_vic_output.to_netcdf(save_path)
+        
 
     def disagg_results(self, rout_input_state_file):
         log.log(NOTIFICATION, "Started disaggregating VIC results")
diff --git a/src/rat/core/sarea/TMS.py b/src/rat/core/sarea/TMS.py
index 031c2bdc..fdc658c5 100644
--- a/src/rat/core/sarea/TMS.py
+++ b/src/rat/core/sarea/TMS.py
@@ -15,7 +15,7 @@ def __init__(self, reservoir_name, area=None, AREA_DEVIATION_THRESHOLD_PCNT=5):
         Args:
             reservoir_name (_type_): _description_
             area (_type_, optional): _description_
-            AREA_DEVIATION_THRESHOLD_PCNT (float, optional): _description_. Defaults to 5.
+            AREA_DEVIATION_THRESHOLD_PCNT (float, optional): _description_. Defaults to 25% for area<10 sq. km,  10% for area<100 sq. km, and 5% otherwise.
         Raises:
             Exception: _description_
         Returns:
@@ -23,6 +23,10 @@ def __init__(self, reservoir_name, area=None, AREA_DEVIATION_THRESHOLD_PCNT=5):
         """
         self.reservoir_name = reservoir_name
         self.area = area
+        if self.area < 100:
+            AREA_DEVIATION_THRESHOLD_PCNT=10
+        elif self.area < 10:
+            AREA_DEVIATION_THRESHOLD_PCNT=25
 
         self.AREA_DEVIATION_THRESHOLD = self.area * AREA_DEVIATION_THRESHOLD_PCNT/100
 
@@ -311,9 +315,10 @@ def deviation_correction(area, DEVIATION_THRESHOLD, AREA_COL_NAME='area'):
 
     inner_area.loc[:, 'corrected_trend'] = inner_area['trend']
     inner_area.loc[inner_area['erroneous'], 'corrected_trend'] = inner_area['sar_trend']
-
-    areas = backcalculate(inner_area[AREA_COL_NAME], inner_area['corrected_trend'], inner_area['erroneous'])
-    inner_area[AREA_COL_NAME] = areas
+    print(inner_area)
+    if(not inner_area['erroneous'].empty):
+        areas = backcalculate(inner_area[AREA_COL_NAME], inner_area['corrected_trend'], inner_area['erroneous'])
+        inner_area[AREA_COL_NAME] = areas
 
     return inner_area
 
@@ -364,7 +369,7 @@ def trend_based_correction(area, sar, AREA_DEVIATION_THRESHOLD=25, TREND_DEVIATI
 
     area_filtered.loc[:, 'days_passed'] = area_filtered.index.to_series().diff().dt.days
     area_filtered.loc[:, 'trend'] = area_filtered['filtered_area'].diff()/area_filtered['days_passed']
-
+    
     sar, area_filtered = clip_ts(sar, area_filtered, which='left')
     # sometimes the sar time-series has duplicate values which have to be removed
     sar = sar[~sar.index.duplicated(keep='first')]   # https://stackoverflow.com/a/34297689/4091712
diff --git a/src/rat/core/sarea/sarea_cli_l9.py b/src/rat/core/sarea/sarea_cli_l9.py
new file mode 100644
index 00000000..787d25d9
--- /dev/null
+++ b/src/rat/core/sarea/sarea_cli_l9.py
@@ -0,0 +1,492 @@
+# import geemap as gee
+import ee
+from datetime import datetime, timedelta, timezone
+import pandas as pd
+import time
+import os
+from random import randint
+import argparse
+from itertools import zip_longest
+from rat.ee_utils.ee_utils import poly2feature
+
+from rat.ee_utils.ee_utils import poly2feature
+from rat.utils.logging import LOG_NAME, NOTIFICATION
+from rat.utils.utils import days_between
+from logging import getLogger
+
+log = getLogger(f"{LOG_NAME}.{__name__}")
+
+def grouper(iterable, n, *, incomplete='fill', fillvalue=None):
+    "Collect data into non-overlapping fixed-length chunks or blocks"
+    # grouper('ABCDEFG', 3, fillvalue='x') --> ABC DEF Gxx
+    # grouper('ABCDEFG', 3, incomplete='strict') --> ABC DEF ValueError
+    # grouper('ABCDEFG', 3, incomplete='ignore') --> ABC DEF
+    args = [iter(iterable)] * n
+    if incomplete == 'fill':
+        return zip_longest(*args, fillvalue=fillvalue)
+    if incomplete == 'strict':
+        return zip(*args, strict=True)
+    if incomplete == 'ignore':
+        return zip(*args)
+    else:
+        raise ValueError('Expected fill, strict, or ignore')
+
+
+# l9 = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")
+# gswd = ee.Image("JRC/GSW1_3/GlobalSurfaceWater").select("occurrence")
+# Threshold = 0.0
+
+# NEW STUFF
+l9 = ee.ImageCollection("LANDSAT/LC09/C02/T1_L2")
+gswd = ee.Image("JRC/GSW1_3/GlobalSurfaceWater")
+rgb_vis_params = {"bands":["B4","B3","B2"],"min":0,"max":0.4}
+
+NDWI_THRESHOLD = 0.3
+SMALL_SCALE = 30
+MEDIUM_SCALE = 120
+LARGE_SCALE = 500
+BUFFER_DIST = 500
+CLOUD_COVER_LIMIT = 90
+start_date = ee.Date('2019-01-01')
+end_date = ee.Date('2019-02-01')
+TEMPORAL_RESOLUTION = 16
+RESULTS_PER_ITER = 5
+
+
+# s2_subset = s2.filterBounds(aoi).filterDate(start_date, end_date)
+
+##############################################
+##       Defining necessary functions       ##
+##############################################
+
+
+##########################
+## Dealing with Cloud   ##
+##########################
+
+def preprocess(im):
+    clipped = im   # clipping adds processing overhead, setting clipped = im
+    
+    # clipped = ee.Image(ee.Algorithms.If('BQA' in clipped.bandNames(), preprocess_image_mask(im), clipped.updateMask(ee.Image.constant(1))))
+    # Mask appropriate QA bits
+    QA = im.select(['QA_PIXEL'])
+    
+    cloudShadowBitMask = 1 << 3
+    cloudsBitMask = 1 << 4
+    
+    cloud = (QA.bitwiseAnd(cloudsBitMask).neq(0)).Or(QA.bitwiseAnd(cloudShadowBitMask).neq(0)).rename("cloud")
+    
+    clipped = clipped.updateMask(cloud.eq(0).select("cloud"))
+    
+    # SR scaling
+    clipped = clipped.select('SR_B.').multiply(0.0000275).add(-0.2)
+    clipped = clipped.addBands(cloud)
+    
+    clipped = clipped.set('system:time_start', im.get('system:time_start'))
+    clipped = clipped.set('system:time_end', im.get('system:time_end'))
+    # clipped = clipped.set('cloud_area', cloud_area)
+    
+    return clipped
+
+
+##########################################################/
+## Processing individual images - water classification   ##
+##########################################################/
+
+def identify_water_cluster(im):
+    im = ee.Image(im)
+    mbwi = im.select('MBWI')
+
+    max_cluster_value = ee.Number(im.select('cluster').reduceRegion(
+        reducer = ee.Reducer.max(),
+        geometry = aoi,
+        scale = LARGE_SCALE,
+        maxPixels =  1e10
+    ).get('cluster'))
+
+    clusters = ee.List.sequence(0, max_cluster_value)
+
+    def calc_avg_mbwi(cluster_val):
+        cluster_val = ee.Number(cluster_val)
+        avg_mbwi = mbwi.updateMask(im.select('cluster').eq(ee.Image(cluster_val))).reduceRegion(
+            reducer = ee.Reducer.mean(),
+            scale = MEDIUM_SCALE,
+            geometry = aoi,
+            maxPixels = 1e10
+        ).get('MBWI')
+        return avg_mbwi
+
+    avg_mbwis = ee.Array(clusters.map(calc_avg_mbwi))
+
+    max_mbwi_index = avg_mbwis.argmax().get(0)
+
+    water_cluster = clusters.get(max_mbwi_index)
+
+    return water_cluster
+
+
+def cordeiro(im):
+    band_subset = ee.List(['NDWI', 'MNDWI', 'SR_B7'])   # using NDWI, MNDWI and B7 (SWIR)
+    sampled_pts = im.select(band_subset).sample(
+        region = aoi,
+        scale = SMALL_SCALE,
+        numPixels = 5e3-1  ## limit of 5k points
+    )
+    
+    ## Agglomerative Clustering isn't available, using Cascade K-Means Clustering based on
+    ##  calinski harabasz's work
+    ## https:##developers.google.com/earth-engine/apidocs/ee-clusterer-wekacascadekmeans
+    clusterer = ee.Clusterer.wekaCascadeKMeans(
+        minClusters = 2,
+        maxClusters = 7,
+        init = True
+    ).train(sampled_pts)
+    
+    classified = im.select(band_subset).cluster(clusterer)
+    im = im.addBands(classified)
+    
+    water_cluster = identify_water_cluster(im)
+    water_map = classified.select('cluster').eq(ee.Image.constant(water_cluster)).select(['cluster'], ['water_map_cordeiro'])
+    im = im.addBands(water_map)
+
+    return im
+
+
+def process_image(im):
+    ## FOR SENTINEL-2
+    # im = im.select(['B3', 'B12', 'B2', 'B4', 'B8', 'B11', 'cloud'])
+    ndwi = im.normalizedDifference(['SR_B5', 'SR_B6']).rename('NDWI')
+    im = im.addBands(ndwi)
+    mndwi = im.normalizedDifference(['SR_B3', 'SR_B6']).rename('MNDWI')
+    im = im.addBands(mndwi)
+    mbwi = im.expression("MBWI = 3*SR_B3-SR_B4-SR_B5-SR_B6-SR_B7", {
+        'SR_B3': im.select('SR_B3'),
+        'SR_B4': im.select('SR_B4'),
+        'SR_B5': im.select('SR_B5'),
+        'SR_B6': im.select('SR_B6'),
+        'SR_B7': im.select('SR_B7')
+    })
+    im = im.addBands(mbwi)
+
+    cloud_area = aoi.area().subtract(im.select('cloud').Not().multiply(ee.Image.pixelArea()).reduceRegion(
+        reducer = ee.Reducer.sum(),
+        geometry = aoi,
+        scale = SMALL_SCALE,
+        maxPixels = 1e10
+    ).get('cloud'))
+    cloud_percent = cloud_area.multiply(100).divide(aoi.area())
+    
+    cordeiro_will_run_when = cloud_percent.lt(CLOUD_COVER_LIMIT)
+
+    # Clusting based
+    im = im.addBands(ee.Image(ee.Algorithms.If(cordeiro_will_run_when, cordeiro(im), ee.Image.constant(-1e6))))  # run cordeiro only if cloud percent is < 90%
+
+    water_area_cordeiro = ee.Number(ee.Algorithms.If(cordeiro_will_run_when,
+        ee.Number(im.select('water_map_cordeiro').eq(1).multiply(ee.Image.pixelArea()).reduceRegion(
+            reducer = ee.Reducer.sum(), 
+            geometry = aoi, 
+            scale = SMALL_SCALE, 
+            maxPixels = 1e10
+        ).get('water_map_cordeiro')),
+        ee.Number(-1e6)
+    ))
+    non_water_area_cordeiro = ee.Number(ee.Algorithms.If(cordeiro_will_run_when,
+        ee.Number(im.select('water_map_cordeiro').neq(1).multiply(ee.Image.pixelArea()).reduceRegion(
+            reducer = ee.Reducer.sum(), 
+            geometry = aoi, 
+            scale = SMALL_SCALE, 
+            maxPixels = 1e10
+        ).get('water_map_cordeiro')),
+        ee.Number(-1e6)
+    ))
+
+    # NDWI based
+    im = im.addBands(ndwi.gte(NDWI_THRESHOLD).select(['NDWI'], ['water_map_NDWI']))
+    water_area_NDWI = ee.Number(im.select('water_map_NDWI').eq(1).multiply(ee.Image.pixelArea()).reduceRegion(
+        reducer = ee.Reducer.sum(), 
+        geometry = aoi, 
+        scale = SMALL_SCALE, 
+        maxPixels = 1e10
+    ).get('water_map_NDWI'))
+    non_water_area_NDWI = ee.Number(im.select('water_map_NDWI').neq(1).multiply(ee.Image.pixelArea()).reduceRegion(
+        reducer = ee.Reducer.sum(), 
+        geometry = aoi, 
+        scale = SMALL_SCALE, 
+        maxPixels = 1e10
+    ).get('water_map_NDWI'))
+    
+    im = im.set('cloud_area', cloud_area.multiply(1e-6))
+    im = im.set('cloud_percent', cloud_percent)
+    im = im.set('water_area_cordeiro', water_area_cordeiro.multiply(1e-6))
+    im = im.set('non_water_area_cordeiro', non_water_area_cordeiro.multiply(1e-6))
+    im = im.set('water_area_NDWI', water_area_NDWI.multiply(1e-6))
+    im = im.set('non_water_area_NDWI', non_water_area_NDWI.multiply(1e-6))
+    
+    return im
+
+
+def postprocess_wrapper(im, bandName, raw_area):
+    im = ee.Image(im)
+    bandName = ee.String(bandName)
+    date = im.get('to_date')
+
+    def postprocess():
+        gswd_masked = gswd.updateMask(im.select(bandName).eq(1))
+        
+        hist = ee.List(gswd_masked.reduceRegion(
+            reducer = ee.Reducer.autoHistogram(minBucketWidth = 1),
+            geometry = aoi,
+            scale = SMALL_SCALE,
+            maxPixels = 1e10
+        ).get('occurrence'))
+        
+        counts = ee.Array(hist).transpose().toList()
+        
+        omega = ee.Number(0.17)
+        count_thresh = ee.Number(counts.map(lambda lis: ee.List(lis).reduce(ee.Reducer.mean())).get(1)).multiply(omega)
+        
+        count_thresh_index = ee.Array(counts.get(1)).gt(count_thresh).toList().indexOf(1)
+        occurrence_thresh = ee.Number(ee.List(counts.get(0)).get(count_thresh_index))
+
+        water_map = im.select([bandName], ['water_map'])
+        gswd_improvement = gswd.clip(aoi).gte(occurrence_thresh).updateMask(water_map.mask().Not()).select(["occurrence"], ["water_map"])
+        
+        improved = ee.ImageCollection([water_map, gswd_improvement]).mosaic().select(['water_map'], ['water_map_zhao_gao'])
+        
+        corrected_area = ee.Number(improved.select('water_map_zhao_gao').multiply(ee.Image.pixelArea()).reduceRegion(
+            reducer = ee.Reducer.sum(), 
+            geometry = aoi, 
+            scale = SMALL_SCALE, 
+            maxPixels = 1e10
+        ).get('water_map_zhao_gao'))
+        
+        improved = improved.set("corrected_area", corrected_area.multiply(1e-6))
+        return improved
+
+    def dont_post_process():
+        improved = ee.Image.constant(-1)
+        improved = improved.set("corrected_area", -1)
+        return improved
+
+    condition = ee.Number(im.get('cloud_percent')).lt(CLOUD_COVER_LIMIT).And(ee.Number(raw_area).gt(0))
+    improved = ee.Image(ee.Algorithms.If(condition, postprocess(), dont_post_process()))
+    
+    improved = improved.set("to_date", date)
+    
+    return improved
+
+
+############################################################/
+## Code from here takes care of the time-series generation   ##
+############################################################/
+def calc_ndwi(im):
+    return im.addBands(im.normalizedDifference(['SR_B5', 'SR_B6']).rename('NDWI'))
+
+def process_date(date):
+    date = ee.Date(date)
+    from_date = date
+    to_date = date.advance(TEMPORAL_RESOLUTION - 1, 'day')
+    l9_subset = l9.filterDate(from_date, to_date).filterBounds(aoi).map(preprocess)
+
+    # im = ee.Image(ee.Algorithms.If(s2_subset.size().neq(0), s2_subset.map(process_image).qualityMosaic('NDWI'), ee.Image.constant(0)))
+    im = ee.Image(ee.Algorithms.If(l9_subset.size().neq(0), l9_subset.map(calc_ndwi).qualityMosaic('NDWI'), ee.Image.constant(0)))
+    
+    im = ee.Image(ee.Algorithms.If(l9_subset.size().neq(0), process_image(im), ee.Image.constant(0)))
+
+    ## im = ee.Algorithms.If(im.bandNames process_image(im))
+    im = im.set('from_date', from_date.format("YYYY-MM-dd"))
+    im = im.set('to_date', to_date.format("YYYY-MM-dd"))
+    im = im.set('l9_images', l9_subset.size())
+    
+    # im = ee.Algorithms.If(im.bandNames().size().eq(1), ee.Number(0), im)
+    
+    return ee.Image(im)
+
+
+def generate_timeseries(dates):
+    raw_ts = dates.map(process_date)
+    # raw_ts = raw_ts.removeAll([0])
+    imcoll = ee.ImageCollection.fromImages(raw_ts)
+
+    return imcoll
+
+def get_first_obs(start_date, end_date):
+    first_im = l9.filterBounds(aoi).filterDate(start_date, end_date).first()
+    str_fmt = 'YYYY-MM-dd'
+    return ee.Date.parse(str_fmt, ee.Date(first_im.get('system:time_start')).format(str_fmt))
+
+def run_process_long(res_name, res_polygon, start, end, datadir):
+    fo = start
+    enddate = end
+
+    # Extracting reservoir geometry 
+    global aoi
+    aoi = poly2feature(res_polygon,BUFFER_DIST).geometry()
+    ## Checking if time interval is small then the image collection should not be empty in GEE
+    if (days_between(start,end) < 30):     # less than a month difference
+        number_of_images = l9.filterBounds(aoi).filterDate(start, end).size().getInfo()
+    else:
+        number_of_images = 1     # more than a month difference simply run, so no need to calculate number_of_images
+    
+    if(number_of_images):
+        fo = get_first_obs(start, end).format('YYYY-MM-dd').getInfo()
+        first_obs = datetime.strptime(fo, '%Y-%m-%d')
+
+        scratchdir = os.path.join(datadir, "_scratch")
+
+        # flag = True
+        num_runs = 0
+
+        # If data already exists, only get new data starting from the last one
+        savepath = os.path.join(datadir, f"{res_name}.csv")
+        
+        if os.path.isfile(savepath):
+            temp_df = pd.read_csv(savepath, parse_dates=['mosaic_enddate']).set_index('mosaic_enddate')
+
+            last_date = temp_df.index[-1].to_pydatetime()
+            fo = (last_date - timedelta(days=TEMPORAL_RESOLUTION*2)).strftime("%Y-%m-%d")
+            to_combine = [savepath]
+            print(f"Existing file found - Last observation ({TEMPORAL_RESOLUTION*2} day lag): {last_date}")
+
+            # If 16 days have not passed since last observation, skip the processing
+            days_passed = (datetime.strptime(end, "%Y-%m-%d") - last_date).days
+            print(f"No. of days passed since: {days_passed}")
+            if days_passed < TEMPORAL_RESOLUTION:
+                print(f"No new observation expected. Quitting early")
+                return None
+        else:
+            to_combine = []
+
+        savedir = os.path.join(scratchdir, f"{res_name}_l9_cordeiro_zhao_gao_{fo}_{enddate}")
+        if not os.path.isdir(savedir):
+            os.makedirs(savedir)
+        
+        print(f"Extracting SA for the period {fo} -> {enddate}")
+
+        dates = pd.date_range(fo, enddate, freq=f'{TEMPORAL_RESOLUTION}D')
+        grouped_dates = grouper(dates, RESULTS_PER_ITER)
+
+        # # redo the calculations part and see where it is complaining about too many aggregations
+        # subset_dates = next(grouped_dates)
+        # dates = ee.List([ee.Date(d) for d in subset_dates if d is not None])
+
+        # print(subset_dates)
+        # res = generate_timeseries(dates).filterMetadata('s2_images', 'greater_than', 0)
+        # pprint.pprint(res.aggregate_array('s2_images').getInfo())
+
+        # uncorrected_columns_to_extract = ['from_date', 'to_date', 'water_area_cordeiro', 'non_water_area_cordeiro', 'water_area_NDWI', 'non_water_area_NDWI', 'cloud_area', 's2_images']
+        # uncorrected_final_data_ee = res.reduceColumns(ee.Reducer.toList(len(uncorrected_columns_to_extract)), uncorrected_columns_to_extract).get('list')
+        # uncorrected_final_data = uncorrected_final_data_ee.getInfo()
+        
+        for subset_dates in grouped_dates:
+            try:
+                print(subset_dates)
+                dates = ee.List([ee.Date(d) for d in subset_dates if d is not None])
+                
+                res = generate_timeseries(dates).filterMetadata('l9_images', 'greater_than', 0)
+                # pprint.pprint(res.getInfo())
+
+                uncorrected_columns_to_extract = ['from_date', 'to_date', 'water_area_cordeiro', 'non_water_area_cordeiro', 'cloud_area', 'l9_images']
+                uncorrected_final_data_ee = res.reduceColumns(ee.Reducer.toList(len(uncorrected_columns_to_extract)), uncorrected_columns_to_extract).get('list')
+                uncorrected_final_data = uncorrected_final_data_ee.getInfo()
+                print("Uncorrected", uncorrected_final_data)
+
+                res_corrected_cordeiro = res.map(lambda im: postprocess_wrapper(im, 'water_map_cordeiro', im.get('water_area_cordeiro')))
+                corrected_columns_to_extract = ['to_date', 'corrected_area']
+                corrected_final_data_cordeiro_ee = res_corrected_cordeiro \
+                                                    .filterMetadata('corrected_area', 'not_equals', None) \
+                                                    .reduceColumns(
+                                                        ee.Reducer.toList(
+                                                            len(corrected_columns_to_extract)), 
+                                                            corrected_columns_to_extract
+                                                            ).get('list')
+                corrected_final_data_cordeiro = corrected_final_data_cordeiro_ee.getInfo()
+                print("Corrected - Cordeiro", corrected_final_data_cordeiro)
+
+                # res_corrected_NDWI = res.map(lambda im: postprocess_wrapper(im, 'water_map_NDWI', im.get('water_area_NDWI')))
+                # corrected_final_data_NDWI_ee = res_corrected_NDWI \
+                #                                     .filterMetadata('corrected_area', 'not_equals', None) \
+                #                                     .reduceColumns(
+                #                                         ee.Reducer.toList(
+                #                                             len(corrected_columns_to_extract)), 
+                #                                             corrected_columns_to_extract
+                #                                             ).get('list')
+                
+                # corrected_final_data_NDWI = corrected_final_data_NDWI_ee.getInfo()
+                
+                # print(uncorrected_final_data, corrected_final_data_cordeiro)
+                if len(uncorrected_final_data) == 0:
+                    continue
+                uncorrected_df = pd.DataFrame(uncorrected_final_data, columns=uncorrected_columns_to_extract)
+                corrected_cordeiro_df = pd.DataFrame(corrected_final_data_cordeiro, columns=corrected_columns_to_extract).rename({'corrected_area': 'corrected_area_cordeiro'}, axis=1)
+                # corrected_NDWI_df = pd.DataFrame(corrected_final_data_NDWI, columns=corrected_columns_to_extract).rename({'corrected_area': 'corrected_area_NDWI'}, axis=1)
+                # corrected_df = pd.merge(corrected_cordeiro_df, corrected_NDWI_df, 'left', 'to_date')
+                df = pd.merge(uncorrected_df, corrected_cordeiro_df, 'left', 'to_date')
+
+                df['from_date'] = pd.to_datetime(df['from_date'], format="%Y-%m-%d")
+                df['to_date'] = pd.to_datetime(df['to_date'], format="%Y-%m-%d")
+                df['mosaic_enddate'] = df['to_date'] - pd.Timedelta(1, unit='day')
+                df = df.set_index('mosaic_enddate')
+                print(df.head(2))
+
+                fname = os.path.join(savedir, f"{df.index[0].strftime('%Y%m%d')}_{df.index[-1].strftime('%Y%m%d')}_{res_name}.csv")
+                df.to_csv(fname)
+
+                s_time = randint(20, 30)
+                print(f"Sleeping for {s_time} seconds")
+                time.sleep(randint(20, 30))
+
+                if (datetime.strptime(enddate, "%Y-%m-%d")-df.index[-1]).days < TEMPORAL_RESOLUTION:
+                    print(f"Quitting: Reached enddate {enddate}")
+                    break
+                elif df.index[-1].strftime('%Y-%m-%d') == fo:
+                    print(f"Reached last available observation - {fo}")
+                    break
+                elif num_runs > 1000:
+                    print("Quitting: Reached 1000 iterations")
+                    break
+            except Exception as e:
+                log.error(e)
+                continue
+        
+
+
+        # Combine the files into one database
+        to_combine.extend([os.path.join(savedir, f) for f in os.listdir(savedir) if f.endswith(".csv")])
+
+        files = [pd.read_csv(f, parse_dates=["mosaic_enddate"]).set_index("mosaic_enddate") for f in to_combine]
+        data = pd.concat(files).drop_duplicates().sort_values("mosaic_enddate")
+
+        data.to_csv(savepath)
+
+        return savepath
+    
+    else:
+        print(f"No observation observed between {start} and {end}. Quitting!")
+        return None
+
+# User-facing wrapper function
+def sarea_l9(res_name,res_polygon, start, end, datadir):
+    return run_process_long(res_name,res_polygon, start, end, datadir)
+
+def main():
+    # Setup argument parser
+    parser = argparse.ArgumentParser(description="Generate Reservoir Surface area time series")
+
+    parser.add_argument("reservoir")
+    parser.add_argument("start_date")
+    parser.add_argument("end_date")
+
+    args = parser.parse_args()
+
+    reservoir = args.reservoir
+    start = args.start_date
+    end = args.end_date
+
+    # print(f"{reservoir = }, {start = }, {end = }")
+    run_process_long(reservoir, start, end, "/houston2/pritam/rat_mekong_v3/backend/data/sarea_tmsos/l9")
+
+
+
+if __name__ == '__main__':
+    main()
diff --git a/src/rat/data_processing/metsim_input_processing.py b/src/rat/data_processing/metsim_input_processing.py
index 184df82c..e4607f7c 100644
--- a/src/rat/data_processing/metsim_input_processing.py
+++ b/src/rat/data_processing/metsim_input_processing.py
@@ -14,7 +14,7 @@
 log = getLogger(LOG_NAME)
 
 class CombinedNC:
-    def __init__(self, start, end, datadir, basingridpath, outputdir):
+    def __init__(self, start, end, datadir, basingridpath, outputdir, use_previous):
         """
         Parameters:
             start: Start date in YYYY-MM-DD format
@@ -28,6 +28,7 @@ def __init__(self, start, end, datadir, basingridpath, outputdir):
         self._total_days = (self._end - self._start).days
         self._datadir = datadir
         self._outputpath = outputdir
+        self._use_previous = use_previous
 
         self._rast = rio.open(basingridpath)
         self._ar = self._rast.read(1, masked=True)
@@ -150,19 +151,22 @@ def _write(self):
             )
         )
 
-        if os.path.isfile(self._outputpath):
-            log.debug(f"Found existing file at {self._outputpath} -- Updating in-place")
-            # Assuming the existing file structure is same as the one generated now. Basically
-            #   assuming that the previous file was also created by MetSimRunner
-            existing = xr.open_dataset(self._outputpath).load()
-            existing.close()
-            last_existing_time = existing.time[-1]
-            log.debug("Existing data: %s", last_existing_time)
-            existing_to_append = existing.sel(time=slice(ds.time[0] - np.timedelta64(120,'D') , last_existing_time))
-            ds = ds.sel(time=slice(last_existing_time + np.timedelta64(1,'D') , ds.time[-1]))
-            # ds = ds.isel(time=slice(1, None))
-            write_ds = xr.merge([existing_to_append, ds])
-            ds = self._impute_basin_missing_data(write_ds)
+        if self._use_previous:
+            if os.path.isfile(self._outputpath):
+                log.debug(f"Found existing file at {self._outputpath} -- Updating in-place")
+                # Assuming the existing file structure is same as the one generated now. Basically
+                #   assuming that the previous file was also created by MetSimRunner
+                existing = xr.open_dataset(self._outputpath).load()
+                existing.close()
+                last_existing_time = existing.time[-1]
+                log.debug("Existing data: %s", last_existing_time)
+                existing_to_append = existing.sel(time=slice(ds.time[0] - np.timedelta64(120,'D') , last_existing_time))
+                ds = ds.sel(time=slice(last_existing_time + np.timedelta64(1,'D') , ds.time[-1]))
+                # ds = ds.isel(time=slice(1, None))
+                write_ds = xr.merge([existing_to_append, ds])
+                ds = self._impute_basin_missing_data(write_ds)
+            else:
+                raise Exception('Previous combined dataset not found. Please run RAT without state files first.')
         else:
             log.debug(f"Creating new file at {self._outputpath}")
             ds = self._impute_basin_missing_data(ds)
diff --git a/src/rat/rat_basin.py b/src/rat/rat_basin.py
index 5fa411c5..c04d32e8 100644
--- a/src/rat/rat_basin.py
+++ b/src/rat/rat_basin.py
@@ -24,14 +24,14 @@
 from rat.utils.route_param_reader import RouteParameterFile
 from rat.core.run_routing import RoutingRunner
 
-from rat.core.run_sarea import run_sarea
 from rat.core.run_altimetry import run_altimetry
 
-from rat.ee_utils.ee_aec_file_creator import aec_file_creator
 from rat.core.run_postprocessing import run_postprocessing
 
-from rat.utils.convert_for_website import convert_sarea, convert_inflow, convert_dels, convert_evaporation, convert_outflow, convert_altimeter
+from rat.utils.convert_to_final_outputs import convert_sarea, convert_inflow, convert_dels, convert_evaporation, convert_outflow, convert_altimeter
 
+# Step-(-1): Reading Configuration settings to run RAT
+# Step-0: Creating required directory structure for RAT
 # Step-1: Downloading and Pre-processing of meteorolgical data
 # Step-2: Pre-processing of data and preparation of MetSim Input
 # Step-3: Preparation of MetSim Parameter Files
@@ -45,98 +45,137 @@
 # Step-11: Elevation extraction from Altimeter
 # Step-12: Generating Area Elevation Curves for reservoirs
 # Step-13: Calculation of Outflow, Evaporation and Storage change
+# Step-14: Conversion of output data to final format as time series
 
+#TODO: Converting steps to separate modules to make RAT more robust and generalized
 #module-1 step-1,2 data_preparation_vic
 #module-2 step-3to8 inflow_vic
 #module-3 step- NA evaporation
 #module-4 step-10to12 storage_change
-#module-5 step-13 outflow
+#module-5 step-13 outflow (mass-balance)
+#RAT using all modules and step-14 to produce final outputs
 
-def rat(config, rat_logger, steps=[1,2,3,4,5,6,7,8,9,10,11,12,13]):
+def rat_basin(config, rat_logger):
+    """Runs RAT as per configuration defined in `config_fn` for one single basin.
 
-    rat_logger = getLogger('run_rat')
-    ##--------------------- Reading and initialising global parameters ----------------------##
-    # Defining resolution to run RAT
-    xres=0.0625
-    yres=0.0625 
-
-    # Obtaining basin related information from RAT_Runner.yml
-    basins_shapefile_path = config['GLOBAL']['basin_shpfile'] # Shapefile containg information of basin(s)- geometry and attributes
-    basins_shapefile = gpd.read_file(basins_shapefile_path)  # Reading basins_shapefile_path to get basin polygons and their attributes
-    basins_shapefile_column_dict = config['GLOBAL']['basin_shpfile_column_dict'] # Dictionary of column names in basins_shapefile, Must contain 'id' field
-    region_name = config['BASIN']['region_name']  # Major basin name used to cluster multiple basins data in data-directory
-    basin_name = config['BASIN']['basin_name']              # Basin name used to save basin related data
-    basin_id = config['BASIN']['basin_id']                  # Unique identifier for each basin used to map basin polygon in basins_shapefile
-    basin_data = basins_shapefile[basins_shapefile[basins_shapefile_column_dict['id']]==basin_id] # Getting the particular basin related information corresponding to basin_id
-    basin_bounds = basin_data.bounds                          # Obtaining bounds of the particular basin
-    basin_bounds = np.array(basin_bounds)[0]
-    basin_geometry = basin_data.geometry                      # Obtaining geometry of the particular basin
-
-    # Defining paths for RAT processing
-    project_dir = config['GLOBAL']['project_dir']    # Directory of RAT 
-    data_dir = config['GLOBAL']['data_dir']          # Data-Directory of RAT
-    region_data_dir = create_directory(os.path.join(data_dir,region_name), True) # Major Basin data-directory within the data-directory of RAT
-    basin_data_dir = create_directory(os.path.join(region_data_dir,'basins',basin_name), True)  # Basin data-directory within the major basin's data-directory 
-    log_dir = create_directory(os.path.join(region_data_dir,'logs',basin_name,''), True)  # Log directory within the major basin's data-directory
-
-    # Change datetimes format
-    #config['BASIN']['begin'] = datetime.datetime.combine(config['BASIN']['begin'], datetime.time.min)
-    config['BASIN']['start'] = datetime.datetime.combine(config['BASIN']['start'], datetime.time.min)
-    config['BASIN']['end'] = datetime.datetime.combine(config['BASIN']['end'], datetime.time.min)
-    rout_init_state_save_date = config['BASIN']['end']+datetime.timedelta(days=1)
-
-    if(not config['BASIN']['first_run']):
-        if(config['BASIN'].get('vic_init_state_date')):
-            config['BASIN']['vic_init_state_date'] = datetime.datetime.combine(config['BASIN']['vic_init_state_date'], datetime.time.min)
+    parameters:
+        config_fn (str): Path to the configuration file for one basin and not for multiple basins. To run rat for multiple basins use run_rat()
     
-    # Changing start date if running RAT for first time for the particular basin to give VIC and MetSim to have their spin off periods
-    if(config['BASIN']['first_run']):
-        user_given_start = config['BASIN']['start']
-        config['BASIN']['start'] = user_given_start-datetime.timedelta(days=800)  # Running RAT for extra 800 days before the user-given start date for VIC to give reliable results starting from user-given start date
-        data_download_start = config['BASIN']['start']-datetime.timedelta(days=90)    # Downloading 90 days of extra meteorological data for MetSim to prepare it's initial state
-        vic_init_state_date = None    # No initial state of VIC is present as running RAT for first time in this basin
-    elif(config['BASIN'].get('vic_init_state_date')):
-        data_download_start = config['BASIN']['start']    # Downloading data from the same date as we want to run RAT from
-        vic_init_state_date = config['BASIN']['vic_init_state_date'] # Date of which initial state of VIC for the particular basin exists
+    returns:
+        no_errors (int): Number of errors/exceptions occurred while running rat
+        latest_altimetry_cycle (int): Latest altimetry jason3 cycle number if rat runs step 11
+    """
+    rat_logger = getLogger('run_rat')
+
+    # Tracking number of errors to output
+    no_errors = 0
+
+    # Tracking latest altimetry cycle number if step 11 is run, otherwise return None
+    latest_altimetry_cycle = None
+
+    ######### Step -1 Mandatory Step
+    try:
+        rat_logger.info("Reading Configuration settings to run RAT")
+        ##--------------------- Reading and initialising global parameters ----------------------##
+
+        # Reading steps to be executed in RAT otherwise seting default value
+        if config['GLOBAL'].get('steps'):
+            steps = config['GLOBAL']['steps']
+        else:
+            steps = [1,2,3,4,5,6,7,8,9,10,11,12,13,14]
+
+        # Defining resolution to run RAT
+        xres = 0.0625
+        yres = 0.0625 
+
+        # Obtaining basin related information from RAT_Runner.yml
+        basins_shapefile_path = config['GLOBAL']['basin_shpfile'] # Shapefile containg information of basin(s)- geometry and attributes
+        basins_shapefile = gpd.read_file(basins_shapefile_path)  # Reading basins_shapefile_path to get basin polygons and their attributes
+        basins_shapefile_column_dict = config['GLOBAL']['basin_shpfile_column_dict'] # Dictionary of column names in basins_shapefile, Must contain 'id' field
+        region_name = config['BASIN']['region_name']  # Major basin name used to cluster multiple basins data in data-directory
+        basin_name = config['BASIN']['basin_name']              # Basin name used to save basin related data
+        basin_id = config['BASIN']['basin_id']                  # Unique identifier for each basin used to map basin polygon in basins_shapefile
+        basin_data = basins_shapefile[basins_shapefile[basins_shapefile_column_dict['id']]==basin_id] # Getting the particular basin related information corresponding to basin_id
+        basin_bounds = basin_data.bounds                          # Obtaining bounds of the particular basin
+        basin_bounds = np.array(basin_bounds)[0]
+        basin_geometry = basin_data.geometry                      # Obtaining geometry of the particular basin
+
+        # Defining paths for RAT processing
+        project_dir = config['GLOBAL']['project_dir']    # Directory of RAT 
+        data_dir = config['GLOBAL']['data_dir']          # Data-Directory of RAT
+        region_data_dir = create_directory(os.path.join(data_dir,region_name), True) # Major Basin data-directory within the data-directory of RAT
+        basin_data_dir = create_directory(os.path.join(region_data_dir,'basins',basin_name), True)  # Basin data-directory within the major basin's data-directory 
+        log_dir = create_directory(os.path.join(region_data_dir,'logs',basin_name,''), True)  # Log directory within the major basin's data-directory
+
+        # Change datetimes format
+        #config['BASIN']['begin'] = datetime.datetime.combine(config['BASIN']['begin'], datetime.time.min)
+        config['BASIN']['start'] = datetime.datetime.combine(config['BASIN']['start'], datetime.time.min)
+        config['BASIN']['end'] = datetime.datetime.combine(config['BASIN']['end'], datetime.time.min)
+        rout_init_state_save_date = config['BASIN']['end']
+
+        if(not config['BASIN']['spin_up']):
+            if(config['BASIN'].get('vic_init_state_date')):
+                config['BASIN']['vic_init_state_date'] = datetime.datetime.combine(config['BASIN']['vic_init_state_date'], datetime.time.min)
+        
+        # Changing start date if running RAT for first time for the particular basin to give VIC and MetSim to have their spin off periods
+        if(config['BASIN']['spin_up']):
+            user_given_start = config['BASIN']['start']
+            config['BASIN']['start'] = user_given_start-datetime.timedelta(days=800)  # Running RAT for extra 800 days before the user-given start date for VIC to give reliable results starting from user-given start date
+            data_download_start = config['BASIN']['start']-datetime.timedelta(days=90)    # Downloading 90 days of extra meteorological data for MetSim to prepare it's initial state
+            vic_init_state_date = None    # No initial state of VIC is present as running RAT for first time in this basin
+            use_state = False            # Routing state file won't be used
+        elif(config['BASIN'].get('vic_init_state_date')):
+            data_download_start = config['BASIN']['start']    # Downloading data from the same date as we want to run RAT from
+            vic_init_state_date = config['BASIN']['vic_init_state_date'] # Date of which initial state of VIC for the particular basin exists
+            use_state = True           # Routing state file will be used
+        else:
+            data_download_start = config['BASIN']['start']-datetime.timedelta(days=90)    # Downloading 90 days of extra meteorological data for MetSim to prepare it's initial state
+            vic_init_state_date = None    # No initial state of VIC is present as running RAT for first time in this basin
+            use_state = False            # Routing state file won't be used
+
+        # Defining logger
+        log = init_logger(
+            log_dir= log_dir,
+            verbose= False,
+            # notify= True,
+            notify= False,
+            log_level= 'DEBUG'
+        )
+
+        # Cleaning class object for removing/deleting unwanted data
+        cleaner = Clean(basin_data_dir)
+
+        # Clearing out previous rat outputs so that the new data does not gets appended.
+        if(config['CLEAN_UP']['clean_previous_outputs']):
+            rat_logger.info("Clearing up memory space: Removal of previous rat outputs, routing inflow, extracted altimetry data and gee extracted surface area time series")
+            cleaner.clean_previous_outputs()
+
+        # Initializing Status for different models & tasks (1 for successful run & 0 for failed run)
+        NEW_DATA_STATUS = 1
+        METSIM_STATUS = 1
+        VIC_STATUS = 1
+        ROUTING_STATUS = 1
+        GEE_STATUS = 1
+        ALTIMETER_STATUS = 1
+        DELS_STATUS = 0
+        EVAP_STATUS = 0
+        OUTFLOW_STATUS = 0
+    except:
+        no_errors = -1
+        rat_logger.exception("Error in Configuration parameters defined to run RAT.")
+        return (no_errors, latest_altimetry_cycle)
     else:
-        data_download_start = config['BASIN']['start']-datetime.timedelta(days=90)    # Downloading 90 days of extra meteorological data for MetSim to prepare it's initial state
-        vic_init_state_date = None    # No initial state of VIC is present as running RAT for first time in this basin
-    
-    # Defining logger
-    log = init_logger(
-        log_dir= log_dir,
-        verbose= False,
-        # notify= True,
-        notify= False,
-        log_level= 'DEBUG'
-    )
-
-    # Cleaning class object for removing/deleting unwanted data
-    cleaner = Clean(basin_data_dir)
-
-    # Clearing out previous rat outputs so that the new data does not gets appended.
-    if(config['CLEAN_UP']['clean_previous_outputs']):
-        rat_logger.info("Clearing up memory space: Removal of previous rat outputs, routing inflow, extracted altimetry data and gee extracted surface area time series")
-        cleaner.clean_previous_outputs()
-
-    # Initializing Status for different models & tasks (1 for successful run & 0 for failed run)
-    NEW_DATA_STATUS = 1
-    METSIM_STATUS = 1
-    VIC_STATUS = 1
-    ROUTING_STATUS = 1
-    GEE_STATUS = 1
-    ALTIMETER_STATUS = 1
-    DELS_STATUS = 0
-    EVAP_STATUS = 0
-    OUTFLOW_STATUS = 0
-
-    ##--------------------- Read and initialised global parameters ----------------------##
-
-    rat_logger.info(f"Running RAT from {config['BASIN']['start']} to {config['BASIN']['end']}")
+        rat_logger.info("Read Configuration settings to run RAT.")
+        ##--------------------- Read and initialised global parameters ----------------------##
+
+    rat_logger.info(f"Running RAT from {config['BASIN']['start'] } to {config['BASIN']['end']}")
 
     ######### Step-0 Mandatory Step
     try:
-        rat_logger.info("Starting Step-0: Creating required directory structure for RAT")
+        rat_logger.info("Creating required directory structure for RAT")
+        
+        ##--------------------- Defining global paths and variables ----------------------##
+
         #----------- Paths Necessary for running of METSIM  -----------#
         # Path of directory which will contain the combined data in nc format.
         combined_datapath = create_directory(os.path.join(basin_data_dir,'pre_processing','nc', ''), True)
@@ -184,6 +223,11 @@ def rat(config, rat_logger, steps=[1,2,3,4,5,6,7,8,9,10,11,12,13]):
         #----------- Paths Necessary for running of Routing  -----------#
 
         #----------- Paths Necessary for running of Surface Area Calculation and Altimetry-----------#
+        # Defining routing station file
+        if(config['ROUTING PARAMETERS'].get('station_file')):
+            basin_station_xy_path = config['ROUTING PARAMETERS'].get('station_file')
+        else:
+            basin_station_xy_path = os.path.join(basin_data_dir,'ro', 'pars','sta_xy.txt')
         # Defining path for basin reservoir shapefile
         basin_reservoir_shpfile_path = create_directory(os.path.join(basin_data_dir,'gee','gee_basin_params',''), True)
         basin_reservoir_shpfile_path = os.path.join(basin_reservoir_shpfile_path,'basin_reservoirs.shp')
@@ -212,13 +256,13 @@ def rat(config, rat_logger, steps=[1,2,3,4,5,6,7,8,9,10,11,12,13]):
         final_output_path = create_directory(os.path.join(basin_data_dir,'final_outputs',''),True)
         ## End of defining paths for storing post-processed data and webformat data
         #----------- Paths Necessary for running of Post-Processing-----------#
-
-
     except:
-        rat_logger.exception("Error Executing Step-0: Downloading and Pre-processing of meteorological data")
+        no_errors = -1
+        rat_logger.exception("Error in creating required directory structure for RAT")
+        return (no_errors, latest_altimetry_cycle)
     else:
-        rat_logger.info("Finished Step-0: Downloading and Pre-processing of meteorological data")
-        
+        rat_logger.info("Finished creating required directory structure for RAT")
+        ##--------------------- Definied global paths and variables ----------------------##
 
     ######### Step-1
     if(1 in steps):
@@ -240,6 +284,7 @@ def rat(config, rat_logger, steps=[1,2,3,4,5,6,7,8,9,10,11,12,13]):
             ##---------- Download Data End -----------# 
             NEW_DATA_STATUS = 1   #Data downloaded successfully
         except:
+            no_errors = no_errors+1
             rat_logger.exception("Error Executing Step-1: Downloading and Pre-processing of meteorological data")
         else:
             rat_logger.info("Finished Step-1: Downloading and Pre-processing of meteorological data")
@@ -263,7 +308,8 @@ def rat(config, rat_logger, steps=[1,2,3,4,5,6,7,8,9,10,11,12,13]):
                 end= config['BASIN']['end'],
                 datadir= processed_datadir,
                 basingridpath= basingridfile_path,
-                outputdir= combined_datapath
+                outputdir= combined_datapath,
+                use_previous= use_state,
             )
             #----------- Process Data End and combined data created -----------#
 
@@ -277,6 +323,7 @@ def rat(config, rat_logger, steps=[1,2,3,4,5,6,7,8,9,10,11,12,13]):
             )
             #------- MetSim Input Data Preparation End -------#
         except:
+            no_errors = no_errors+1
             rat_logger.exception("Error Executing Step-2: Pre-processing of data and preparation of MetSim Input")
         else:
             rat_logger.info("Finished Step-2: Pre-processing of data and preparation of MetSim Input")
@@ -288,11 +335,15 @@ def rat(config, rat_logger, steps=[1,2,3,4,5,6,7,8,9,10,11,12,13]):
             rat_logger.info("Starting Step-3: Preparation of MetSim Parameter Files")
             ##----------------------- Preparing parameter files for METSIM ------------------------##
             ## ---------- Creating domain.nc file for basin if not exist---------#
-            if not os.path.exists(domain_nc_path):
-                elevation_tif_filepath = config['GLOBAL']['elevation_tif_file']
-                create_basin_domain_nc_file(elevation_tif_filepath,basingridfile_path,domain_nc_path)
+            if(config['GLOBAL'].get('elevation_tif_file')):
+                if not os.path.exists(domain_nc_path):
+                    elevation_tif_filepath = config['GLOBAL']['elevation_tif_file']
+                    create_basin_domain_nc_file(elevation_tif_filepath,basingridfile_path,domain_nc_path)
+            else:
+                domain_nc_path = config['METSIM']['metsim_domain_file']
             #----------domain.nc file created for basin if not exist ---------#
         except:
+            no_errors = no_errors+1
             rat_logger.exception("Error Executing Step-3: Preparation of MetSim Parameter Files")
         else:
             rat_logger.info("Finished Step-3: Preparation of MetSim Parameter Files")
@@ -327,6 +378,7 @@ def rat(config, rat_logger, steps=[1,2,3,4,5,6,7,8,9,10,11,12,13]):
             else:
                 rat_logger.info("New Data Download Failed. Skipping Step-4: Running MetSim & preparation of VIC input")
         except:
+            no_errors = no_errors+1
             rat_logger.exception("Error Executing Step-4: Running MetSim & preparation of VIC input")
         else:
             rat_logger.info("Finished Step-4: Running MetSim & preparation of VIC input")
@@ -355,6 +407,7 @@ def rat(config, rat_logger, steps=[1,2,3,4,5,6,7,8,9,10,11,12,13]):
                 config['VIC']['vic_soil_param_file']=os.path.join(vic_param_dir,'vic_soil_param.nc')
                 config['VIC']['vic_domain_file']=os.path.join(vic_param_dir,'vic_domain.nc')
         except:
+            no_errors = no_errors+1
             rat_logger.exception("Error Executing Step-5: Preparation of VIC Parameter Files")
         else:
             rat_logger.info("Finished Step-5: Preparation of VIC Parameter Files")
@@ -382,21 +435,21 @@ def rat(config, rat_logger, steps=[1,2,3,4,5,6,7,8,9,10,11,12,13]):
                         rout_input_dir= rout_input_path
                     )
                     vic.run_vic(np=config['GLOBAL']['multiprocessing'])
-                    rout_input_state_start_date = vic.generate_routing_input_state(ndays=365, rout_input_state_file=rout_input_state_file, save_path=rout_init_state_save_file) # Start date of routing state file will be returned
-                    if(config['BASIN']['first_run']):
-                        vic.disagg_results(rout_input_state_file=p.vic_result_file)       # If first run, use vic result file
-                        rout_input_state_start_date = config['BASIN']['start']            # Start date will become same as vic result file 
+                    vic.generate_routing_input_state(ndays=365, rout_input_state_file=rout_input_state_file, 
+                                                                                        save_path=rout_init_state_save_file, use_rout_state=use_state) # Start date of routing state file will be returned
+                    if(config['BASIN']['spin_up']):
+                        vic.disagg_results(rout_input_state_file=p.vic_result_file)       # If spin_up, use vic result file
                     elif(config['BASIN'].get('vic_init_state_date')):                      # If vic_state file exists
-                        vic.disagg_results(rout_input_state_file=rout_input_state_file)    # If not first run, use rout input state file
+                        vic.disagg_results(rout_input_state_file=rout_init_state_save_file)    # If not spin_up, use rout_init_state_save file just creates as it contains the recent vic outputs.
                     else:
-                        vic.disagg_results(rout_input_state_file=p.vic_result_file)       # If first run is false and vic state file does not exist, use vic result file
-                        rout_input_state_start_date = config['BASIN']['start']            # Start date will become same as vic result file 
+                        vic.disagg_results(rout_input_state_file=p.vic_result_file)       # If spin_up is false and vic state file does not exist, use vic result file
                     vic_startdate = p.vic_startdate
                     vic_enddate = p.vic_enddate
                 VIC_STATUS=1         #Vic run successfully
             else:
                 rat_logger.info("MetSim Run Failed. Skipping Step-6: Running of VIC and preparation of Routing input")
         except:
+            no_errors = no_errors+1
             rat_logger.exception("Error Executing Step-6: Running of VIC and preparation of Routing input")
         else:
             rat_logger.info("Finished Step-6: Running of VIC and preparation of Routing input")
@@ -415,12 +468,12 @@ def rat(config, rat_logger, steps=[1,2,3,4,5,6,7,8,9,10,11,12,13]):
                     create_basin_grid_flow_asc(config['ROUTING']['global_flow_dir_tif_file'], basingridfile_path, basin_flow_dir_file[:-4],
                                                                     config['ROUTING'].get('replace_flow_directions'))
             ### Basin Station File
-            basin_station_latlon_file = os.path.join(rout_param_dir,'basin_station_latlon.csv')
             if (config['ROUTING']['station_global_data']):
                 if not os.path.exists(basin_station_latlon_file):
                     create_basin_station_latlon_csv(region_name,basin_name, config['ROUTING']['stations_vector_file'], basin_data, 
                                                         config['ROUTING']['stations_vector_file_columns_dict'], basin_station_latlon_file)
         except:
+            no_errors = no_errors+1
             rat_logger.exception("Error Executing Step-7: Preparation of Routing Parameter Files")
         else:    
             rat_logger.info("Finished Step-7: Preparation of Routing Parameter Files")
@@ -432,10 +485,10 @@ def rat(config, rat_logger, steps=[1,2,3,4,5,6,7,8,9,10,11,12,13]):
             if(VIC_STATUS):
                 rat_logger.info("Starting Step-8: Runnning Routing and generating Inflow")
                 ### Extracting routing start date ###
-                if(config['BASIN']['first_run']):
+                if(config['BASIN']['spin_up']):
                     rout_input_state_start_date = config['BASIN']['start']
                 elif(config['BASIN'].get('vic_init_state_date')): 
-                    rout_state_data = xr.open_dataset(rout_input_state_file).load()
+                    rout_state_data = xr.open_dataset(rout_init_state_save_file).load()
                     rout_input_state_start_date = rout_state_data.time[0].values.astype('datetime64[us]').astype(datetime.datetime)
                     rout_state_data.close()
                 else:
@@ -469,6 +522,7 @@ def rat(config, rat_logger, steps=[1,2,3,4,5,6,7,8,9,10,11,12,13]):
             else:
                 rat_logger.info("VIC Run Failed. Skipping Step-8: Runnning Routing and generating Inflow")
         except:
+            no_errors = no_errors+1
             rat_logger.exception("Error Executing Step-8: Runnning Routing and generating Inflow")
         else:
             rat_logger.info("Finished Step-8: Runnning Routing and generating Inflow")
@@ -482,10 +536,12 @@ def rat(config, rat_logger, steps=[1,2,3,4,5,6,7,8,9,10,11,12,13]):
             ###### Preparing basin's reservoir shapefile and it's associated column dictionary for calculating surface area #####
             ### Creating Basin Reservoir Shapefile, if not exists ###
             if not os.path.exists(basin_reservoir_shpfile_path):
-                create_basin_reservoir_shpfile(config['GEE']['reservoir_vector_file'], reservoirs_gdf_column_dict, basin_station_xy_path,
+                if os.path.exists(basin_station_xy_path):
+                    create_basin_reservoir_shpfile(config['GEE']['reservoir_vector_file'], reservoirs_gdf_column_dict, basin_station_xy_path,
                                                                                 config['ROUTING']['station_global_data'], basin_reservoir_shpfile_path)
             ###### Prepared basin's reservoir shapefile and it's associated column dictionary #####
         except:
+            no_errors = no_errors+1
             rat_logger.exception("Error Executing Step-9: Preparation of parameter files for Surface Area Calculation")
         else:
             rat_logger.info("Finished Step-9: Preparation of parameter files for Surface Area Calculation")
@@ -494,13 +550,20 @@ def rat(config, rat_logger, steps=[1,2,3,4,5,6,7,8,9,10,11,12,13]):
     ######### Step-10
     if(10 in steps):
         try:
+            from rat.core.run_sarea import run_sarea
             rat_logger.info("Starting Step-10: TMS-OS Surface Area Calculation from GEE")
             ##----------- Remote Sensing to estimate surface area begins -----------##
+            if (not os.path.exists(basin_reservoir_shpfile_path)):
+                if (not config['ROUTING']['station_global_data']):
+                    basin_reservoir_shpfile_path = config['GEE']['reservoir_vector_file']
+                else: 
+                    raise Exception('There was an error in creating reservoir shapefile using spatial join for this basin from the global reservoir vector file.')
             # Get Sarea
             run_sarea(config['BASIN']['start'].strftime("%Y-%m-%d"), config['BASIN']['end'].strftime("%Y-%m-%d"), sarea_savepath, 
                                                                                     basin_reservoir_shpfile_path, reservoirs_gdf_column_dict)
             GEE_STATUS = 1
         except:
+            no_errors = no_errors+1
             rat_logger.exception("Error Executing Step-10: TMS-OS Surface Area Calculation from GEE")
         else:
             rat_logger.info("Finished Step-10: TMS-OS Surface Area Calculation from GEE")                                                                        
@@ -510,20 +573,22 @@ def rat(config, rat_logger, steps=[1,2,3,4,5,6,7,8,9,10,11,12,13]):
     if(11 in steps):
         try:
             rat_logger.info("Starting Step-11: Elevation extraction from Altimeter")
-            ##----------- Altimeter height ectraction begins -----------##
+            ##----------- Altimeter height extraction begins -----------##
             # Altimeter
             latest_altimetry_cycle = run_altimetry(config, 'ALTIMETER', basin_reservoir_shpfile_path, reservoirs_gdf_column_dict, 
                                                                                     basin_name, basin_data_dir, altimetry_savepath)
             ALTIMETER_STATUS = 1
         except:
+            no_errors = no_errors+1
             rat_logger.exception("Error Executing Step-11: Elevation extraction from Altimeter")
         else:
             rat_logger.info("Finished Step-11: Elevation extraction from Altimeter")                                                                        
-            ##----------- Altimeter height ectraction ends -----------##
+            ##----------- Altimeter height extraction ends -----------##
 
     ######### Step-12
     if(12 in steps):
         try:
+            from rat.ee_utils.ee_aec_file_creator import aec_file_creator
             rat_logger.info("Starting Step-12: Generating Area Elevation Curves for reservoirs")
             ##--------------------------------Area Elevation Curves Extraction begins ------------------- ##
             ## Creating AEC files if not present for post-processing dels calculation
@@ -539,11 +604,23 @@ def rat(config, rat_logger, steps=[1,2,3,4,5,6,7,8,9,10,11,12,13]):
         try:
             rat_logger.info("Starting Step-13: Calculation of Outflow, Evaporation and Storage change")
             
-            ##---------- Mass-balance Approach begins and then post-processing ----------##
+            ##---------- Mass-balance Approach begins and then post-processing ----------## 
             DELS_STATUS, EVAP_STATUS, OUTFLOW_STATUS = run_postprocessing(basin_name, basin_data_dir, basin_reservoir_shpfile_path, reservoirs_gdf_column_dict,
-                                aec_dir_path, config['BASIN']['start'], config['BASIN']['end'], evap_savedir, dels_savedir, outflow_savedir, VIC_STATUS, ROUTING_STATUS, GEE_STATUS)
+                                aec_dir_path, config['BASIN']['start'], config['BASIN']['end'], rout_init_state_save_file, use_state, evap_savedir, dels_savedir, outflow_savedir, VIC_STATUS, ROUTING_STATUS, GEE_STATUS)
 
-            # Convert to format that is expected by the website and save in web_dir
+        except:
+            no_errors = no_errors+1
+            rat_logger.exception("Error Executing Step-13: Calculation of Outflow, Evaporation and Storage change")
+        else:
+            rat_logger.info("Finished Step-13: Calculation of Outflow, Evaporation and Storage change")
+            ##---------- Mass-balance Approach ends and then post-processed outputs to obtain timeseries  -----------------##
+    
+    ######### Step-14
+    if(14 in steps):
+        try:
+            rat_logger.info("Starting Step-14: Creating final outputs in a timeseries format and cleaning up.")
+
+            ##---------- Convert all time-series to final output csv format and clean up----------## 
             ## Surface Area
             if(GEE_STATUS):
                 convert_sarea(sarea_savepath,final_output_path)
@@ -585,27 +662,29 @@ def rat(config, rat_logger, steps=[1,2,3,4,5,6,7,8,9,10,11,12,13]):
                 rat_logger.info("Converted Extracted Height from Altimeter to the Output Format.")
             else:
                 rat_logger.info("Could not convert Extracted Height from Altimeter to the Output Format as Altimeter Run failed.")
+            
+            # Clearing out memory space as per user input 
+            if(config['CLEAN_UP'].get('clean_metsim')):
+                rat_logger.info("Clearing up memory space: Removal of metsim output files")
+                cleaner.clean_metsim()
+            if(config['CLEAN_UP'].get('clean_vic')):
+                rat_logger.info("Clearing up memory space: Removal of vic input, output files and previous init_state_files")
+                cleaner.clean_vic()
+            if(config['CLEAN_UP'].get('clean_routing')):
+                rat_logger.info("Clearing up memory space: Removal of routing input and output files")
+                cleaner.clean_routing()
+            if(config['CLEAN_UP'].get('clean_gee')):
+                rat_logger.info("Clearing up memory space: Removal of unwanted gee extracted small chunk files")
+                cleaner.clean_gee()
+            if(config['CLEAN_UP'].get('clean_altimetry')):
+                rat_logger.info("Clearing up memory space: Removal of raw altimetry downloaded data files.")
+                cleaner.clean_altimetry()
         except:
-            rat_logger.exception("Error Executing Step-13: Calculation of Outflow, Evaporation and Storage change")
+            no_errors = no_errors+1
+            rat_logger.exception("Error Executing Step-14: Creating final outputs in a timeseries format and cleaning up.")
         else:
-            rat_logger.info("Finished Step-13: Calculation of Outflow, Evaporation and Storage change")
-            ##---------- Mass-balance Approach ends and then post-processed outputs to obtain timeseries  -----------------##
-        
-        # Clearing out memory space as per user input 
-        if(config['CLEAN_UP']['clean_metsim']):
-            rat_logger.info("Clearing up memory space: Removal of metsim output files")
-            cleaner.clean_metsim()
-        if(config['CLEAN_UP']['clean_vic']):
-            rat_logger.info("Clearing up memory space: Removal of vic input, output files and previous init_state_files")
-            cleaner.clean_vic()
-        if(config['CLEAN_UP']['clean_routing']):
-            rat_logger.info("Clearing up memory space: Removal of routing input and output files")
-            cleaner.clean_routing()
-        if(config['CLEAN_UP']['clean_gee']):
-            rat_logger.info("Clearing up memory space: Removal of unwanted gee extracted small chunk files")
-            cleaner.clean_gee()
-        if(config['CLEAN_UP']['clean_altimetry']):
-            rat_logger.info("Clearing up memory space: Removal of raw altimetry downloaded data files.")
-            cleaner.clean_altimetry()
-
-    close_logger()
\ No newline at end of file
+            rat_logger.info("Finished Step-14: Creating final outputs in a timeseries format and cleaning up.")
+            ##----------Converted all time-series to final output csv format and cleaned up----------## 
+
+    close_logger()
+    return (no_errors, latest_altimetry_cycle)
\ No newline at end of file
diff --git a/src/rat/run_rat.py b/src/rat/run_rat.py
index 87d9af42..758c15e6 100644
--- a/src/rat/run_rat.py
+++ b/src/rat/run_rat.py
@@ -1,22 +1,34 @@
 import pandas as pd
+import numpy as np
 import yaml
 import os
 import argparse
 import configparser
 import ee
+import ruamel_yaml as ryaml
+from pathlib import Path 
+import datetime
+import copy
 
 from rat.utils.logging import init_logger,close_logger
 import rat.ee_utils.ee_config as ee_configuration
+from rat.rat_basin import rat_basin
 
 #------------ Define Variables ------------#
-def run_rat(config_fn):
+def run_rat(config_fn, operational_latency=None):
     """Runs RAT as per configuration defined in `config_fn`.
 
     parameters:
         config_fn (str): Path to the configuration file
+        operational_latency (int): Number of days in the past from today to end RAT operational run . RAT won't run operationally if operational_latency is None. Default is None.
     """
-    config = yaml.safe_load(open(config_fn, 'r'))
 
+    # Reading config with comments
+    config_fn = Path(config_fn).resolve()
+    ryaml_client = ryaml.YAML()
+    config = ryaml_client.load(config_fn.read_text())
+
+    # Logging this run
     log_dir = os.path.join(config['GLOBAL']['data_dir'],'runs','logs','')
     log = init_logger(
         log_dir,
@@ -27,8 +39,8 @@ def run_rat(config_fn):
         logger_name='run_rat',
         for_basin=False
     )
-    # 
-    # Tring the ee credentials given by user
+     
+    # Trying the ee credentials given by user
     try:
         log.info("Checking earth engine credentials")
         secrets = configparser.ConfigParser()
@@ -41,23 +53,100 @@ def run_rat(config_fn):
     except:
         log.info("Failed to connect to Earth engine. Wrong credentials. If you want to use Surface Area Estimations from RAT, please update the EE credentials.")
 
-    from rat.rat_basin import rat
-
+    ############ ----------- Single basin run ---------------- ################
     if(not config['GLOBAL']['multiple_basin_run']):
         log.info('############## Starting RAT for '+config['BASIN']['basin_name']+' #################')
-        rat(config, log)
-        log.info('############## Succesfully run RAT for '+config['BASIN']['basin_name']+' #################')
+        
+        # Checking if Rat is running operationally with some latency. If yes, update start, end and vic_init_state dates.
+        if operational_latency:
+            try:
+                log.info(f'Running RAT operationally at a latency of {operational_latency}. Updating start and end date.')
+                config['BASIN']['start'] = copy.deepcopy(config['BASIN']['end'])
+                config['BASIN']['vic_init_state_date'] = copy.deepcopy(config['BASIN']['end'])
+                config['BASIN']['end'] = datetime.datetime.now().date() - datetime.timedelta(days=int(operational_latency))
+            except:
+                log.exception('Failed to update start and end date for RAT to run operationally. Please make sure RAT has been run atleast once before.')
+                log.info('polo')
+                return None
+        # Running RAT if start < end date
+        if (config['BASIN']['start'] >= config['BASIN']['end']):
+            log.error('Sorry, RAT operational run for '+config['BASIN']['basin_name']+' failed.')
+            log.error(f"Start date - {config['BASIN']['start']} is before the end date - {config['BASIN']['end']}")
+            return None
+        else:
+            ryaml_client.dump(config, config_fn.open('w'))
+            no_errors, latest_altimetry_cycle = rat_basin(config, log)
+        # Displaying and storing RAT function outputs
+        if(latest_altimetry_cycle):
+            config['ALTIMETER']['last_cycle_number'] = latest_altimetry_cycle
+            ryaml_client.dump(config, config_fn.open('w'))
+        if(no_errors>0):
+            log.info('############## RAT run finished for '+config['BASIN']['basin_name']+ 'with '+str(no_errors)+' errors. #################')
+        elif(no_errors==0):
+            log.info('############## Succesfully run RAT for '+config['BASIN']['basin_name']+' #################')
+        else:
+            log.error('############## RAT run failed for '+config['BASIN']['basin_name']+' #################')
+
+    ############ ----------- Multiple basin run ---------------- ################
     else:
         basins_to_process = config['GLOBAL']['basins_to_process']
-        basins_metadata = pd.read_csv(config['GLOBAL']['basins_metadata'],header=[0,1])
+        # For each basin
         for basin in basins_to_process:
             log.info('############## Starting RAT for '+basin+' #################')
+            # Reading basins metadata
+            basins_metadata = pd.read_csv(config['GLOBAL']['basins_metadata'],header=[0,1])
+            # Checking if Rat is running operationally with some latency. If yes, update start, end and vic_init_state dates.
+            if operational_latency:
+                try:
+                    log.info(f'Running RAT operationally at a latency of {operational_latency}. Updating start and end date.')
+                    if ('BASIN','end') not in basins_metadata.columns:
+                        config['BASIN']['start'] = copy.deepcopy(config['BASIN']['end'])
+                        config['BASIN']['vic_init_state_date'] = copy.deepcopy(config['BASIN']['end'])
+                        config['BASIN']['end'] = datetime.datetime.now().date() - datetime.timedelta(days=int(operational_latency))
+                    else:
+                        # Updating start
+                        if ('BASIN','start') not in basins_metadata.columns:
+                            basins_metadata['BASIN','start'] = None    
+                        basins_metadata['BASIN','start'].where(basins_metadata['BASIN','basin_name']!= basin, basins_metadata['BASIN','end'], inplace=True)
+                        # Updating end
+                        operational_end_date = datetime.datetime.now().date() - datetime.timedelta(days=int(operational_latency))
+                        basins_metadata['BASIN','end'].where(basins_metadata['BASIN','basin_name']!= basin, operational_end_date, inplace=True)
+                        # Updating vic_init_state_date
+                        if ('BASIN','vic_init_state_date') not in basins_metadata.columns:
+                            basins_metadata['BASIN','vic_init_state_date'] = None    
+                        basins_metadata['BASIN','vic_init_state_date'].where(basins_metadata['BASIN','basin_name']!= basin, basins_metadata['BASIN','start'], inplace=True)
+                except:
+                    log.exception('Failed to update start and end date for RAT to run operationally. Please make sure RAT has been run atleast once before.')
+                    return None
+            # Extracting basin information and populating it in config if it's not NaN
             basin_info = basins_metadata[basins_metadata['BASIN']['basin_name']==basin]
-            # Extracting basin information and populating it in config
+            config_copy = copy.deepcopy(config)
             for col in basin_info.columns:
-                config[col[0]][col[1]] = basin_info[col[0]][col[1]].values[0]
-            rat(config, log)
-            log.info('############## Succesfully run RAT for '+basin+' #################')
+                if(not np.isnan(basin_info[col[0]][col[1]].values[0])):
+                    config_copy[col[0]][col[1]] = basin_info[col[0]][col[1]].values[0]
+            # Running RAT if start < end date
+            if (config_copy['BASIN']['start'] >= config_copy['BASIN']['end']):
+                log.error('Sorry, RAT operational run for '+config_copy['BASIN']['basin_name']+' failed.')
+                log.error(f"Start date - {config_copy['BASIN']['start']} is before the end date - {config_copy['BASIN']['end']}")
+                continue
+            else:
+                basins_metadata.to_csv(config['GLOBAL']['basins_metadata'], index=False)
+                ryaml_client.dump(config, config_fn.open('w'))
+                no_errors, latest_altimetry_cycle = rat_basin(config_copy, log)
+            # Displaying and storing RAT function outputs
+            if(latest_altimetry_cycle):
+                # If column doesn't exist in basins_metadata, create one
+                if ('ALTIMETER','last_cycle_number') not in basins_metadata.columns:
+                    basins_metadata['ALTIMETER','last_cycle_number'] = None    
+                basins_metadata['ALTIMETER','last_cycle_number'].where(basins_metadata['BASIN','basin_name']!= basin, latest_altimetry_cycle, inplace=True)
+                basins_metadata.to_csv(config_copy['GLOBAL']['basins_metadata'], index=False)
+            if(no_errors>0):
+                log.info('############## RAT run finished for '+config_copy['BASIN']['basin_name']+ 'with '+str(no_errors)+' errors. #################')
+            elif(no_errors==0):
+                log.info('############## Succesfully run RAT for '+config_copy['BASIN']['basin_name']+' #################')
+            else:
+                log.error('############## RAT run failed for '+config_copy['BASIN']['basin_name']+' #################')
+
     # Clsoing logger
     close_logger('rat_run')
 
diff --git a/src/rat/utils/clean.py b/src/rat/utils/clean.py
index 3add5960..5a20c630 100644
--- a/src/rat/utils/clean.py
+++ b/src/rat/utils/clean.py
@@ -37,7 +37,7 @@ def clean_vic(self):
         
         try:
             vic_init_states_dir_path = os.path.join(self.basin_data_dir,'vic','vic_init_states','') 
-            days_old = 2 #n max of days
+            days_old = 15 #n max of days
 
             time_interval = datetime.now() - timedelta(days_old)
             file_namelist = os.listdir(vic_init_states_dir_path)
@@ -46,7 +46,7 @@ def clean_vic(self):
             for file in file_namelist:
                 path = os.path.join(vic_init_states_dir_path, file)
                 filetime = datetime.fromtimestamp(os.path.getctime(path))
-                if filetime > time_interval:
+                if filetime < time_interval:
                     os.remove(path)
 
         except:
@@ -67,7 +67,7 @@ def clean_routing(self):
 
         try:
             rout_init_states_dir_path = os.path.join(self.basin_data_dir,'ro','rout_state_file','') 
-            days_old = 2 #n max of days
+            days_old = 15 #n max of days
 
             time_interval = datetime.now() - timedelta(days_old)
             file_namelist = os.listdir(rout_init_states_dir_path)
@@ -76,7 +76,7 @@ def clean_routing(self):
             for file in file_namelist:
                 path = os.path.join(rout_init_states_dir_path, file)
                 filetime = datetime.fromtimestamp(os.path.getctime(path))
-                if filetime > time_interval:
+                if filetime < time_interval:
                     os.remove(path)
 
         except:
@@ -134,7 +134,7 @@ def clean_previous_outputs(self):
             print("No previous altimetry_extracted folder to delete")
         
         try:
-            web_format_data_path = os.path.join(self.basin_data_dir,'web_format_data')
-            shutil.rmtree(web_format_data_path)
+            final_outputs_path = os.path.join(self.basin_data_dir,'final_outputs')
+            shutil.rmtree(final_outputs_path)
         except:
-            print("No web_format_data folder to delete with previous outputs")
+            print("No final_outputs folder to delete with previous outputs")
diff --git a/src/rat/utils/convert_for_website.py b/src/rat/utils/convert_to_final_outputs.py
similarity index 99%
rename from src/rat/utils/convert_for_website.py
rename to src/rat/utils/convert_to_final_outputs.py
index 2054ddc2..8cc90570 100644
--- a/src/rat/utils/convert_for_website.py
+++ b/src/rat/utils/convert_to_final_outputs.py
@@ -7,7 +7,7 @@
 def convert_sarea(sarea_dir, website_v_dir):
     # Surface Area
     sarea_paths = [os.path.join(sarea_dir, f) for f in os.listdir(sarea_dir) if f.endswith(".csv")]
-    sarea_web_dir = create_directory(os.path.join(website_v_dir,'sarea' ),True)
+    sarea_web_dir = create_directory(os.path.join(website_v_dir,'sarea_tmsos' ),True)
     for sarea_path in sarea_paths:
         res_name = os.path.splitext(os.path.split(sarea_path)[-1])[0]
 
diff --git a/src/rat/utils/files_creator.py b/src/rat/utils/files_creator.py
index 9c98ca99..d119e074 100644
--- a/src/rat/utils/files_creator.py
+++ b/src/rat/utils/files_creator.py
@@ -138,5 +138,7 @@ def create_basin_reservoir_shpfile(reservoir_shpfile,reservoir_shpfile_column_di
         reservoirs_gdf = reservoirs.merge(stations_df[reservoirs_gdf_column_dict['dam_name_column']], 
                                         how='inner', on=reservoirs_gdf_column_dict['dam_name_column'])
     
-    reservoirs_gdf.to_file(savepath)
+    if(reservoirs_gdf.empty):
+        raise Exception('Reservoir names in reservoir shapefile are not matching with the station names in the station file used for routing.')
+    reservoirs_gdf.to_file(savepath, index=False)
 
diff --git a/src/rat/utils/utils.py b/src/rat/utils/utils.py
index 98a54c71..1c676bc7 100644
--- a/src/rat/utils/utils.py
+++ b/src/rat/utils/utils.py
@@ -45,6 +45,9 @@ def clip_ts(*tss, which='left'):
     """
     mint = max([min(ts.index) for ts in tss])
     maxt = min([max(ts.index) for ts in tss])
+    
+    if mint > maxt:
+        raise Exception('No overlapping time period between the time series.')
 
     if which == 'both':
         clipped_tss = [ts.loc[(ts.index>=mint)&(ts.index<=maxt)] for ts in tss]
diff --git a/src/rat/utils/vic_param_reader.py b/src/rat/utils/vic_param_reader.py
index c3abe6f7..bd17afea 100644
--- a/src/rat/utils/vic_param_reader.py
+++ b/src/rat/utils/vic_param_reader.py
@@ -108,16 +108,18 @@ def __init__(self, config, basin_name, startdate=None, enddate=None, vic_output_
         self.straight_from_metsim = False
         self.save_init_state = save_init_state
         #VIC State Save Date
-        self.vic_init_state_save_date = config['BASIN']['end']+datetime.timedelta(days=1)
+        self.vic_init_state_save_date = config['BASIN']['end']
 
         if runname is None:
             self.runname = datetime.datetime.now().strftime("%Y%m%d_%H%M%S")
         else:
             self.runname = str(runname)
-        if self.workspace:
-            self.workspace = create_directory(os.path.join(self.workspace, f'run_{self.runname}'),True)
-        else :
-            self.workspace = create_directory(os.path.join(config['GLOBAL']['data_dir'],config['BASIN']['region_name'],'basins',basin_name,'vic','vic_workspace',f'run_{self.runname}'),True)
+        
+        if (self.intermediate_files):
+            if self.workspace:
+                self.workspace = create_directory(os.path.join(self.workspace, f'run_{self.runname}'),True)
+            else :
+                self.workspace = create_directory(os.path.join(config['GLOBAL']['data_dir'],config['BASIN']['region_name'],'basins',basin_name,'vic','vic_workspace',f'run_{self.runname}'),True)
 
         if self.init_param_file:
             self._load_from_vic_param()


Test patch:
diff --git a/.github/workflows/run_tests.yaml b/.github/workflows/run_tests.yaml
index d16a33b3..a7a83f93 100644
--- a/.github/workflows/run_tests.yaml
+++ b/.github/workflows/run_tests.yaml
@@ -4,7 +4,6 @@ on: [push]
 
 jobs:
   build:
-
     runs-on: ubuntu-latest
     strategy:
       matrix:
@@ -15,24 +14,43 @@ jobs:
         uses: actions/setup-python@v4
         with:
           python-version: ${{ matrix.python-version }}
-      - name: Setup environment with mamba
+      - name: Setup environment with conda
         uses: conda-incubator/setup-miniconda@v2
         with:
-          mamba-version: "*"
           auto-update-conda: true
-          channels: conda-forge,defaults
-          environment-file: environment.yml
-      - name: Install test modules
-        shell: bash -l {0}
-        run: |
-          conda info
-          conda install pytest conda-build -y
-      - name: Build RAT
+          python-version: ${{ matrix.python-version }}
+          environment-file: test_env.yml
+          activate-environment: rat_test_env
+      # - name: Test conda
+      #   shell: bash -l {0}
+      #   run: |
+      #     which conda
+      #     conda info
+      #     conda env list
+      #     conda env export --from-history
+      # - name: Build and install RAT
+      #   shell: bash -l {0}
+      #   run: |
+      #     conda build .
+      #     conda install -c file://${CONDA_PREFIX}/conda-bld/ rat
+      # - name: Test RAT
+      #   shell: bash -l {0}
+      #   run: |
+      #     mkdir ~/rat_test
+      #     rat init -d ~/rat_test
+      #     rat test -d ~/rat_test
+      - name: Install RAT (develop mode)
         shell: bash -l {0}
         run: |
+          mamba env update --file environment.yml
           conda develop src/
       - name: Test RAT
         shell: bash -l {0}
         run: |
-          conda activate test
-          pytest
+          mkdir ~/rat_test
+          which conda
+          which python
+          conda list
+          mamba install ruamel.yaml --force-reinstall
+          python src/rat/cli/rat_cli.py init -d ~/rat_test -s ${{ secrets.SECRETS }}
+          python src/rat/cli/rat_cli.py test -d ~/rat_test -s ${{ secrets.SECRETS }}
diff --git a/src/rat/cli/rat_test_config.py b/src/rat/cli/rat_test_config.py
new file mode 100644
index 00000000..1f47dc26
--- /dev/null
+++ b/src/rat/cli/rat_test_config.py
@@ -0,0 +1,151 @@
+from datetime import datetime
+
+DOWNLOAD_LINK = {
+    'test_data': "https://www.dropbox.com/s/huen14qwzphya61/test_data.zip?dl=1"
+}
+
+PATHS = {
+    'GUNNISON' : {
+        'GLOBAL': {
+            'data_dir': 'data/test_output',
+            'basin_shpfile': 'data/test_data/gunnison/gunnison_boundary/gunnison_boundary.shp',
+        },
+
+        'METSIM':{
+            'metsim_env': 'models/metsim',
+            'metsim_param_file': 'params/metsim/params.yaml',
+            'metsim_domain_file': 'data/test_data/gunnison/metsim_inputs/domain.nc'
+        },
+
+        'VIC': {
+            'vic_env': 'models/vic',
+            'vic_param_file': 'params/vic/vic_params.txt',
+            'vic_soil_param_file': 'data/test_data/gunnison/vic_basin_params/vic_soil_param.nc',
+            'vic_domain_file': 'data/test_data/gunnison/vic_basin_params/vic_domain.nc'
+        },
+
+        'ROUTING': {
+            'route_model': 'models/routing/rout',
+            'route_param_file': 'params/routing/route_param.txt',
+            'station_latlon_path': 'data/test_data/gunnison/gunnison_reservoirs/gunnison_reservoirs_locations.csv'
+        },
+
+        'ROUTING PARAMETERS':{
+            'flow_direction_file': 'data/test_data/gunnison/fl/fl.asc',
+            'uh': 'params/routing/uh.txt'
+        },
+
+        'GEE': {
+            'reservoir_vector_file': 'data/test_data/gunnison/gunnison_reservoirs/gunnison_reservoirs_named.geojson',
+        },
+
+        'ALTIMETER': {
+            'altimeter_tracks': 'global_data/global_altimetry/j3_tracks.geojson',
+            'geoid_grid': 'global_data/global_altimetry/geoidegm2008grid.mat'
+        },
+    },
+    'NUECES':{
+        'GLOBAL': {
+            'data_dir': 'data/test_output',
+            'basin_shpfile': 'data/test_data/Nueces/basin_shapefile/Nueces_polygon.json',
+        },
+
+        'METSIM':{
+            'metsim_env': 'models/metsim',
+            'metsim_param_file': 'params/metsim/params.yaml',
+            'metsim_domain_file': 'data/test_data/Nueces/metsim_inputs/domain.nc'
+        },
+
+        'VIC': {
+            'vic_env': 'models/vic',
+            'vic_param_file': 'params/vic/vic_params.txt',
+            'vic_soil_param_file': 'data/test_data/Nueces/vic_basin_params/vic_soil_param.nc',
+            'vic_domain_file': 'data/test_data/Nueces/vic_basin_params/vic_domain.nc'
+        },
+
+        'ROUTING': {
+            'route_model': 'models/routing/rout',
+            'route_param_file': 'params/routing/route_param.txt',
+            'station_latlon_path': 'data/test_data/Nueces/reservoirs/basin_station_latlon.csv'
+        },
+
+        'ROUTING PARAMETERS':{
+            'flow_direction_file': 'data/test_data/Nueces/fl/fl.asc',
+            'uh': 'params/routing/uh.txt'
+        },
+
+        'GEE': {
+            'reservoir_vector_file': 'data/test_data/Nueces/reservoirs/basin_reservoirs.shp',
+        },
+
+        'ALTIMETER': {
+            'altimeter_tracks': 'global_data/global_altimetry/j3_tracks.geojson',
+            'geoid_grid': 'global_data/global_altimetry/geoidegm2008grid.mat'
+        },
+        }
+}
+
+PARAMS = {
+    'GUNNISON':{
+        'GLOBAL': {
+            'steps':[1,2,3,4,5,6,7,8,9,10,12,13,14],
+            'basin_shpfile_column_dict': {'id': 'gridcode'},
+            'multiple_basin_run': False,
+        },
+        
+        'BASIN': {
+            'region_name': 'colorado',
+            'basin_name': 'gunnison',
+            'basin_id': 0,
+            'spin_up': False,
+            'start': datetime(2022, 1, 1),
+            'end': datetime(2022, 1, 31),
+        },
+
+        'VIC': {
+            'vic_global_data': False,
+        },
+
+        'GEE': {
+            'reservoir_vector_file_columns_dict': {
+                'id_column': None, 
+                'dam_name_column': 'DAM_NAME',
+                'area_column': 'area'}
+        },
+
+        'ROUTING': {
+            'station_global_data': False
+        }
+    },
+    'NUECES':{
+        'GLOBAL': {
+            'steps':[1,2,3,4,5,6,7,8,9,10,12,13,14],
+            'basin_shpfile_column_dict': {'id': 'MRBID'},
+            'multiple_basin_run': False,
+        },
+        
+        'BASIN': {
+            'region_name': 'Texas',
+            'basin_name': 'Nueces',
+            'basin_id': 4223,
+            'spin_up': False,
+            'start': datetime(2022, 8, 1),
+            'end': datetime(2022, 8, 31),
+        },
+
+        'VIC': {
+            'vic_global_data': False,
+        },
+
+        'GEE': {
+            'reservoir_vector_file_columns_dict': 
+                                    {'id_column' : 'GRAND_ID',
+                                    'dam_name_column' : 'uniq_id',
+                                    'area_column'     : 'AREA_SKM'} 
+        },
+
+        'ROUTING': {
+            'station_global_data': False
+        }
+    }
+}
\ No newline at end of file
diff --git a/test_env.yml b/test_env.yml
new file mode 100644
index 00000000..07e38b8f
--- /dev/null
+++ b/test_env.yml
@@ -0,0 +1,7 @@
+name: rat_test_env
+channels:
+  - conda-forge
+  - defaults
+dependencies:
+  - conda-build
+  - mamba
\ No newline at end of file
diff --git a/tests/test_imports.py b/tests/test_imports.py
index 098a620a..65e130ad 100644
--- a/tests/test_imports.py
+++ b/tests/test_imports.py
@@ -1,5 +1,4 @@
 import pytest
 
 def test_imports():
-    from rat.utils.logging import init_logger,close_logger
-    # from rat.rat_basin import rat
\ No newline at end of file
+    import rat
\ No newline at end of file