From de69cd56a8784d0ad73ff85215ee4e3d33db8bdd Mon Sep 17 00:00:00 2001
From: cuill <cui86ll@gmail.com>
Date: Thu, 21 Sep 2023 16:32:34 -0400
Subject: [PATCH 01/15] added link to online docs

---
 README.md | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/README.md b/README.md
index 0052f486..7eb39c04 100644
--- a/README.md
+++ b/README.md
@@ -59,6 +59,10 @@ hgrid = Hgrid.open('hgrid.gr3')
 hgrid.write("/path/to/output/file.2dm", fmt='2dm')
 ```
 ---
+## Online manual
+
+https://schism-dev.github.io/schism/master/getting-started/pre-processing-with-pyschism/overview.html
+
 ## References
 
 If you used this software as part of your work, please use the following citation format.

From cebb988d9c5692c4a4ceaafbfcaa0bd3a01afc55 Mon Sep 17 00:00:00 2001
From: cuill <cui86ll@gmail.com>
Date: Thu, 21 Sep 2023 16:33:04 -0400
Subject: [PATCH 02/15] added package pylib-essentials

---
 setup.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/setup.py b/setup.py
index 1e301d49..be27b4ac 100755
--- a/setup.py
+++ b/setup.py
@@ -130,6 +130,7 @@ def run(self):
         'xmltodict',
         'zarr',
         'appdirs',
+        'pylib-essentials',
     ],
     entry_points={'console_scripts': ['pyschism = pyschism.__main__:main']},
     tests_require=['nose'],

From 45fd49a032ad60acec8809fe3fd94e9d05a2b6ce Mon Sep 17 00:00:00 2001
From: cuill <cui86ll@gmail.com>
Date: Thu, 21 Sep 2023 16:35:54 -0400
Subject: [PATCH 03/15] removed unused scripts and update an example script

---
 examples/Bctides/gen_bctides.py     |  25 +++---
 pyschism/forcing/bctides/bctides.py | 108 -----------------------
 pyschism/forcing/bctides/elev2d.py  |  74 ----------------
 pyschism/forcing/bctides/mod3d.py   | 130 ----------------------------
 pyschism/forcing/bctides/nudge.py   | 111 ------------------------
 pyschism/forcing/bctides/tides.py   |   4 -
 pyschism/forcing/bctides/uv3d.py    |  83 ------------------
 tests/test_bctides.py               |   1 -
 8 files changed, 14 insertions(+), 522 deletions(-)
 delete mode 100644 pyschism/forcing/bctides/elev2d.py
 delete mode 100644 pyschism/forcing/bctides/mod3d.py
 delete mode 100644 pyschism/forcing/bctides/nudge.py
 delete mode 100644 pyschism/forcing/bctides/uv3d.py

diff --git a/examples/Bctides/gen_bctides.py b/examples/Bctides/gen_bctides.py
index 39097a3c..c735169c 100644
--- a/examples/Bctides/gen_bctides.py
+++ b/examples/Bctides/gen_bctides.py
@@ -29,34 +29,37 @@
     logging.getLogger("pyschism").setLevel(logging.DEBUG)
 
     start_date = datetime(2017, 12, 1)
-    rnday = 396
-    end_date = start_date + timedelta(days=61)
+    rnday = 10
     outdir = './'
     hgrid = Hgrid.open("./hgrid.gr3", crs="epsg:4326")
 
     #elevation
+    #constituents: major [M2, S2, N2, K2, K1, O1, P1, Q1]
     iet3 = iettype.Iettype3(constituents='major', database='fes2014')
     iet4 = iettype.Iettype4()
     iet5 = iettype.Iettype5(iettype3=iet3, iettype4=iet4)
 
     #velocity
-    ifl1 = ifltype.Ifltype1()
     ifl3 = ifltype.Ifltype3(constituents='major', database='fes2014')
     ifl4 = ifltype.Ifltype4()
     ifl5 = ifltype.Ifltype5(ifltype3=ifl3, ifltype4=ifl4)
 
-    #salinity
-    isa4 = isatype.Isatype4()
-    isa2 = isatype.Isatype2(5, 1)
+    #salinity 
+    isa4 = isatype.Isatype4() #time-varying salinity
+    isa2 = isatype.Isatype2(0, 1) # constant salinity (=0)
 
     #temperature
-    ite4 = itetype.Itetype4()
-    ite2 = itetype.Itetype2(0, 1)
+    ite4 = itetype.Itetype4() #time-varying temperature
+    ite2 = itetype.Itetype2(0, 1) #constant temperature (=0)
 
+    #example1 - only one ocean open boundary and type is [5, 5, 4, 4]
+    #bctides = Bctides(hgrid, iettype=iet5, ifltype=ifl5, isatype=isa4, itetype=ite4)
+
+    #example2 - two ocean open boundaries, one river boundary and type is [[5, 5, 4, 4], [5, 5, 4, 4], [0, 0, 3, 3]]
     bctides=Bctides(
         hgrid, 
         iettype={'1': iet5, '2': iet5}, 
-        ifltype={'1': ifl5, '2': ifl5, '3': ifl1}, 
+        ifltype={'1': ifl5, '2': ifl5}, 
         isatype={'1': isa4, '2':isa4, '3': isa2}, 
         itetype={'1': ite4, '2':ite4, '3': ite2}
     )
@@ -64,7 +67,7 @@
     bctides.write(
         outdir, 
         start_date=start_date, 
-        end_date=end_date, 
+        end_date=rnday, 
         bctides=True, 
-        overwrite=True
+        overwrite=True,
     )
diff --git a/pyschism/forcing/bctides/bctides.py b/pyschism/forcing/bctides/bctides.py
index c92c2c93..56834e2e 100644
--- a/pyschism/forcing/bctides/bctides.py
+++ b/pyschism/forcing/bctides/bctides.py
@@ -10,9 +10,6 @@
 from pyschism import dates
 from pyschism.mesh.vgrid import Vgrid
 from pyschism.forcing.bctides import iettype, ifltype, isatype, itetype, itrtype, Tides
-from pyschism.forcing.bctides.elev2d import Elev2D
-from pyschism.forcing.bctides.uv3d import UV3D
-from pyschism.forcing.bctides.mod3d import TEM_3D, SAL_3D
 
 logger = logging.getLogger(__name__)
 
@@ -130,13 +127,7 @@ def write(
         start_date: datetime = None,
         end_date: Union[datetime, timedelta] = None,
         bctides: Union[bool, str] = True,
-        elev2D: Union[bool, str] = False,
-        uv3D: Union[bool, str] = False,
-        tem3D: Union[bool, str] = False,
-        sal3D: Union[bool, str] = False,
         overwrite: bool = False,
-        parallel_download=False,
-        progress_bar=True,
     ):
         if start_date is not None:
             self.start_date = start_date
@@ -150,105 +141,6 @@ def write(
             raise IOError("path exists and overwrite is False")
         with open(bctides, "w") as f:
             f.write(str(self))
-        # write nudge
-        for bctype, tracer in {"itetype": "TEM", "isatype": "SAL"}.items():
-            for boundary in self.gdf.itertuples():
-                data_source = getattr(boundary, bctype)
-                if data_source is not None:
-                    import importlib
-                    if hasattr(data_source, 'rlmax'):
-                        # This generates all the nudges and writes the nudge files.
-                        nudgemod = importlib.import_module('pyschism.forcing.bctides.nudge')
-                        nudgeclass = getattr(nudgemod, f'{tracer}_Nudge')
-                        _tracerfile = locals()[f'{tracer.lower()}3D']
-                        if _tracerfile is False:
-                            continue
-                        elif _tracerfile is True:
-                            _tracerfile = output_directory / f'{tracer}_nudge.gr3'
-                        nudgeclass(self, data_source, rlmax=data_source.rlmax, rnu_day=data_source.rnu_day
-                                   ).write(_tracerfile, overwrite=overwrite)
-                        break
-
-        def write_elev2D():
-            _elev2D = output_directory / "elev2D.th.nc" if elev2D is True else elev2D
-            Elev2D(self).write(
-                _elev2D,
-                self.start_date,
-                self.rnday,
-                timedelta(days=1),
-                overwrite,
-                progress_bar=progress_bar,
-            )
-
-        def write_uv3D():
-            # write uv3D.th.nc
-            _uv3D = output_directory / "uv3D.th.nc" if uv3D is True else uv3D
-            UV3D(self).write(
-                _uv3D,
-                self.start_date,
-                self.rnday,
-                timedelta(days=1),
-                overwrite,
-                progress_bar=progress_bar,
-            )
-
-        def write_tem3D():
-            # write TEM_3D.th.nc
-            _tem3D = output_directory / "TEM_3D.th.nc" if tem3D is True else tem3D
-            TEM_3D(self).write(
-                _tem3D,
-                self.start_date,
-                self.rnday,
-                timedelta(days=1),
-                overwrite,
-                progress_bar=progress_bar,
-            )
-
-        def write_sal3D():
-            _sal3D = output_directory / "SAL_3D.th.nc" if sal3D is True else sal3D
-            SAL_3D(self).write(
-                _sal3D,
-                self.start_date,
-                self.rnday,
-                timedelta(days=1),
-                overwrite,
-                progress_bar=progress_bar,
-            )
-
-        if parallel_download is True:
-            from multiprocessing import Process
-
-            jobs = [
-                Process(target=fn)
-                for fn, fl in (
-                    (write_elev2D, elev2D),
-                    (write_uv3D, uv3D),
-                    (write_tem3D, tem3D),
-                    (write_sal3D, sal3D)
-                ) if fl
-            ]
-            for job in jobs:
-                job.start()
-            for job in jobs:
-                job.join()
-            if any(j.exitcode != 0 for j in jobs):
-                raise RuntimeError("Some parallel writer jobs failed!")
-        else:
-            if elev2D:
-                write_elev2D()
-            if uv3D:
-                write_uv3D()
-            if tem3D:
-                write_tem3D()
-            if sal3D:
-                write_sal3D()
-
-        # def write_tracer(tracer):
-        #     tracer.write()
-
-        # for tracer in [self.temperature, self.salinity, *self.tracers]:
-        #     if tracer is not None:
-        #         write_tracer(tracer)
 
     def get_forcing_string(self, boundary, global_constituents):
 
diff --git a/pyschism/forcing/bctides/elev2d.py b/pyschism/forcing/bctides/elev2d.py
deleted file mode 100644
index 8b6a7992..00000000
--- a/pyschism/forcing/bctides/elev2d.py
+++ /dev/null
@@ -1,74 +0,0 @@
-from datetime import timedelta
-import pathlib
-
-# import geopandas as gpd
-from netCDF4 import Dataset
-
-
-class Elev2D:
-
-    def __init__(self, bctides):
-        self.bctides = bctides
-
-    def write(
-            self,
-            elev2D,
-            start_date,
-            rnday,
-            output_interval=timedelta(days=1),
-            overwrite: bool = False,
-            progress_bar: bool = True,
-    ):
-        elev2D = pathlib.Path(elev2D)
-        if elev2D.exists() and overwrite is not True:
-            raise IOError(f'File {elev2D} exists and overwrite is not True.')
-
-        # file_is_not_needed = True
-        timevec = None
-        for boundary in self.bctides.gdf.itertuples():
-            if boundary.iettype is not None:
-                if boundary.iettype.iettype in [4, 5]:
-                    ds = boundary.iettype.data_component
-                    datasets = ds.get_datasets(
-                                start_date,
-                                rnday,
-                                output_interval
-                                )
-                    timevec = range(len(datasets))
-
-        if timevec is None:
-            return
-
-        nOpenBndNodes = 0
-        for boundary in self.bctides.gdf.itertuples():
-            nOpenBndNodes += len(boundary.indexes)
-
-        dst = Dataset(elev2D, 'w', format='NETCDF4')
-        # dimensions
-        dst.createDimension('nOpenBndNodes', nOpenBndNodes)
-        dst.createDimension('one', 1)
-        dst.createDimension('time', None)
-        dst.createDimension('nLevels', 1)
-        dst.createDimension('nComponents', 1)
-
-        # variables
-        dst.createVariable('time_step', 'f', ('one',))
-        dst['time_step'][:] = int(output_interval.total_seconds())
-        dst.createVariable('time', 'f', ('time',))
-        dst['time'][:] = timevec
-        dst.createVariable('time_series', 'f',
-                           ('time', 'nOpenBndNodes', 'nLevels', 'nComponents'))
-        offset = 0
-        for boundary in self.bctides.gdf.itertuples():
-            if boundary.iettype is not None:
-                if boundary.iettype.iettype in [4, 5]:
-                    boundary.iettype.data_component.put_boundary_ncdata(
-                        boundary, dst, start_date, rnday, overwrite=overwrite,
-                        offset=offset, output_interval=output_interval,
-                        pixel_buffer=10, progress_bar=progress_bar)
-            else:
-                self.put_null_boundary_data(dst, len(boundary.indexes))
-            offset += len(boundary.indexes)
-
-    def put_null_boundary_data(self, dst, np):
-        raise NotImplementedError('Must write null data.')
diff --git a/pyschism/forcing/bctides/mod3d.py b/pyschism/forcing/bctides/mod3d.py
deleted file mode 100644
index dfb14452..00000000
--- a/pyschism/forcing/bctides/mod3d.py
+++ /dev/null
@@ -1,130 +0,0 @@
-from abc import ABC, abstractmethod
-from datetime import timedelta
-import pathlib
-
-from netCDF4 import Dataset
-
-
-class MOD_3D(ABC):
-
-    def __init__(self, bctides):
-        self.bctides = bctides
-
-    def write(
-            self,
-            path,
-            start_date,
-            rnday,
-            output_interval=timedelta(days=1),
-            overwrite: bool = False,
-            progress_bar: bool = True,
-    ):  
-        path = pathlib.Path(path)
-        if path.exists() and overwrite is not True:
-            raise IOError(f'File {path} exists and overwrite is not True.')
-
-        timevec = None
-        for boundary in self.bctides.gdf.itertuples():
-            obj = getattr(boundary, self.bctype)
-            if obj is not None:
-                bctype = getattr(obj, self.bctype)
-                if bctype == 4:
-                    datasets = obj.data_component.get_datasets(
-                                start_date,
-                                rnday,
-                                output_interval
-                                )
-                    timevec = range(len(datasets))
-
-        if timevec is None:
-            return
-
-        nOpenBndNodes = 0
-        for boundary in self.bctides.gdf.itertuples():
-            nOpenBndNodes += len(boundary.indexes)
-
-        dst = Dataset(path, 'w', format='NETCDF4')
-        # dimensions
-        dst.createDimension('nOpenBndNodes', nOpenBndNodes)
-        dst.createDimension('one', 1)
-        dst.createDimension('time', None)
-        dst.createDimension('nLevels', self.bctides.vgrid.nvrt)
-        dst.createDimension('nComponents', self.nComponents)
-
-        # variables
-        dst.createVariable('time_step', 'f', ('one',))
-        dst['time_step'][:] = int(output_interval.total_seconds())
-        dst.createVariable('time', 'f', ('time',))
-        dst['time'][:] = timevec
-        dst.createVariable('time_series', 'f',
-                           ('time', 'nOpenBndNodes', 'nLevels', 'nComponents'))
-        offset = 0
-        for boundary in self.bctides.gdf.itertuples():
-            obj = getattr(boundary, self.bctype)
-            if obj is not None:
-                bctype = getattr(obj, self.bctype)
-                if bctype == 4:
-                    obj.data_component.put_boundary_ncdata(
-                        self.bctides.hgrid,
-                        self.bctides.vgrid,
-                        boundary,
-                        dst,
-                        start_date,
-                        rnday,
-                        overwrite=overwrite,
-                        offset=offset,
-                        output_interval=output_interval,
-                        pixel_buffer=10,
-                        progress_bar=progress_bar
-                    )
-            else:
-                self.put_null_boundary_data(dst, len(boundary.indexes))
-            offset += len(boundary.indexes)
-
-    def put_null_boundary_data(self, dst, np):
-        raise NotImplementedError('Must write null data.')
-
-    @property
-    @abstractmethod
-    def nComponents(self):
-        pass
-
-    @property
-    @abstractmethod
-    def bctype(self):
-        pass
-
-    # @property
-    # @abstractmethod
-    # def name(self):
-    #     pass
-
-
-class TEM_3D(MOD_3D):
-
-    @property
-    def nComponents(self):
-        return 1
-
-    @property
-    def bctype(self):
-        return 'itetype'
-
-    # @property
-    # def name(self):
-    #     return 'temperature'
-
-
-class SAL_3D(MOD_3D):
-
-    @property
-    def nComponents(self):
-        return 1
-
-    @property
-    def bctype(self):
-        return 'isatype'
-
-    # @property
-    # def name(self):
-    #     return 'salinity'
diff --git a/pyschism/forcing/bctides/nudge.py b/pyschism/forcing/bctides/nudge.py
deleted file mode 100644
index 864b9c70..00000000
--- a/pyschism/forcing/bctides/nudge.py
+++ /dev/null
@@ -1,111 +0,0 @@
-from abc import abstractmethod
-import logging
-from time import time
-
-
-from numba import prange, jit
-import numpy as np
-
-from pyschism.mesh.gridgr3 import Gr3Field
-
-
-logger = logging.getLogger(__name__)
-
-
-class Nudge(Gr3Field):
-    def __init__(self, bctides, data_source, rlmax=1.5, rnu_day=0.25):
-
-        @jit(nopython=True, parallel=True)
-        def compute_nudge(lon, lat, opbd, out):
-
-            rnu_max = 1.0 / rnu_day / 86400.0
-
-            for idn in prange(lon.shape[0]):
-                if idn in opbd:
-                    rnu = rnu_max
-                    distmin = 0.0
-                else:
-                    distmin = np.finfo(np.float64).max
-                    for j in opbd:
-                        rl2 = np.sqrt(
-                            np.square(lon[idn] - lon[j]) + np.square(lat[idn] - lat[j])
-                        )
-                        if rl2 < distmin:
-                            distmin = rl2
-                rnu = 0.0
-                if distmin <= rlmax:
-                    rnu = (1 - distmin / rlmax) * rnu_max
-                out[idn] = rnu
-
-        out = np.zeros(bctides.hgrid.values.shape)
-        elnode = bctides.hgrid.elements.array
-        xy = bctides.hgrid.get_xy(crs="epsg:4326")
-        opbd = []
-        for boundary in bctides.gdf.itertuples():
-            forcing = getattr(boundary, self.bctype)
-            if forcing.nudge is True:
-                opbd.extend(list(boundary.indexes))
-        opbd = np.array(opbd)
-
-        if len(opbd) > 0:
-            logger.info(f"Begin compute_nudge for {self.name}.")
-            start = time()
-            compute_nudge(xy[:, 0], xy[:, 1], opbd, out)
-            
-            logger.info(f"Get the indexes for nudge.")
-            idxs_nudge = np.zeros(out.shape, dtype=int)
-            idxs = np.where(out > 0)[0]
-            idxs_nudge[idxs] = 1
-            idxs = np.where(np.max(out[elnode], axis=1) > 0)[0]
-            fp = elnode[idxs, -1] < 0
-            idxs_nudge[elnode[idxs[fp], :3]] = 1
-            idxs_nudge[elnode[idxs[~fp], :]] = 1
-            idxs = np.where(idxs_nudge == 1)[0]
-            include = idxs
-            logger.info(f'The shape of include is {len(include)}')
- 
-            logger.info(f"compute_nudge took {time()-start} seconds.")
-        super().__init__(
-            nodes={i: (coord, out[i]) for i, coord in enumerate(bctides.hgrid.coords)},
-            elements=bctides.hgrid.elements.elements,
-            description=f"{rlmax}, {rnu_day}",
-            crs=bctides.hgrid.crs,
-        )
-
-        self.data_source = data_source
-
-    #def write()
-
-    @property
-    @abstractmethod
-    def name(self):
-        """ """
-
-    @property
-    @abstractmethod
-    def bctype(self):
-        """ """
-
-
-class TEM_Nudge(Nudge):
-    @property
-    def name(self):
-        return "temperature"
-
-    @property
-    def bctype(self):
-        return "itetype"
-
-
-class SAL_Nudge(Nudge):
-    @property
-    def name(self):
-        return "salinity"
-
-    @property
-    def bctype(self):
-        return "isatype"
-
-
-TempNudge = TEM_Nudge
-SaltNudge = SAL_Nudge
diff --git a/pyschism/forcing/bctides/tides.py b/pyschism/forcing/bctides/tides.py
index 0fa07ac6..3692e65f 100644
--- a/pyschism/forcing/bctides/tides.py
+++ b/pyschism/forcing/bctides/tides.py
@@ -2,17 +2,13 @@
 from datetime import datetime, timedelta, timezone
 from enum import Enum
 import logging
-# import os
 from typing import Callable, Union
 
 import numpy as np
 
-# from pyschism.forcing.tides import bctypes
 from pyschism.forcing.bctides.tpxo import TPXO
 from pyschism.forcing.bctides.fes2014 import FES2014
 from pyschism.forcing.bctides.hamtide import HAMTIDE
-# from pyschism.forcing.baroclinic.gofs import GOFS
-# from pyschism.forcing.baroclinic.rtofs import RTOFS
 
 _logger = logging.getLogger(__name__)
 
diff --git a/pyschism/forcing/bctides/uv3d.py b/pyschism/forcing/bctides/uv3d.py
deleted file mode 100644
index 9124977a..00000000
--- a/pyschism/forcing/bctides/uv3d.py
+++ /dev/null
@@ -1,83 +0,0 @@
-from datetime import timedelta
-import pathlib
-
-# import geopandas as gpd
-from netCDF4 import Dataset
-
-
-class UV3D:
-
-    def __init__(self, bctides):
-        self.bctides = bctides
-
-    def write(
-            self,
-            uv3d,
-            start_date,
-            rnday,
-            output_interval=timedelta(days=1),
-            overwrite: bool = False,
-            progress_bar: bool = True,
-    ):
-        uv3d = pathlib.Path(uv3d)
-        if uv3d.exists() and overwrite is not True:
-            raise IOError(f'File {uv3d} exists and overwrite is not True.')
-
-        # file_is_not_needed = True
-        timevec = None
-        for boundary in self.bctides.gdf.itertuples():
-            if boundary.ifltype is not None:
-                if boundary.ifltype.ifltype in [4, 5]:
-                    ds = boundary.ifltype.data_component
-                    datasets = ds.get_datasets(
-                                start_date,
-                                rnday,
-                                output_interval
-                                )
-                    timevec = range(len(datasets))
-
-        if timevec is None:
-            return
-
-        nOpenBndNodes = 0
-        for boundary in self.bctides.gdf.itertuples():
-            nOpenBndNodes += len(boundary.indexes)
-
-        dst = Dataset(uv3d, 'w', format='NETCDF4')
-        # dimensions
-        dst.createDimension('nOpenBndNodes', nOpenBndNodes)
-        dst.createDimension('one', 1)
-        dst.createDimension('time', None)
-        dst.createDimension('nLevels', self.bctides.vgrid.nvrt)
-        dst.createDimension('nComponents', 2)
-
-        # variables
-        dst.createVariable('time_step', 'f', ('one',))
-        dst['time_step'][:] = int(output_interval.total_seconds())
-        dst.createVariable('time', 'f', ('time',))
-        dst['time'][:] = timevec
-        dst.createVariable('time_series', 'f',
-                           ('time', 'nOpenBndNodes', 'nLevels', 'nComponents'))
-        offset = 0
-        for boundary in self.bctides.gdf.itertuples():
-            if boundary.ifltype is not None:
-                if boundary.ifltype.ifltype in [4, 5]:
-                    boundary.ifltype.data_component.put_boundary_ncdata(
-                        self.bctides.hgrid,
-                        self.bctides.vgrid,
-                        boundary,
-                        dst,
-                        start_date,
-                        rnday,
-                        overwrite=overwrite,
-                        offset=offset,
-                        output_interval=output_interval,
-                        pixel_buffer=10,
-                        progress_bar=progress_bar
-                    )
-            else:
-                self.put_null_boundary_data(dst, len(boundary.indexes))
-            offset += len(boundary.indexes)
-
-    def put_null_boundary_data(self, dst, np):
-        raise NotImplementedError('Must write null data.')
diff --git a/tests/test_bctides.py b/tests/test_bctides.py
index 360d1330..50cc09f5 100755
--- a/tests/test_bctides.py
+++ b/tests/test_bctides.py
@@ -148,6 +148,5 @@ def test_simple_tidal_setup(self):
             tmpdir = pathlib.Path(dn)
             coldstart.write(tmpdir, overwrite=True)
 
-
 if __name__ == '__main__':
     unittest.main()

From 8734745a83d7241b93c825d6597d2ef059f59fa5 Mon Sep 17 00:00:00 2001
From: cuill <cui86ll@gmail.com>
Date: Wed, 27 Sep 2023 17:47:11 -0400
Subject: [PATCH 04/15] modified bctides per Joseph's request

---
 examples/Bctides/gen_bctides.py     |  56 +---
 pyschism/driver.py                  |   8 +-
 pyschism/forcing/bctides/bctides.py | 465 +++++++++++++---------------
 3 files changed, 225 insertions(+), 304 deletions(-)

diff --git a/examples/Bctides/gen_bctides.py b/examples/Bctides/gen_bctides.py
index c735169c..f005bb7f 100644
--- a/examples/Bctides/gen_bctides.py
+++ b/examples/Bctides/gen_bctides.py
@@ -3,12 +3,15 @@
 from datetime import datetime, timedelta
 import logging
 import pathlib
-import f90nml
 
 from pyschism.mesh import Hgrid
-from pyschism.forcing.bctides import Bctides, iettype, ifltype, isatype, itetype
-from pyschism.forcing.source_sink.nwm import NationalWaterModel, NWMElementPairings
-from pyschism.forcing.nws import NWS2, GFS, HRRR
+from pyschism.forcing.bctides import Bctides
+
+logging.basicConfig( 
+    format = "[%(asctime)s] %(name)s %(levelname)s: %(message)s",
+    force = True,
+)
+logging.getLogger("pyschism").setLevel(logging.DEBUG)
 
 if __name__ == "__main__":
 
@@ -21,53 +24,24 @@
     database = 'tpxo'
         ~/.local/share/tpxo/
     '''
-    #setup logging
-    logging.basicConfig(
-        format = "[%(asctime)s] %(name)s %(levelname)s: %(message)s",
-        force=True,
-    )
-    logging.getLogger("pyschism").setLevel(logging.DEBUG)
 
-    start_date = datetime(2017, 12, 1)
-    rnday = 10
+    start_date = datetime(2014, 12, 1)
+    rnday = 397
     outdir = './'
     hgrid = Hgrid.open("./hgrid.gr3", crs="epsg:4326")
 
-    #elevation
-    #constituents: major [M2, S2, N2, K2, K1, O1, P1, Q1]
-    iet3 = iettype.Iettype3(constituents='major', database='fes2014')
-    iet4 = iettype.Iettype4()
-    iet5 = iettype.Iettype5(iettype3=iet3, iettype4=iet4)
-
-    #velocity
-    ifl3 = ifltype.Ifltype3(constituents='major', database='fes2014')
-    ifl4 = ifltype.Ifltype4()
-    ifl5 = ifltype.Ifltype5(ifltype3=ifl3, ifltype4=ifl4)
-
-    #salinity 
-    isa4 = isatype.Isatype4() #time-varying salinity
-    isa2 = isatype.Isatype2(0, 1) # constant salinity (=0)
-
-    #temperature
-    ite4 = itetype.Itetype4() #time-varying temperature
-    ite2 = itetype.Itetype2(0, 1) #constant temperature (=0)
-
-    #example1 - only one ocean open boundary and type is [5, 5, 4, 4]
-    #bctides = Bctides(hgrid, iettype=iet5, ifltype=ifl5, isatype=isa4, itetype=ite4)
-
-    #example2 - two ocean open boundaries, one river boundary and type is [[5, 5, 4, 4], [5, 5, 4, 4], [0, 0, 3, 3]]
     bctides=Bctides(
         hgrid, 
-        iettype={'1': iet5, '2': iet5}, 
-        ifltype={'1': ifl5, '2': ifl5}, 
-        isatype={'1': isa4, '2':isa4, '3': isa2}, 
-        itetype={'1': ite4, '2':ite4, '3': ite2}
+        flags = [[5, 5, 4, 4], [5, 5, 4, 4], [0, 1, 2, 2]],
+        constituents = 'major',
+        database = 'fes2014',
+        tthconst = 10.0,
+        sthconst = 0.0,
     )
 
     bctides.write(
         outdir, 
         start_date=start_date, 
-        end_date=rnday, 
-        bctides=True, 
+        rnday=rnday, 
         overwrite=True,
     )
diff --git a/pyschism/driver.py b/pyschism/driver.py
index 77fc3430..5a4ed469 100644
--- a/pyschism/driver.py
+++ b/pyschism/driver.py
@@ -46,8 +46,6 @@ def write(
         driver: "ModelDriver",
         output_directory,
         overwrite,
-        parallel_download,
-        progress_bar,
     ):
 
         if self.bctides is not None:
@@ -56,8 +54,6 @@ def write(
                 start_date=driver.param.opt.start_date,
                 end_date=driver.param.core.rnday,
                 overwrite=True,
-                parallel_download=parallel_download,
-                progress_bar=progress_bar,
             )
 
         if self.nws is not None:
@@ -487,8 +483,6 @@ def write(
         salt_ic=True,
         # rtofs=True,
         hydrology=True,
-        parallel_download=True,
-        progress_bar=False,
     ):
         """Writes to disk the full set of input files necessary to run SCHISM."""
 
@@ -548,7 +542,7 @@ def obj_write(var, obj, default_filename, overwrite):
         obj_write(stations, self.stations, "station.in", overwrite)
 
         self.config.forcings.write(
-            self, output_directory, overwrite, parallel_download, progress_bar
+            self, output_directory, overwrite,
         )
 
         MakefileDriver(self.server_config, hotstart=self.hotstart).write(
diff --git a/pyschism/forcing/bctides/bctides.py b/pyschism/forcing/bctides/bctides.py
index 56834e2e..ac6b5a7e 100644
--- a/pyschism/forcing/bctides/bctides.py
+++ b/pyschism/forcing/bctides/bctides.py
@@ -1,306 +1,259 @@
 from datetime import datetime, timedelta
-from functools import cached_property, lru_cache
+from functools import cached_property, lru_cache 
 import pathlib
-from typing import Dict, Union
 import logging
 
-from ordered_set import OrderedSet
-import numpy as np
-
-from pyschism import dates
-from pyschism.mesh.vgrid import Vgrid
-from pyschism.forcing.bctides import iettype, ifltype, isatype, itetype, itrtype, Tides
+from pyschism.forcing.bctides import Tides
 
 logger = logging.getLogger(__name__)
-
-
-class IbctypeDescriptor:
-    def __init__(self, name, bctype):
-        self.name = name
-        self.bctype = bctype
-
-    def __get__(self, obj, val):
-        return obj.gdf[self.name]
-
-    def __set__(self, obj, val):
-        if val is not None:
-            if isinstance(val, dict):
-                for bnd_id, ibctype in val.items():
-                    if not isinstance(ibctype, (self.bctype, type(None))):
-                        raise TypeError(
-                            f"Argument {ibctype} for boundary {bnd_id} must be of type {self.bctype} "
-                            f" or None, not type {type(ibctype)}."
-                        )
-                    obj.gdf.at[np.where(obj.gdf["id"] == bnd_id)[0][0], self.bctype.__name__.lower()] = ibctype
-            else:
-                if not isinstance(val, self.bctype):
-                    raise TypeError(
-                        f"Argument {self.name} must be of type "
-                        f"{self.bctype}, not type {type(val)}."
-                    )
-                obj.gdf[self.name] = val
-
-
-class BctidesMeta(type):
-    def __new__(meta, name, bases, attrs):
-        bctypes = {
-            "iettype": iettype.Iettype,
-            "ifltype": ifltype.Ifltype,
-            "isatype": isatype.Isatype,
-            "itetype": itetype.Itetype,
-            "itrtype": itrtype.Itrtype,
-        }
-        for name, ibctype in bctypes.items():
-            attrs[name] = IbctypeDescriptor(name, ibctype)
-        return type(name, bases, attrs)
-
-
-class Bctides(metaclass=BctidesMeta):
-
-    start_date = dates.StartDate()
-    end_date = dates.EndDate()
-    run_days = dates.RunDays()
+class Bctides:
 
     def __init__(
         self,
         hgrid,
-        vgrid=None,
-        iettype: Union[Dict, iettype.Iettype] = None,
-        ifltype: Union[Dict, ifltype.Ifltype] = None,
-        isatype: Union[Dict, isatype.Isatype] = None,
-        itetype: Union[Dict, itetype.Itetype] = None,
-        itrtype: Union[Dict, itrtype.Itrtype] = None,
+        flags,
+        constituents = 'major',
+        database = 'tpxo',
         cutoff_depth: float = 50.0,
+        ethconst = None,
+        vthconst = None,
+        tthconst = None,
+        sthconst = None,
+        tobc = None,
+        sobc = None,
     ):
+
         self.hgrid = hgrid
-        self.vgrid = Vgrid.default() if vgrid is None else vgrid
         self.cutoff_depth = cutoff_depth
-        self.iettype = iettype
-        self.ifltype = ifltype
-        self.isatype = isatype
-        self.itetype = itetype
-        self.itrtype = itrtype
+        self.tides = Tides(constituents=constituents, tidal_database=database)
+        self.flags = flags
+        self.ethconst = ethconst 
+        self.vthconst = vthconst 
+        self.tthconst = tthconst 
+        self.sthconst = sthconst 
+        self.tobc = tobc if tobc is not None else 1.0
+        self.sobc = sobc if sobc is not None else 1.0
 
     def __str__(self):
+ 
+        #first line in the bctides.in is a note, not used in the schism code
         f = [
-            f"{str(self.start_date)}",
-            f"{self.ntip} {self.cutoff_depth}",
+            f"!{str(self.start_date)} UTC",
         ]
+
+        #get earth tidal potential and frequency
+        f.append(f"{self.ntip} {self.cutoff_depth} !number of earth tidal potential, cut-off depth for applying tidal potential")
         if self.ntip > 0:
-            for constituent in self.tides.get_active_potential_constituents():
+            for constituent in self.tides.get_active_potential_constituents(): 
                 forcing = self.tides(self.start_date, self.rnday, constituent)
                 f.append(
-                    " ".join(
-                        [
-                            f"{constituent}\n",
-                            f"{forcing[0]:G}",
-                            f"{forcing[1]:G}",
-                            f"{forcing[2]:G}",
-                            f"{forcing[3]:G}",
-                            f"{forcing[4]:G}",
-                        ]
-                    )
+                    " ".join([
+                        f"{constituent}\n",
+                        f"{forcing[0]:G}",
+                        f"{forcing[1]:G}",
+                        f"{forcing[2]:G}",
+                        f"{forcing[3]:G}",
+                        f"{forcing[4]:G}",
+                    ])
                 )
-        f.append(f"{self.nbfr:d}")
+
+        #get tidal boundary 
+        f.append(f"{self.nbfr:d} !nbfr")
         if self.nbfr > 0:
             for constituent in self.tides.get_active_forcing_constituents():
                 forcing = self.tides(self.start_date, self.rnday, constituent)
                 f.append(
-                    " ".join(
-                        [
-                            f"{constituent}\n",
-                            f"{forcing[2]:G}",
-                            f"{forcing[3]:G}",
-                            f"{forcing[4]:G}",
-                        ]
-                    )
+                    " ". join([
+                        f"{constituent}\n", 
+                        f"{forcing[2]:G}",
+                        f"{forcing[3]:G}",
+                        f"{forcing[4]:G}",
+
+                    ])
                 )
-        global_constituents = self.tides.get_active_constituents()
-        f.append(f"{len(self.gdf)}")
-        for boundary in self.gdf.itertuples():
-            f.append(self.get_forcing_string(boundary, global_constituents))
+ 
+        #get amplitude and phase for each open boundary         
+        f.append(f"{len(self.gdf)} !nope")
+        for ibnd, (boundary, flag) in enumerate(zip(self.gdf.itertuples(), self.flags)):
+            logger.info(f"Processing boundary {ibnd}:")
+            #number of nodes and flags
+            line = [
+                f"{len(boundary.indexes)}",
+                *[str(digit) for digit in flag],
+                f"!open bnd {ibnd+1}",
+            ]
+            f.append(" ".join(line))
+
+            #It only accounts for elevation, velocity, temperature, salinity
+            #TODO: add information for tracers
+            if len(flag) > 4:
+                raise NotImplementedError(f"Tracer module is not implemented yet!") 
+            iettype, ifltype, itetype, isatype = [i for i in flag]
+
+            #elevation boundary
+            logger.info(f"Elevation type: {iettype}")
+            if iettype == 1:
+                logger.warning(f'time history of elevation is read in from elev.th (ASCII)!')
+            elif iettype == 2:
+                logger.info("This boundary is forced by a constant elevation!")
+                if self.ethconst is None:
+                    raise IOError("You are choosing type 2 for elevation, please specify a value!")
+                else:
+                    f.append("{str(self.ethconst)} !T")
+            elif iettype == 4:
+                logger.warning('time history of elevation is read in from elev2D.th.nc (netcdf)')
+            elif iettype == 3 or iettype == 5:
+                if iettype == 5:
+                    logger.warning(f'Combination of 3 and 4, time history of elevation is read in from elev2D.th.nc!')
+                for constituent in self.tides.get_active_forcing_constituents():
+                    f.append(f"{constituent}")
+                    vertices = self.hgrid.get_xy(crs=self.hgrid.crs)[boundary.indexes, :]
+                    amp, phase = self.tides.get_elevation(constituent, vertices)
+                    for i in range(len(boundary.indexes)):
+                        f.append(f"{amp[i]: .6f} {phase[i]: .6f}")
+            elif iettype == 0:
+                logger.warning(f"elevations are not specified for this boundary (in this case the discharge must be specified)")
+            else:
+                raise IOError(f"Invalid type {iettype} for elevation!")
+
+            #velocity
+            logger.info(f"Velocity type: {ifltype}")
+            if ifltype == 0:
+                logger.info("Velocity is not sepcified, not input needed!")
+            elif ifltype == 1:
+                logger.warning("time history of discharge is read in from flux.th (ASCII)!")
+            elif ifltype == 2:
+                logger.info("This boundary is forced by a constant discharge!")
+                if self.vthconst is None:
+                    raise IOError("You are choosing type 2 for velocity, please specify dischage value!")
+                else:
+                    f.append("{str(self.vthconst)}")
+            elif ifltype == 3 or ifltype == 5:
+                if ifltype == 5:
+                    logger.warning(f'Combination of 3 and 4, time history of velocity is read in from uv.3D.th.nc!')
+                for constituent in self.tides.get_active_forcing_constituents():
+                    f.append(f"{constituent}")
+                    vertices = self.hgrid.get_xy(crs=self.hgrid.crs)[boundary.indexes, :]
+                    uamp, uphase, vamp, vphase = self.tides.get_velocity(constituent, vertices)
+                    for i in range(len(boundary.indexes)):
+                        f.append(f"{uamp[i]: .6f} {uphase[i]: .6f} {vamp[i]: .6f} {vphase[i]: 6f}")
+            elif ifltype == 4 or -4:
+                logger.warning("time history of velocity (not discharge!) is read in from uv3D.th.nc (netcdf)")
+                if ifltype == -4:
+                    if self.rel1 is None or self.rel2 is None:
+                        raise IOError(f"Please specify relaxation constants for inflow and outflow (between 0 and 1 with 1 being strongest nudging)")
+                    else:
+                        logger.info(f"relax value 1: {self.rel1}, relax value 2: {self.rel2}")
+            elif ifltype == -1: 
+                raise NotImplementedError(f"Velocity type {ifltype} not implemented yet!")
+                #logger.info(f"Flather type radiation b.c. (iettype must be 0 in this case)!")
+                #f.append(['vn_mean'])
+                
+                #TODO: add mean normal velocity at the node (at all levels)[
+            else:
+                raise IOError(f"Invalid type {ifltype} for veloicty!")
+
+            #temperature
+            logger.info(f"Temperature type: {itetype}")
+            if itetype == 0:
+                logger.warning("Temperature is not sepcified, not input needed!")
+            elif itetype == 1:
+                logger.warning("time history of temperature will be read in from TEM_1.th!")
+                logger.info(f"Using nudging factor: {self.tobc}!")
+                f.append(f"{self.tobc} !nudging factor for T")
+            elif itetype == 2:
+                logger.info("This boundary is forced by a constant temperature!")
+                if self.tthconst is None:
+                    raise IOError("You are choosing type 2 for temperature, please specify temperature value!")
+                else:
+                    f.append(f"{str(self.tthconst)} !T")
+
+                logger.info(f"Using nudging factor: {self.tobc}!")
+                f.append(f"{self.tobc} !nudging factor for T")
+            elif itetype == 3:
+                logger.info("Using initial temperature profile for inflow")
+                logger.info(f"Using nudging factor: {self.tobc}!")
+                f.append(f"{self.tobc} !nudging factor for T")
+            elif itetype == 4:
+                logger.warning("time history of temperature is read in from TEM_3D.th.nc (netcdf)!")
+                logger.info(f"Using nudging factor: {self.tobc}!")
+                f.append(f"{self.tobc} !nudging factor for T")
+            else:
+                raise IOError(f"Invalid type {itetype} for salinity!")        
+
+            #salinity
+            logger.info(f"Salinity type: {isatype}")
+            if isatype == 0:
+                logger.info("Salinity is not sepcified, not input needed!")
+            elif isatype == 1:
+                logger.warning("time history of salinity will be read in from SAL_1.th!")
+                logger.info(f"Using nudging factor: {self.sobc}!")
+                f.append(f"{self.sobc} !nudging factor for S")
+            elif isatype == 2:
+                logger.info("This boundary is forced by a constant salinity!")
+                if self.sthconst is None:
+                    raise IOError("You are choosing type 2 for salinity, please specify salinity value!")
+                else:
+                    f.append(f"{str(self.sthconst)} !S")
+
+                logger.info(f"Using nudging factor: {self.sobc}!")
+                f.append(f"{self.sobc} !nudging factor for S")
+            elif isatype == 3:
+                logger.info("Using initial salinity profile for inflow")
+                logger.info(f"Using nudging factor: {self.sobc}!")
+                f.append(f"{self.sobc} !nudging factor for S")
+            elif isatype == 4:
+                logger.warning("time history of salinity is read in from SAL_3D.th.nc (netcdf)!")
+                logger.info(f"Using nudging factor: {self.sobc}!")
+                f.append(f"{self.sobc} !nudging factor for S")
+            else:
+                raise IOError(f"Invalid type {isatype} for salinity!") 
+                
+
         return "\n".join(f)
 
+
     def write(
         self,
         output_directory,
         start_date: datetime = None,
-        end_date: Union[datetime, timedelta] = None,
-        bctides: Union[bool, str] = True,
-        overwrite: bool = False,
+        rnday = None,
+        constituents = ['K1', 'O1'],
+        overwrite: bool = True, 
     ):
+
         if start_date is not None:
             self.start_date = start_date
-        if end_date is not None:
-            self.run_days = end_date
-        # self.tidal_database.write(path, )
+        else:
+            raise IOError("Please specify the start_date!")
+
+        if rnday is not None:
+            self.rnday = rnday
+        else:
+            raise IOError("Please specify the number of days!")
+
+        if constituents is not None:
+            self.constituents = constituents
+        else:
+            raise IOError("Please specify tidal constituents!")
+
         output_directory = pathlib.Path(output_directory)
-        output_directory.mkdir(exist_ok=overwrite, parents=True)
-        bctides = output_directory / "bctides.in" if bctides is True else bctides
+        output_directory.mkdir(exist_ok=overwrite, parents=True) 
+        bctides = output_directory / "bctides.in"
         if bctides.exists() and not overwrite:
-            raise IOError("path exists and overwrite is False")
+            raise IOError("path exists and overwrite is False") 
         with open(bctides, "w") as f:
             f.write(str(self))
 
-    def get_forcing_string(self, boundary, global_constituents):
-
-        bctypes = [
-            boundary.iettype,
-            boundary.ifltype,
-            boundary.itetype,
-            boundary.isatype,
-        ]
-
-        def get_focing_digit(bctype):
-            if bctype is not None:
-                return bctype.forcing_digit
-            return "0"
-
-        line = [
-            f"{len(boundary.indexes)}",
-            *[str(digit) for digit in map(get_focing_digit, bctypes)],
-        ]
-        f = [" ".join(line)]
-        for bctype in bctypes:
-            if bctype is not None:
-                f.append(
-                    bctype.get_boundary_string(
-                        self.hgrid, boundary, global_constituents=global_constituents
-                    )
-                )
-        return "\n".join(f)
-
-    @property
-    def rnday(self):
-        return self.run_days
-
     @cached_property
     def gdf(self):
-        gdf = self.hgrid.boundaries.open.copy()
-        gdf["iettype"] = None
-        gdf["ifltype"] = None
-        gdf["isatype"] = None
-        gdf["itetype"] = None
-        gdf["itrtype"] = None
-        return gdf
-
+       return self.hgrid.boundaries.open.copy()
+        
     @property
     def ntip(self):
         return len(self.tides.get_active_potential_constituents())
 
     @property
     def nbfr(self):
-        return len(self.tides.get_active_forcing_constituents())
-
-    @cached_property
-    def tides(self):
-        return TidalConstituentCombiner(self.gdf)
-
-
-class TidalConstituentCombiner(Tides):
-    def __init__(self, gdf):
-        self.gdf = gdf
-        afc = self.get_active_forcing_constituents()
-        apc = self.get_active_potential_constituents()
-        for constituent in set([*afc, *apc]):
-            self.use_constituent(
-                constituent,
-                forcing=True if constituent in afc else False,
-                potential=True if constituent in apc else False,
-            )
-
-    def get_active_forcing_constituents(self):
-        active_constituents = OrderedSet()
-        for row in self.gdf.itertuples():
-            if row.iettype is not None:
-                if row.iettype.iettype in [3, 5]:
-                    [
-                        active_constituents.add(x)
-                        for x in row.iettype.tides.get_active_constituents()
-                    ]
-            if row.ifltype is not None:
-                if row.ifltype.ifltype in [3, 5]:
-                    [
-                        active_constituents.add(x)
-                        for x in row.ifltype.tides.get_active_constituents()
-                    ]
-
-        return list(active_constituents)
-
-    def get_active_potential_constituents(self):
-        active_constituents = OrderedSet()
-        for row in self.gdf.itertuples():
-            if row.iettype is not None:
-                if row.iettype.iettype in [3, 5]:
-                    [
-                        active_constituents.add(x)
-                        for x in row.iettype.tides.get_active_potential_constituents()
-                    ]
-            if row.ifltype is not None:
-                if row.ifltype.ifltype in [3, 5]:
-                    [
-                        active_constituents.add(x)
-                        for x in row.ifltype.tides.get_active_potential_constituents()
-                    ]
-
-        return list(active_constituents)
-
-    @lru_cache
-    def get_active_constituents(self):
-        return list(
-            OrderedSet(
-                [
-                    *self.get_active_potential_constituents(),
-                    *self.get_active_forcing_constituents(),
-                ]
-            )
-        )
-
-
-def ad_hoc_test():
-    from datetime import datetime
-    import logging
-
-    from pyschism.mesh import Hgrid
-    from pyschism.forcing.bctides import Bctides, iettype, ifltype, isatype, itetype
-
-    # setup logging
-    logging.basicConfig(
-        format="[%(asctime)s] %(name)s %(levelname)s: %(message)s",
-        force=True,
-    )
-    logging.getLogger("pyschism").setLevel(logging.DEBUG)
-
-    startdate = datetime(2018, 8, 17)
-    print(startdate)
-    rnday = 61
-    hgrid = Hgrid.open("./hgrid.gr3", crs="epsg:4326")
-
-    # Bctides
-    iet3 = iettype.Iettype3(constituents='major', database='tpxo')
-    iet4 = iettype.Iettype4()
-    iet5 = iettype.Iettype5(iettype3=iet3, iettype4=iet4)
-    ifl3 = ifltype.Ifltype3(constituents='major', database='tpxo')
-    ifl4 = ifltype.Ifltype4()
-    ifl5 = ifltype.Ifltype5(ifltype3=ifl3, ifltype4=ifl4)
-    isa3 = isatype.Isatype4()
-    # ite3 = itetype.Itetype4()
-    bctides = Bctides(hgrid, iettype={'1': iet5}, ifltype={'1': ifl5},
-                      isatype=isa3,
-                      itetype={'1': itetype.Itetype2(10, 1)}
-                      )
-    bctides.write(
-            './',
-            startdate,
-            rnday,
-            bctides=True,
-            elev2D=False,
-            uv3D=False,
-            tem3D=False,
-            sal3D=False,
-            overwrite=True
-            )
-
+        return len(self.tides.get_active_potential_constituents())
 
-if __name__ == "__main__":
-    ad_hoc_test()

From e205512a17af89a141c27bf414d15b4610d8bc81 Mon Sep 17 00:00:00 2001
From: cuill <cui86ll@gmail.com>
Date: Thu, 28 Sep 2023 13:13:36 -0400
Subject: [PATCH 05/15] modified gen_bctides.py to allow command-line inputs

---
 examples/Bctides/gen_bctides.py     | 49 +++++++++++++++++++++++------
 pyschism/forcing/bctides/bctides.py |  8 +++--
 2 files changed, 46 insertions(+), 11 deletions(-)

diff --git a/examples/Bctides/gen_bctides.py b/examples/Bctides/gen_bctides.py
index f005bb7f..db8e14ef 100644
--- a/examples/Bctides/gen_bctides.py
+++ b/examples/Bctides/gen_bctides.py
@@ -1,18 +1,23 @@
 from time import time
 import os
+import argparse
 from datetime import datetime, timedelta
 import logging
 import pathlib
+import json
 
 from pyschism.mesh import Hgrid
 from pyschism.forcing.bctides import Bctides
 
-logging.basicConfig( 
+logging.basicConfig(
     format = "[%(asctime)s] %(name)s %(levelname)s: %(message)s",
-    force = True,
+    force=True,
 )
 logging.getLogger("pyschism").setLevel(logging.DEBUG)
 
+def list_of_strings(arg):
+    return arg.split(',')
+
 if __name__ == "__main__":
 
     '''
@@ -24,17 +29,43 @@
     database = 'tpxo'
         ~/.local/share/tpxo/
     '''
+    parser = argparse.ArgumentParser(description="Create bctides.in for SCHISM with command-line arguments! e.g. python test_bctides.py hgrid.gr3 2014-12-01 397 '[[5,5,4,4],[5,5,4,4],[0,1,1,2]]' major fes2014")
+    
+    #Add arguments
+    parser.add_argument('hgrid', type=str, help='hgrid (lon/lat) file')
+    parser.add_argument('start_date', type=datetime.fromisoformat, help='model startdate')
+    parser.add_argument('rnday', type=float, help='model rnday')
+    parser.add_argument('bctypes', type=str, help="JSON format for Flags for each open boundary, '[5,5,4,4],[5,5,4,4],[0,1,1,2]'")
+    parser.add_argument('constituents', type=list_of_strings, help="Choose tidal constituents to be included, major, minor, or list of constituents ['K1', 'O1', 'M2']")
+    parser.add_argument('database', type=str, help='Tidal datbase: tpxo or fes2014')
+
+    #Parse the command-line arguments
+    args = parser.parse_args()
+    hgrid_filename = args.hgrid
+    start_date = args.start_date
+    rnday = args.rnday
+    bctypes = args.bctypes
+    constituents = args.constituents
+    database = args.database
+
+    # Parse the JSON string into a Python data structure
+    try:
+        flags = json.loads(bctypes)
+        print("Parsed bctype list:", flags)
+    except json.JSONDecodeError:
+        raise TypeError("Invalid JSON format for bctype list.")
 
-    start_date = datetime(2014, 12, 1)
-    rnday = 397
+    #start_date = datetime(2014, 12, 1)
+    #rnday = 397
     outdir = './'
-    hgrid = Hgrid.open("./hgrid.gr3", crs="epsg:4326")
+    hgrid = Hgrid.open(hgrid_filename, crs="epsg:4326")
 
     bctides=Bctides(
-        hgrid, 
-        flags = [[5, 5, 4, 4], [5, 5, 4, 4], [0, 1, 2, 2]],
-        constituents = 'major',
-        database = 'fes2014',
+        hgrid = hgrid, 
+        #flags = [[5, 5, 4, 4], [5, 5, 4, 4], [0, 1, 2, 2]],
+        flags = flags,
+        constituents = constituents, 
+        database = database,
         tthconst = 10.0,
         sthconst = 0.0,
     )
diff --git a/pyschism/forcing/bctides/bctides.py b/pyschism/forcing/bctides/bctides.py
index ac6b5a7e..beb7c13c 100644
--- a/pyschism/forcing/bctides/bctides.py
+++ b/pyschism/forcing/bctides/bctides.py
@@ -2,17 +2,19 @@
 from functools import cached_property, lru_cache 
 import pathlib
 import logging
+from typing import Union
 
 from pyschism.forcing.bctides import Tides
 
 logger = logging.getLogger(__name__)
+
 class Bctides:
 
     def __init__(
         self,
         hgrid,
         flags,
-        constituents = 'major',
+        constituents: Union[str, list] = 'major',
         database = 'tpxo',
         cutoff_depth: float = 50.0,
         ethconst = None,
@@ -74,6 +76,8 @@ def __str__(self):
  
         #get amplitude and phase for each open boundary         
         f.append(f"{len(self.gdf)} !nope")
+        if len(self.gdf) != len(self.flags):
+            raise ValueError(f'Number of open boundary {len(self.gdf)} is not consistent with number of given bctypes {len(self.flags)}!')
         for ibnd, (boundary, flag) in enumerate(zip(self.gdf.itertuples(), self.flags)):
             logger.info(f"Processing boundary {ibnd}:")
             #number of nodes and flags
@@ -218,7 +222,7 @@ def write(
         output_directory,
         start_date: datetime = None,
         rnday = None,
-        constituents = ['K1', 'O1'],
+        constituents = 'major',
         overwrite: bool = True, 
     ):
 

From d75f27b1820a6a70101a59cf800e56f2264a188e Mon Sep 17 00:00:00 2001
From: cuill <cui86ll@gmail.com>
Date: Thu, 28 Sep 2023 16:03:38 -0400
Subject: [PATCH 06/15] added prompt in gen_bctides.py

---
 examples/Bctides/gen_bctides.py     | 68 ++++++++++++++++++++-
 pyschism/forcing/bctides/bctides.py | 94 +++++++++++++----------------
 2 files changed, 106 insertions(+), 56 deletions(-)

diff --git a/examples/Bctides/gen_bctides.py b/examples/Bctides/gen_bctides.py
index db8e14ef..ed6f4955 100644
--- a/examples/Bctides/gen_bctides.py
+++ b/examples/Bctides/gen_bctides.py
@@ -3,9 +3,10 @@
 import argparse
 from datetime import datetime, timedelta
 import logging
-import pathlib
 import json
 
+import numpy as np
+
 from pyschism.mesh import Hgrid
 from pyschism.forcing.bctides import Bctides
 
@@ -55,6 +56,62 @@ def list_of_strings(arg):
     except json.JSONDecodeError:
         raise TypeError("Invalid JSON format for bctype list.")
 
+    #earth tidal potential
+    add_earth_tidal_potential = input("Would you like to add earth tidal potential? Y/N: ")
+    if add_earth_tidal_potential == "Y":
+        earth_tidal_potential = True
+    else:
+        earth_tidal_potential = False
+
+    #Check if constant values needed
+    ethconst = []
+    vthconst = []
+    tthconst = []
+    sthconst = []
+    tobc = []
+    sobc = []
+
+    for ibnd, flag in enumerate(flags):
+        iettype, ifltype, itetype, isatype = [i for i in flag]
+        if iettype == 2:
+            val = input(f"Elevation value at boundary {ibnd+1}: ")
+            ethconst.append(float(val))
+        else:
+            ethconst.append(np.nan)
+
+        if ifltype == 2:
+            val = input(f"Discharge value at boundary {ibnd+1}: ")
+            vthconst.append(float(val))
+        else:
+            vthconst.append(np.nan)
+
+        if itetype == 2:
+            val = input(f"Temperature value at boundary {ibnd+1}: ")
+            tthconst.append(float(val))
+            val = input(f"Nuding factor for temperature at boundary {ibnd+1}: ")
+            tobc.append(float(val))
+        elif itetype == 1 or itetype == 3 or itetype == 4:
+            tthconst.append(np.nan)
+            val = input(f"Nuding factor for temperature at boundary {ibnd+1}: ")
+            tobc.append(float(val))
+        else:
+            tthconst.append(np.nan)
+            tobc.append(np.nan)
+
+        if isatype == 2:
+            val = input(f"Salinity value at boundary {ibnd+1}: ")
+            sthconst.append(float(val))
+            val = input(f"Nuding factor for salinity at boundary {ibnd+1}: ")
+            sobc.append(float(val))
+        elif isatype == 1 or isatype == 3 or isatype == 4:
+            sthconst.append(np.nan)
+            val = input(f"Nuding factor for salinity at boundary {ibnd+1}: ")
+            sobc.append(float(val))
+        else:
+            sthconst.append(np.nan)
+            sobc.append(np.nan)
+
+
     #start_date = datetime(2014, 12, 1)
     #rnday = 397
     outdir = './'
@@ -66,8 +123,13 @@ def list_of_strings(arg):
         flags = flags,
         constituents = constituents, 
         database = database,
-        tthconst = 10.0,
-        sthconst = 0.0,
+        add_earth_tidal = earth_tidal_potential,
+        ethconst = ethconst,
+        vthconst = vthconst,
+        tthconst = tthconst,
+        sthconst = sthconst,
+        tobc = tobc,
+        sobc = sobc,
     )
 
     bctides.write(
diff --git a/pyschism/forcing/bctides/bctides.py b/pyschism/forcing/bctides/bctides.py
index beb7c13c..698e3283 100644
--- a/pyschism/forcing/bctides/bctides.py
+++ b/pyschism/forcing/bctides/bctides.py
@@ -16,16 +16,18 @@ def __init__(
         flags,
         constituents: Union[str, list] = 'major',
         database = 'tpxo',
+        add_earth_tidal: bool = True,
         cutoff_depth: float = 50.0,
-        ethconst = None,
-        vthconst = None,
-        tthconst = None,
-        sthconst = None,
-        tobc = None,
-        sobc = None,
+        ethconst: list = None,
+        vthconst: list = None,
+        tthconst: list = None,
+        sthconst: list = None,
+        tobc: list = None,
+        sobc: list = None,
     ):
 
         self.hgrid = hgrid
+        self.add_earth_tidal = add_earth_tidal
         self.cutoff_depth = cutoff_depth
         self.tides = Tides(constituents=constituents, tidal_database=database)
         self.flags = flags
@@ -33,8 +35,8 @@ def __init__(
         self.vthconst = vthconst 
         self.tthconst = tthconst 
         self.sthconst = sthconst 
-        self.tobc = tobc if tobc is not None else 1.0
-        self.sobc = sobc if sobc is not None else 1.0
+        self.tobc = tobc
+        self.sobc = sobc
 
     def __str__(self):
  
@@ -44,8 +46,8 @@ def __str__(self):
         ]
 
         #get earth tidal potential and frequency
-        f.append(f"{self.ntip} {self.cutoff_depth} !number of earth tidal potential, cut-off depth for applying tidal potential")
-        if self.ntip > 0:
+        if self.add_earth_tidal:
+            f.append(f"{self.ntip} {self.cutoff_depth} !number of earth tidal potential, cut-off depth for applying tidal potential")
             for constituent in self.tides.get_active_potential_constituents(): 
                 forcing = self.tides(self.start_date, self.rnday, constituent)
                 f.append(
@@ -79,7 +81,7 @@ def __str__(self):
         if len(self.gdf) != len(self.flags):
             raise ValueError(f'Number of open boundary {len(self.gdf)} is not consistent with number of given bctypes {len(self.flags)}!')
         for ibnd, (boundary, flag) in enumerate(zip(self.gdf.itertuples(), self.flags)):
-            logger.info(f"Processing boundary {ibnd}:")
+            logger.info(f"Processing boundary {ibnd+1}:")
             #number of nodes and flags
             line = [
                 f"{len(boundary.indexes)}",
@@ -99,11 +101,8 @@ def __str__(self):
             if iettype == 1:
                 logger.warning(f'time history of elevation is read in from elev.th (ASCII)!')
             elif iettype == 2:
-                logger.info("This boundary is forced by a constant elevation!")
-                if self.ethconst is None:
-                    raise IOError("You are choosing type 2 for elevation, please specify a value!")
-                else:
-                    f.append("{str(self.ethconst)} !T")
+                logger.info("You are choosing type 2 for elevation, value is {selfethconst[ibnd]} ")
+                f.append(f"{self.ethconst[ibnd]}")
             elif iettype == 4:
                 logger.warning('time history of elevation is read in from elev2D.th.nc (netcdf)')
             elif iettype == 3 or iettype == 5:
@@ -127,11 +126,8 @@ def __str__(self):
             elif ifltype == 1:
                 logger.warning("time history of discharge is read in from flux.th (ASCII)!")
             elif ifltype == 2:
-                logger.info("This boundary is forced by a constant discharge!")
-                if self.vthconst is None:
-                    raise IOError("You are choosing type 2 for velocity, please specify dischage value!")
-                else:
-                    f.append("{str(self.vthconst)}")
+                logger.info("You are choosing type 2 for velocity, value is {self.vthconst[ibnd]} ")
+                f.append(f"{self.vthconst[ibnd]}")
             elif ifltype == 3 or ifltype == 5:
                 if ifltype == 5:
                     logger.warning(f'Combination of 3 and 4, time history of velocity is read in from uv.3D.th.nc!')
@@ -144,10 +140,10 @@ def __str__(self):
             elif ifltype == 4 or -4:
                 logger.warning("time history of velocity (not discharge!) is read in from uv3D.th.nc (netcdf)")
                 if ifltype == -4:
-                    if self.rel1 is None or self.rel2 is None:
-                        raise IOError(f"Please specify relaxation constants for inflow and outflow (between 0 and 1 with 1 being strongest nudging)")
-                    else:
-                        logger.info(f"relax value 1: {self.rel1}, relax value 2: {self.rel2}")
+                    logger.info(f"You are using type 4 , please specify relaxation constants for inflow and outflow (between 0 and 1 with 1 being strongest nudging)")
+                    rel1 = input()
+                    rel2 = input()
+                    f.append(f"{rel1} {rel2}")
             elif ifltype == -1: 
                 raise NotImplementedError(f"Velocity type {ifltype} not implemented yet!")
                 #logger.info(f"Flather type radiation b.c. (iettype must be 0 in this case)!")
@@ -163,25 +159,22 @@ def __str__(self):
                 logger.warning("Temperature is not sepcified, not input needed!")
             elif itetype == 1:
                 logger.warning("time history of temperature will be read in from TEM_1.th!")
-                logger.info(f"Using nudging factor: {self.tobc}!")
-                f.append(f"{self.tobc} !nudging factor for T")
+                logger.info(f"Nudging factor for T at boundary {ibnd+1} is {self.tobc[ibnd]}")
+                f.append(f"{self.tobc[ibnd]} !nudging factor for T")
             elif itetype == 2:
-                logger.info("This boundary is forced by a constant temperature!")
-                if self.tthconst is None:
-                    raise IOError("You are choosing type 2 for temperature, please specify temperature value!")
-                else:
-                    f.append(f"{str(self.tthconst)} !T")
+                logger.info("You are choosing type 2 for temperature, value is {self.tthconst[ibnd]} ")
+                f.append(f"{self.tthconst[ibnd]} !T")
 
-                logger.info(f"Using nudging factor: {self.tobc}!")
-                f.append(f"{self.tobc} !nudging factor for T")
+                logger.info(f"Nudging factor for T at boundary {ibnd+1} is {self.tobc[ibnd]}")
+                f.append(f"{self.tobc[ibnd]} !nudging factor for T")
             elif itetype == 3:
                 logger.info("Using initial temperature profile for inflow")
-                logger.info(f"Using nudging factor: {self.tobc}!")
-                f.append(f"{self.tobc} !nudging factor for T")
+                logger.info(f"Nudging factor for T at boundary {ibnd+1} is{self.tobc[ibnd]}")
+                f.append(f"{self.tobc[ibnd]} !nudging factor for T")
             elif itetype == 4:
                 logger.warning("time history of temperature is read in from TEM_3D.th.nc (netcdf)!")
-                logger.info(f"Using nudging factor: {self.tobc}!")
-                f.append(f"{self.tobc} !nudging factor for T")
+                logger.info(f"Nudging factor for T at boundary {ibnd+1} is{self.tobc[ibnd]}")
+                f.append(f"{self.tobc[ibnd]} !nudging factor for T")
             else:
                 raise IOError(f"Invalid type {itetype} for salinity!")        
 
@@ -191,32 +184,27 @@ def __str__(self):
                 logger.info("Salinity is not sepcified, not input needed!")
             elif isatype == 1:
                 logger.warning("time history of salinity will be read in from SAL_1.th!")
-                logger.info(f"Using nudging factor: {self.sobc}!")
-                f.append(f"{self.sobc} !nudging factor for S")
+                logger.info(f"Nudging factor for salt at boundary {ibnd+1} is {self.sobc[ibnd]}")
+                f.append(f"{self.sobc[ibnd]} !nudging factor for S")
             elif isatype == 2:
-                logger.info("This boundary is forced by a constant salinity!")
-                if self.sthconst is None:
-                    raise IOError("You are choosing type 2 for salinity, please specify salinity value!")
-                else:
-                    f.append(f"{str(self.sthconst)} !S")
+                logger.info("Yor are choosing type 2 for salinity, value is {self.sthconst[ibnd]} ")
+                f.append(f"{self.sthconst[ibnd]} !S")
 
-                logger.info(f"Using nudging factor: {self.sobc}!")
-                f.append(f"{self.sobc} !nudging factor for S")
+                logger.info(f"Nudging factor for salt at boundary {ibnd+1} is {self.sobc[ibnd]}")
+                f.append(f"{self.sobc[ibnd]} !nudging factor for S")
             elif isatype == 3:
                 logger.info("Using initial salinity profile for inflow")
-                logger.info(f"Using nudging factor: {self.sobc}!")
-                f.append(f"{self.sobc} !nudging factor for S")
+                logger.info(f"Nudging factor for salt at boundary {ibnd+1} is {self.sobc[ibnd]}")
+                f.append(f"{self.sobc[ibnd]} !nudging factor for S")
             elif isatype == 4:
                 logger.warning("time history of salinity is read in from SAL_3D.th.nc (netcdf)!")
-                logger.info(f"Using nudging factor: {self.sobc}!")
-                f.append(f"{self.sobc} !nudging factor for S")
+                logger.info(f"Nudging factor for salt at boundary {ibnd+1} is {self.sobc[ibnd]}")
+                f.append(f"{self.sobc[ibnd]} !nudging factor for S")
             else:
                 raise IOError(f"Invalid type {isatype} for salinity!") 
                 
-
         return "\n".join(f)
 
-
     def write(
         self,
         output_directory,

From 46cfa18c034e65137b007e4ada9842aca49da09c Mon Sep 17 00:00:00 2001
From: cuill <cui86ll@gmail.com>
Date: Thu, 28 Sep 2023 17:24:06 -0400
Subject: [PATCH 07/15] more edits on gen_bctides.py

---
 examples/Bctides/gen_bctides.py     | 22 +++++++++++++---------
 pyschism/forcing/bctides/bctides.py | 10 ++++++----
 2 files changed, 19 insertions(+), 13 deletions(-)

diff --git a/examples/Bctides/gen_bctides.py b/examples/Bctides/gen_bctides.py
index ed6f4955..7f45a87b 100644
--- a/examples/Bctides/gen_bctides.py
+++ b/examples/Bctides/gen_bctides.py
@@ -70,6 +70,7 @@ def list_of_strings(arg):
     sthconst = []
     tobc = []
     sobc = []
+    relax = []
 
     for ibnd, flag in enumerate(flags):
         iettype, ifltype, itetype, isatype = [i for i in flag]
@@ -80,19 +81,25 @@ def list_of_strings(arg):
             ethconst.append(np.nan)
 
         if ifltype == 2:
-            val = input(f"Discharge value at boundary {ibnd+1}: ")
+            val = input(f"Discharge value (negative for inflow) at boundary {ibnd+1}: ")
             vthconst.append(float(val))
+        elif ifltype == -4:
+            val = input(f"Relaxation constants (between 0 and 1) for inflow at boundary {ibnd+1}: ")
+            relax.append(float(val))
+            val = input(f"Relaxation constants (between 0 and 1) for outflow at boundary {ibnd+1}: ")
+            relax.append(float(val))
+            
         else:
             vthconst.append(np.nan)
 
         if itetype == 2:
             val = input(f"Temperature value at boundary {ibnd+1}: ")
             tthconst.append(float(val))
-            val = input(f"Nuding factor for temperature at boundary {ibnd+1}: ")
+            val = input(f"Nuding factor (between 0 and 1) for temperature at boundary {ibnd+1}: ")
             tobc.append(float(val))
         elif itetype == 1 or itetype == 3 or itetype == 4:
             tthconst.append(np.nan)
-            val = input(f"Nuding factor for temperature at boundary {ibnd+1}: ")
+            val = input(f"Nuding factor (between 0 and 1) for temperature at boundary {ibnd+1}: ")
             tobc.append(float(val))
         else:
             tthconst.append(np.nan)
@@ -101,25 +108,21 @@ def list_of_strings(arg):
         if isatype == 2:
             val = input(f"Salinity value at boundary {ibnd+1}: ")
             sthconst.append(float(val))
-            val = input(f"Nuding factor for salinity at boundary {ibnd+1}: ")
+            val = input(f"Nuding factor (between 0 and 1) for salinity at boundary {ibnd+1}: ")
             sobc.append(float(val))
         elif isatype == 1 or isatype == 3 or isatype == 4:
             sthconst.append(np.nan)
-            val = input(f"Nuding factor for salinity at boundary {ibnd+1}: ")
+            val = input(f"Nuding factor (between 0 and 1) for salinity at boundary {ibnd+1}: ")
             sobc.append(float(val))
         else:
             sthconst.append(np.nan)
             sobc.append(np.nan)
 
-
-    #start_date = datetime(2014, 12, 1)
-    #rnday = 397
     outdir = './'
     hgrid = Hgrid.open(hgrid_filename, crs="epsg:4326")
 
     bctides=Bctides(
         hgrid = hgrid, 
-        #flags = [[5, 5, 4, 4], [5, 5, 4, 4], [0, 1, 2, 2]],
         flags = flags,
         constituents = constituents, 
         database = database,
@@ -130,6 +133,7 @@ def list_of_strings(arg):
         sthconst = sthconst,
         tobc = tobc,
         sobc = sobc,
+        relax = relax,
     )
 
     bctides.write(
diff --git a/pyschism/forcing/bctides/bctides.py b/pyschism/forcing/bctides/bctides.py
index 698e3283..81bc5b75 100644
--- a/pyschism/forcing/bctides/bctides.py
+++ b/pyschism/forcing/bctides/bctides.py
@@ -24,6 +24,7 @@ def __init__(
         sthconst: list = None,
         tobc: list = None,
         sobc: list = None,
+        relax: list = None,
     ):
 
         self.hgrid = hgrid
@@ -37,6 +38,7 @@ def __init__(
         self.sthconst = sthconst 
         self.tobc = tobc
         self.sobc = sobc
+        self.relax = relax
 
     def __str__(self):
  
@@ -60,6 +62,8 @@ def __str__(self):
                         f"{forcing[4]:G}",
                     ])
                 )
+        else:
+            f.append(f"0 {self.cutoff_depth} !number of earth tidal potential, cut-off depth for applying tidal potential")
 
         #get tidal boundary 
         f.append(f"{self.nbfr:d} !nbfr")
@@ -140,10 +144,8 @@ def __str__(self):
             elif ifltype == 4 or -4:
                 logger.warning("time history of velocity (not discharge!) is read in from uv3D.th.nc (netcdf)")
                 if ifltype == -4:
-                    logger.info(f"You are using type 4 , please specify relaxation constants for inflow and outflow (between 0 and 1 with 1 being strongest nudging)")
-                    rel1 = input()
-                    rel2 = input()
-                    f.append(f"{rel1} {rel2}")
+                    logger.info(f"You are using type -4, relaxation constants for inflow  is {self.relax[0]}, for outflow is {self.relax[1]}")
+                    f.append(f"{self.relax[0]} {self.relax[1]} !relaxation constant")
             elif ifltype == -1: 
                 raise NotImplementedError(f"Velocity type {ifltype} not implemented yet!")
                 #logger.info(f"Flather type radiation b.c. (iettype must be 0 in this case)!")

From f386045011e7f3a78f0b0aa491ec0b4df5c80914 Mon Sep 17 00:00:00 2001
From: cuill <cui86ll@gmail.com>
Date: Fri, 29 Sep 2023 17:29:55 -0400
Subject: [PATCH 08/15] update test scripts

---
 pyschism/forcing/bctides/iettype.py | 192 --------------------
 pyschism/forcing/bctides/ifltype.py | 186 -------------------
 pyschism/forcing/bctides/isatype.py | 106 -----------
 pyschism/forcing/bctides/itetype.py | 122 -------------
 pyschism/forcing/bctides/itrtype.py |  73 --------
 tests/test_barotropic.py            |  19 +-
 tests/test_bctides.py               | 271 ++++++++++++++++++++--------
 tests/test_icogs.py                 |  24 ---
 tests/test_outputs_api.py           |  38 ----
 tests/test_outputs_api_2.py         |  38 ----
 tests/test_param_ncor.py            |   6 +-
 tests/test_paramwind.py             |  84 ++++-----
 12 files changed, 247 insertions(+), 912 deletions(-)
 delete mode 100644 pyschism/forcing/bctides/iettype.py
 delete mode 100644 pyschism/forcing/bctides/ifltype.py
 delete mode 100644 pyschism/forcing/bctides/isatype.py
 delete mode 100644 pyschism/forcing/bctides/itetype.py
 delete mode 100644 pyschism/forcing/bctides/itrtype.py
 delete mode 100644 tests/test_icogs.py
 delete mode 100755 tests/test_outputs_api.py
 delete mode 100755 tests/test_outputs_api_2.py

diff --git a/pyschism/forcing/bctides/iettype.py b/pyschism/forcing/bctides/iettype.py
deleted file mode 100644
index 42598ff9..00000000
--- a/pyschism/forcing/bctides/iettype.py
+++ /dev/null
@@ -1,192 +0,0 @@
-from abc import abstractmethod
-from enum import Enum
-from typing import List
-
-from pyschism.forcing.bctides.bctypes import Bctype
-from pyschism.forcing.bctides.tides import Tides
-from pyschism.forcing import hycom
-
-
-class Iettype(Bctype):
-
-    @property
-    @abstractmethod
-    def iettype(self) -> int:
-        pass
-
-    @property
-    def forcing_digit(self):
-        return self.iettype
-
-
-class UniformTimeHistoryElevation(Iettype):
-    def __init__(self, time_history: List[List[float]], boundaries):
-        self.time_history = time_history
-        self.boundaries = boundaries
-
-    def get_boundary_string(self, hgrid, boundary):
-        return ""
-
-    @property
-    def iettype(self) -> int:
-        return 1
-
-
-class ConstantElevation(Iettype):
-    def __init__(self, values: List[float]):
-        self.values = List[float]
-
-    def get_boundary_string(self, hgrid, boundary):
-        return f"{self.values[boundary.index]:G}"
-
-    @property
-    def iettype(self) -> int:
-        return 2
-
-
-class TidalElevation(Iettype):
-    def __init__(self, constituents='all', database="hamtide"):
-        self.tides = Tides(tidal_database=database, constituents=constituents)
-        
-    def add_constituent(
-            self,
-            name: str,
-            amplitude,
-            angular_frequency,
-            phase=0.,
-            **kwargs,
-    ):
-        self.tides.add_constituent(
-            name=name,
-            angular_frequency=angular_frequency,
-            elevation_amplitude=amplitude,
-            elevation_phase=phase,
-            velocity_amplitude=self.tides._amplitudes['velocity'].get(name, 0.),
-            velocity_phase=self.tides._phases['velocity'].get(name, 0.),
-            **kwargs
-        )
-
-    def get_boundary_string(self, hgrid, boundary, global_constituents: List[str] = None):
-        f = []
-        crs = 'epsg:4326' if hgrid.crs is not None else None
-        if global_constituents is None:
-            for constituent in self.tides.get_active_forcing_constituents():
-                f.append(f"{constituent}")
-                vertices = hgrid.get_xy(crs=crs)[boundary.indexes, :]
-                amp, phase = self.tides.get_elevation(constituent, vertices)
-                for i in range(len(boundary.indexes)):
-                    f.append(f"{amp[i]:.8e} {phase[i]:.8e}")
-        else:
-            required_constituent = self.tides.get_active_forcing_constituents()
-            for constituent in global_constituents:
-                f.append(f"{constituent}")
-                if constituent not in required_constituent:
-                    for i in range(len(boundary.indexes)):
-                        f.append(f"{0:.8e} {0:.8e}")
-                else:
-                    vertices = hgrid.get_xy(crs=crs)[boundary.indexes, :]
-                    amp, phase = self.tides.get_elevation(constituent, vertices)
-                    for i in range(len(boundary.indexes)):
-                        f.append(f"{amp[i]:.8e} {phase[i]:.8e}")
-
-        return "\n".join(f)
-
-    @property
-    def iettype(self):
-        return 3
-
-
-class HarmonicElevation(TidalElevation):
-
-    def __init__(self):
-        self.tides = Tides(constituents=None)
-
-class SpatiallyVaryingTimeHistoryElevation(Iettype):
-    class BaroclinicDatabases(Enum):
-        #RTOFS = hycom.RTOFS
-        GOFS = hycom.GOFS
-
-    def __init__(self, data_source="gofs"):
-
-        if isinstance(data_source, str):
-            data_source = self.BaroclinicDatabases[data_source.upper()].value()
-
-        if not isinstance(data_source, hycom.Hycom):
-            raise TypeError(
-                "Argument data_source must be of type str or type "
-                f"{type(hycom.Hycom)}, not type {type(data_source)}."
-            )
-
-        self.data_component: hycom.HycomComponent = data_source.elevation
-
-    def get_boundary_string(self, hgrid, boundary, global_constituents: List[str] = None):
-        return ""
-
-    def write(self, path, hgrid, start_date, run_days, overwrite: bool = False):
-        self.data_component.write(path, hgrid, start_date, run_days, overwrite)
-
-    @property
-    def iettype(self):
-        return 4
-
-
-class TidalAndSpatiallyVaryingElevationTimeHistory(Iettype):
-    def __init__(
-        self, iettype3: TidalElevation, iettype4: SpatiallyVaryingTimeHistoryElevation
-    ):
-        self.iettype3 = iettype3
-        self.iettype4 = iettype4
-
-    def get_boundary_string(self, hgrid, boundary, global_constituents: List[str] = None):
-        return self.iettype3.get_boundary_string(hgrid, boundary, global_constituents)
-
-    def write(self, path, hgrid, start_date, run_days, overwrite: bool = False):
-        self.data_component.write(path, hgrid, start_date, run_days, overwrite)
-
-    @property
-    def iettype(self):
-        return 5
-
-    @property
-    def tides(self):
-        return self.iettype3.tides
-
-    @property
-    def data_component(self):
-        return self.iettype4.data_component
-
-
-class ZeroElevation(Iettype):
-    def __init__(self):
-        raise NotImplementedError(f"{self.__class__.__name__}")
-
-    def get_boundary_string(self):
-        return ""
-
-    @property
-    def iettype(self):
-        return -1
-
-
-class Iettype1(UniformTimeHistoryElevation):
-    pass
-
-
-class Iettype2(ConstantElevation):
-    pass
-
-
-class Iettype3(TidalElevation):
-    pass
-
-
-class Iettype4(SpatiallyVaryingTimeHistoryElevation):
-    pass
-
-
-class Iettype5(TidalAndSpatiallyVaryingElevationTimeHistory):
-    pass
-
-
-class Iettype_1(ZeroElevation):
-    pass
diff --git a/pyschism/forcing/bctides/ifltype.py b/pyschism/forcing/bctides/ifltype.py
deleted file mode 100644
index cf708039..00000000
--- a/pyschism/forcing/bctides/ifltype.py
+++ /dev/null
@@ -1,186 +0,0 @@
-from abc import abstractmethod
-from enum import Enum
-
-from pyschism.forcing.bctides.bctypes import Bctype
-from pyschism.forcing.bctides.tides import Tides
-from pyschism.forcing import hycom
-
-
-class Ifltype(Bctype):
-
-    @property
-    @abstractmethod
-    def ifltype(self) -> int:
-        """Returns integer representig SCHISM ifltype code for bctides.in"""
-
-    @property
-    def forcing_digit(self):
-        return self.ifltype
-
-
-class UniformTimeHistoryVelocity(Ifltype):
-    def __init__(self, time_history):
-        raise NotImplementedError(f"{self.__class__.__name__}")
-        self.time_history = time_history
-
-    @property
-    def ifltype(self):
-        return 1
-
-
-class ConstantVelocity(Ifltype):
-    def __init__(self, value):
-        self.value = value
-
-    @property
-    def ifltype(self) -> int:
-        return 2
-
-    def get_boundary_string(self, *args, **kwargs):
-        return f'{self.value}'
-
-
-class TidalVelocity(Ifltype):
-    def __init__(self, constituents="all", database="hamtide"):
-        self.tides = Tides(tidal_database=database, constituents=constituents)
-
-    def get_boundary_string(self, hgrid, boundary, global_constituents=None):
-        f = []
-        if global_constituents is None:
-            for constituent in self.tides.get_active_forcing_constituents():
-                f.append(f"{constituent}")
-                vertices = hgrid.get_xy(crs="EPSG:4326")[boundary.indexes, :]
-                uamp, uphase, vamp, vphase = self.tides.get_velocity(constituent, vertices)
-                for i in range(len(vertices)):
-                    f.append(
-                        f"{uamp[i]:.8e} {uphase[i]:.8e} {vamp[i]:.8e} {vphase[i]:.8e}"
-                    )
-        else:
-            required_constituent = self.tides.get_active_forcing_constituents()
-            for constituent in global_constituents:
-                f.append(f"{constituent}")
-                if constituent not in required_constituent:
-                    for i in range(len(boundary.indexes)):
-                        f.append(f"{0:.8e} {0:.8e} {0:.8e} {0:.8e}")
-                else:
-                    vertices = hgrid.get_xy(crs="EPSG:4326")[boundary.indexes, :]
-                    uamp, uphase, vamp, vphase = self.tides.get_velocity(constituent, vertices)
-                    for i in range(len(boundary.indexes)):
-                        f.append(f"{uamp[i]:.8e} {uphase[i]:.8e} {vamp[i]:.8e} {vphase[i]:.8e}")
-
-        return "\n".join(f)
-
-    @property
-    def ifltype(self):
-        return 3
-
-
-class SpatiallyVaryingTimeHistoryVelocity(Ifltype):
-    class BaroclinicDatabases(Enum):
-        #RTOFS = hycom.RTOFS
-        GOFS = hycom.GOFS
-
-    def __init__(self, data_source="gofs"):
-        if isinstance(data_source, str):
-            data_source = self.BaroclinicDatabases[data_source.upper()].value()
-
-        if not isinstance(data_source, hycom.Hycom):
-            raise TypeError(
-                "Argument data_source must be of type str or type "
-                f"{type(hycom.Hycom)}, not type {type(data_source)}."
-            )
-
-        self.data_component: hycom.HycomComponent = data_source.velocity
-
-    def write(self, path, hgrid, vgrid, start_date, run_days, overwrite: bool = False):
-        self.data_component.write(path, hgrid, vgrid, start_date, run_days, overwrite)
-
-    def get_boundary_string(self, *args, **kwargs):
-        return ""
-
-    @property
-    def ifltype(self):
-        return 4
-
-
-class TidalAndSpatiallyVaryingVelocityTimeHistory(Ifltype):
-    def __init__(
-        self,
-        ifltype3: TidalVelocity,
-        ifltype4: SpatiallyVaryingTimeHistoryVelocity,
-    ):
-        self.ifltype3 = ifltype3
-        self.ifltype4 = ifltype4
-
-    def write(self, path, hgrid, vgrid, start_date, run_days, overwrite: bool = False):
-        self.data_component.write(path, hgrid, vgrid, start_date, run_days, overwrite)
-
-    @property
-    def ifltype(self):
-        return 5
-
-    def get_boundary_string(self, *args, **kwargs):
-        return self.ifltype3.get_boundary_string(*args, **kwargs)
-
-    @property
-    def tides(self):
-        return self.ifltype3.tides
-
-    @property
-    def data_component(self):
-        return self.ifltype4.data_component
-
-
-class ZeroVelocity(Ifltype):
-    def __init__(self):
-        raise NotImplementedError(f"{self.__class__.__name__}")
-
-    @property
-    def ifltype(self):
-        return -1
-
-
-class InflowVelocity(Ifltype):
-    def __init__(self):
-        raise NotImplementedError(f"{self.__class__.__name__}")
-
-    @property
-    def ifltype(self):
-        return -4
-
-
-class OutflowVelocity(Ifltype):
-    def __init__(self):
-        raise NotImplementedError(f"{self.__class__.__name__}")
-
-    @property
-    def ifltype(self):
-        return -5
-
-
-class Ifltype1(UniformTimeHistoryVelocity):
-    pass
-
-
-class Ifltype2(ConstantVelocity):
-    pass
-
-
-class Ifltype3(TidalVelocity):
-    pass
-
-
-class Ifltype4(SpatiallyVaryingTimeHistoryVelocity):
-    pass
-
-
-class Ifltype5(TidalAndSpatiallyVaryingVelocityTimeHistory):
-    pass
-
-
-class Ifltype_1(ZeroVelocity):
-    pass
-
-
-class Flather(ZeroVelocity):
-    pass
diff --git a/pyschism/forcing/bctides/isatype.py b/pyschism/forcing/bctides/isatype.py
deleted file mode 100644
index 3ddfcb6d..00000000
--- a/pyschism/forcing/bctides/isatype.py
+++ /dev/null
@@ -1,106 +0,0 @@
-from abc import abstractmethod
-from enum import Enum
-
-from pyschism.forcing.bctides.bctypes import Bctype
-from pyschism.forcing import hycom
-
-
-class Isatype(Bctype):
-
-    @property
-    @abstractmethod
-    def isatype(self) -> int:
-        pass
-
-    @property
-    def forcing_digit(self):
-        return self.isatype
-
-
-class UniformTimeHistorySalinity(Isatype):
-    def __init__(self, time_history):
-        raise NotImplementedError(f"{self.__class__.__name__}")
-        self.time_history = time_history
-
-    @property
-    def isatype(self) -> int:
-        return 1
-
-
-class ConstantSalinity(Isatype):
-    def __init__(self, value: float, nudging_factor: float):
-        self.value = value
-        if not (nudging_factor >= 0) and (nudging_factor <= 1):
-            raise ValueError(f'Argument `nudging_factor` must be get 0 and let 1, but got {nudging_factor}.')
-        self.nudging_factor = nudging_factor
-        super().__init__()
-
-
-    @property
-    def isatype(self) -> int:
-        return 2
-
-    def get_boundary_string(self, *args, **kwargs) -> str:
-        boundary_string = [
-            f'{self.value:0.6f}',
-            f'{self.nudging_factor:0.6f}',
-        ]
-        return '\n'.join(boundary_string)
-
-
-class SalinityInitialConditions(Isatype):
-    def __init__(self):
-        raise NotImplementedError(f"{self.__class__.__name__}")
-
-    @property
-    def isatype(self):
-        return 3
-
-
-class SpatiallyVaryingTimeHistorySalinity(Isatype):
-    class BaroclinicDatabases(Enum):
-        #RTOFS = hycom.RTOFS
-        GOFS = hycom.GOFS
-
-    def __init__(
-        self,
-        data_source="gofs",
-        nudge: bool = True,
-        rlmax=1.5,
-        rnu_day=0.25,
-    ):
-        if isinstance(data_source, str):
-            data_source = self.BaroclinicDatabases[data_source.upper()].value()
-        if not isinstance(data_source, hycom.Hycom):
-            raise TypeError(
-                "Argument data_source must be of type str or type "
-                f"{type(hycom.Hycom)}, not type {type(data_source)}."
-            )
-        self.data_source = data_source
-        self.data_component = data_source.salinity
-        self.nudge = nudge
-        self.rlmax = rlmax
-        self.rnu_day = rnu_day
-
-    def get_boundary_string(self, *args, **kwargs):
-        return "1."
-
-    @property
-    def isatype(self):
-        return 4
-
-
-class Isatype1(UniformTimeHistorySalinity):
-    pass
-
-
-class Isatype2(ConstantSalinity):
-    pass
-
-
-class Isatype3(SalinityInitialConditions):
-    pass
-
-
-class Isatype4(SpatiallyVaryingTimeHistorySalinity):
-    pass
diff --git a/pyschism/forcing/bctides/itetype.py b/pyschism/forcing/bctides/itetype.py
deleted file mode 100644
index c8add4df..00000000
--- a/pyschism/forcing/bctides/itetype.py
+++ /dev/null
@@ -1,122 +0,0 @@
-from abc import abstractmethod
-from enum import Enum
-
-# from typing import Union
-
-from pyschism.forcing.bctides.bctypes import Bctype
-from pyschism.forcing import hycom
-
-# from pyschism.forcing.bctides.nudge import TempNudge
-# from pyschism.mesh.base import Gr3
-
-
-class Itetype(Bctype):
-    @property
-    @abstractmethod
-    def itetype(self) -> int:
-        pass
-
-    @property
-    def forcing_digit(self):
-        return self.itetype
-
-
-class UniformTimeHistoryTemperature(Itetype):
-
-    # def __init__(self, time_history: Dict[str, List[float]], tobc: float = 1.):
-    #     self.time_history = time_history
-    #     super().__init__(tobc)
-
-    @property
-    def itetype(self) -> int:
-        return 1
-
-
-class ConstantTemperature(Itetype):
-
-    def __init__(self, value: float, nudging_factor: float):
-        self.value = value
-        if not (nudging_factor >= 0) and (nudging_factor <= 1):
-            raise ValueError(
-                    'Argument `nudging_factor` must be >= 0 and <= 1,'
-                    f'but got {nudging_factor}.')
-        self.nudging_factor = nudging_factor
-        super().__init__()
-
-    @property
-    def itetype(self) -> int:
-        return 2
-
-    def get_boundary_string(self, *args, **kwargs) -> str:
-        boundary_string = [
-            f'{self.value:0.6f}',
-            f'{self.nudging_factor:0.6f}',
-        ]
-        return '\n'.join(boundary_string)
-
-
-class TemperatureInitialConditions(Itetype):
-
-    # def __init__(self, tobc: float = 1., tth):
-    #     raise NotImplementedError(f'{self.__class__.__name__}')
-
-    @property
-    def itetype(self):
-        return 3
-
-
-class SpatiallyVaryingTimeHistoryTemperature(Itetype):
-    class BaroclinicDatabases(Enum):
-        #RTOFS = hycom.RTOFS
-        GOFS = hycom.GOFS
-
-    def __init__(
-        self,
-        data_source="gofs",
-        nudge: bool = True,
-        rlmax=1.5,
-        rnu_day=0.25,
-    ):
-        if isinstance(data_source, str):
-            data_source = self.BaroclinicDatabases[data_source.upper()].value()
-        if not isinstance(data_source, hycom.Hycom):
-            raise TypeError(
-                "Argument data_source must be of type str or type "
-                f"{type(hycom.Hycom)}, not type {type(data_source)}."
-            )
-        self.data_source = data_source
-        self.data_component = data_source.temperature
-        self.nudge = bool(nudge)
-        self.rlmax = rlmax
-        self.rnu_day = rnu_day
-
-    def get_boundary_string(self, *args, **kwargs):
-        return "1."
-
-    # def write(self, *args, **kwargs):
-    #     if self.nudge is not None:
-    #         self.nudge.write(*args, **kwargs)
-
-    @property
-    def itetype(self):
-        return 4
-
-
-class Itetype0(Itetype):
-    pass
-
-
-class Itetype1(UniformTimeHistoryTemperature):
-    pass
-
-
-class Itetype2(ConstantTemperature):
-    pass
-
-
-class Itetype3(TemperatureInitialConditions):
-    pass
-
-
-class Itetype4(SpatiallyVaryingTimeHistoryTemperature):
-    pass
diff --git a/pyschism/forcing/bctides/itrtype.py b/pyschism/forcing/bctides/itrtype.py
deleted file mode 100644
index f8fbdaf2..00000000
--- a/pyschism/forcing/bctides/itrtype.py
+++ /dev/null
@@ -1,73 +0,0 @@
-from pyschism.forcing.bctides.bctypes import Bctype
-
-
-class Itrtype(Bctype):
-
-    @property
-    def itrtype(self) -> int:
-        '''Returns integer representig SCHISM itrtype code for bctides.in'''
-
-    @property
-    def forcing_digit(self):
-        return self.itrtype
-
-    def get_boundary_string(self, hgrid, boundary):
-        return ''
-
-
-class UniformTimeHistoryTracer(Itrtype):
-
-    def __init__(self, time_history):
-        raise NotImplementedError(f'{self.__class__.__name__}')
-        self.time_history = time_history
-
-    @property
-    def itrtype(self) -> int:
-        return 1
-
-
-class ConstantTracer(Itrtype):
-
-    def __init__(self, value):
-        raise NotImplementedError(f'{self.__class__.__name__}')
-        self.value = value
-
-    @property
-    def itrtype(self) -> int:
-        return 2
-
-
-class TracerInitialConditions(Itrtype):
-
-    def __init__(self):
-        raise NotImplementedError(f'{self.__class__.__name__}')
-
-    @property
-    def itrtype(self):
-        return 3
-
-
-class SpatiallyVaryingTimeHistoryTracer(Itrtype):
-
-    def __init__(self, data_source):
-        self.data_source = data_source
-
-    @property
-    def itrtype(self):
-        return 4
-
-
-class Itrtype1(UniformTimeHistoryTracer):
-    pass
-
-
-class Itrtype2(ConstantTracer):
-    pass
-
-
-class Itrtype3(TracerInitialConditions):
-    pass
-
-
-class Itrtype4(SpatiallyVaryingTimeHistoryTracer):
-    pass
diff --git a/tests/test_barotropic.py b/tests/test_barotropic.py
index 9e6a2f75..a8a72969 100755
--- a/tests/test_barotropic.py
+++ b/tests/test_barotropic.py
@@ -1,12 +1,11 @@
 #! /usr/bin/env python
-from datetime import timedelta
+from datetime import datetime, timedelta
 import logging
 import pathlib
 import shutil
 
 from pyschism import dates
 from pyschism.driver import ModelConfig
-from pyschism.forcing.bctides import iettype, ifltype
 from pyschism.forcing.nws import NWS2, GFS, HRRR
 from pyschism.forcing.source_sink import NWM
 from pyschism.mesh import Hgrid
@@ -25,22 +24,24 @@ def test_barotropic_forecast():
             "https://raw.githubusercontent.com/geomesh/test-data/main/NWM/hgrid.ll",
             crs="epsg:4326",
         ),
-        iettype=iettype.Iettype3(database="tpxo"),
-        ifltype=ifltype.Ifltype3(database="tpxo"),
-        nws=NWS2(GFS(), HRRR()),
+        flags = [[3, 3, 0, 0]],
+        constituents = 'major',
+        database="tpxo",
+        #nws=NWS2(GFS(), HRRR()),
         source_sink=NWM(
             # aggregation_radius=4000.
         ),
     )
 
-    config.forcings.nws.sflux_2.inventory.file_interval = timedelta(hours=6)
+    #config.forcings.nws.sflux_2.inventory.file_interval = timedelta(hours=6)
 
     # create reference dates
-    nearest_cycle = dates.nearest_cycle()
+    #nearest_cycle = dates.nearest_cycle()
+    start_date = datetime(2023, 9, 10)
     spinup_time = timedelta(days=1)
     coldstart = config.coldstart(
-        start_date=nearest_cycle - spinup_time,
-        end_date=nearest_cycle,
+        start_date= start_date,
+        end_date=start_date + timedelta(days=3),
         timestep=300.0,
         dramp=spinup_time,
         dramp_ss=spinup_time,
diff --git a/tests/test_bctides.py b/tests/test_bctides.py
index 50cc09f5..0c09d557 100755
--- a/tests/test_bctides.py
+++ b/tests/test_bctides.py
@@ -13,7 +13,6 @@
 
 from pyschism.cmd.bctides import BctidesCli
 from pyschism.driver import ModelConfig
-from pyschism.forcing.bctides import iettype, ifltype
 from pyschism.mesh import Hgrid
 
 
@@ -37,85 +36,198 @@ def setUp(self):
                 tar.extractall(data_directory / "GulfStreamDevel")
         self.hgrid2D = hgrid2D
 
-    def _test_bctides_cli_2d(self):
-        import tempfile
-        tmpdir = tempfile.TemporaryDirectory()
-        parser = argparse.ArgumentParser()
-        add_bctides(parser.add_subparsers(dest='mode'))
-        args = parser.parse_args(
-            [
-                'bctides',
-                f'{self.hgrid2D}',  # hgrid path
-                '--run-days=5',  # rndays
-                '--hgrid-crs=epsg:4326',
-                f'-o={tmpdir.name}',
-                '--overwrite'
-            ])
-        BctidesCli(args)
-        # with open(tmpdir.name + '/bctides.in') as f:
-        #     for line in f:
-        #         print(line, end='')
-
-    def _test_bctides_cli_2d_w_velo(self):
-        import tempfile
-        tmpdir = tempfile.TemporaryDirectory()
-        parser = argparse.ArgumentParser()
-        add_bctides(parser.add_subparsers(dest='mode'))
-        args = parser.parse_args(
-            [
-                'bctides',
-                f'{self.hgrid2D}',  # hgrid path
-                '--run-days=5',  # rndays,
-                '--include-velocity',
-                '--hgrid-crs=epsg:4326',
-                f'-o={tmpdir.name}',
-                '--overwrite'
-            ])
-        BctidesCli(args)
-        # with open(tmpdir.name + '/bctides.in') as f:
-        #     for line in f:
-        #         print(line, end='')
-
-    def test_bctides_cli_2d_w_elev2d(self):
-        import tempfile
-        tmpdir = tempfile.TemporaryDirectory()
-        parser = argparse.ArgumentParser()
-        add_bctides(parser.add_subparsers(dest='mode'))
-        args = parser.parse_args(
-            [
-                'bctides',
-                f'{self.hgrid2D}',  # hgrid path
-                '--run-days=5',  # rndays,
-                '--iettype-5',
-                # '--ifltype-5',
-                '--hgrid-crs=epsg:4326',
-                f'-o={tmpdir.name}',
-                '--overwrite'
-            ])
-        BctidesCli(args)
-
-    def _test_bctides_cli_3d(self):
-        import tempfile
-        tmpdir = tempfile.TemporaryDirectory()
-        parser = argparse.ArgumentParser()
-        add_bctides(parser.add_subparsers(dest='mode'))
-        args = parser.parse_args(
-            [
-                'bctides',
-                f'{self.hgrid}',  # hgrid path
-                '--run-days=5',  # rndays
-                '--hgrid-crs=epsg:4326',
-                f'-o={tmpdir.name}',
-                '--overwrite'
-            ])
-        BctidesCli(args)
-
-    def test_simple_tidal_setup(self):
+    #def _test_bctides_cli_2d(self):
+    #    import tempfile
+    #    tmpdir = tempfile.TemporaryDirectory()
+    #    parser = argparse.ArgumentParser()
+    #    add_bctides(parser.add_subparsers(dest='mode'))
+    #    args = parser.parse_args(
+    #        [
+    #            'bctides',
+    #            f'{self.hgrid2D}',  # hgrid path
+    #            '--run-days=5',  # rndays
+    #            '--hgrid-crs=epsg:4326',
+    #            f'-o={tmpdir.name}',
+    #            '--overwrite'
+    #        ])
+    #    BctidesCli(args)
+    #    # with open(tmpdir.name + '/bctides.in') as f:
+    #    #     for line in f:
+    #    #         print(line, end='')
+
+    #def _test_bctides_cli_2d_w_velo(self):
+    #    import tempfile
+    #    tmpdir = tempfile.TemporaryDirectory()
+    #    parser = argparse.ArgumentParser()
+    #    add_bctides(parser.add_subparsers(dest='mode'))
+    #    args = parser.parse_args(
+    #        [
+    #            'bctides',
+    #            f'{self.hgrid2D}',  # hgrid path
+    #            '--run-days=5',  # rndays,
+    #            '--include-velocity',
+    #            '--hgrid-crs=epsg:4326',
+    #            f'-o={tmpdir.name}',
+    #            '--overwrite'
+    #        ])
+    #    BctidesCli(args)
+    #    # with open(tmpdir.name + '/bctides.in') as f:
+    #    #     for line in f:
+    #    #         print(line, end='')
+
+    #def test_bctides_cli_2d_w_elev2d(self):
+    #    import tempfile
+    #    tmpdir = tempfile.TemporaryDirectory()
+    #    parser = argparse.ArgumentParser()
+    #    add_bctides(parser.add_subparsers(dest='mode'))
+    #    args = parser.parse_args(
+    #        [
+    #            'bctides',
+    #            f'{self.hgrid2D}',  # hgrid path
+    #            '--run-days=5',  # rndays,
+    #            '--iettype-5',
+    #            # '--ifltype-5',
+    #            '--hgrid-crs=epsg:4326',
+    #            f'-o={tmpdir.name}',
+    #            '--overwrite'
+    #        ])
+    #    BctidesCli(args)
+
+    #def _test_bctides_cli_3d(self):
+    #    import tempfile
+    #    tmpdir = tempfile.TemporaryDirectory()
+    #    parser = argparse.ArgumentParser()
+    #    add_bctides(parser.add_subparsers(dest='mode'))
+    #    args = parser.parse_args(
+    #        [
+    #            'bctides',
+    #            f'{self.hgrid}',  # hgrid path
+    #            '--run-days=5',  # rndays
+    #            '--hgrid-crs=epsg:4326',
+    #            f'-o={tmpdir.name}',
+    #            '--overwrite'
+    #        ])
+    #    BctidesCli(args)
+
+    def test_case1_setup(self):
+        # Inputs: Modify As Needed!
+        start_date = "2018091000" #in YYYYMMDDHH format
+        end_date = "2018091600" #in YYYYMMDDHH format
+
+        dt = timedelta(seconds=150.)
+        nspool = timedelta(minutes=20.)
+
+        start_date = datetime.strptime(start_date, "%Y%m%d%H")
+        end_date = datetime.strptime(end_date, "%Y%m%d%H")
+        rnday = end_date - start_date
+
+        hgrid = Hgrid.open(
+            'https://raw.githubusercontent.com/geomesh/test-data/main/NWM/hgrid.ll',
+            crs=4326
+        )
+
+        config = ModelConfig(
+            hgrid=hgrid,
+            flags=[[3, 3, 0, 0]],
+            constituents = 'major',
+            database = 'tpxo',
+        )
+
+        coldstart = config.coldstart(
+            start_date=start_date,
+            end_date=start_date + rnday,
+            timestep=dt,
+            nspool=timedelta(hours=1),
+            elev=True,
+            dahv=True,
+        )
+
+        with tempfile.TemporaryDirectory() as dn:
+            tmpdir = pathlib.Path(dn)
+            coldstart.write(tmpdir, overwrite=True)
+
+    def test_case2_setup(self):
+        # Inputs: Modify As Needed!
+        start_date = "2018091000" #in YYYYMMDDHH format
+        end_date = "2018091600" #in YYYYMMDDHH format
+
+        dt = timedelta(seconds=150.)
+        nspool = timedelta(minutes=20.)
+
+        start_date = datetime.strptime(start_date, "%Y%m%d%H")
+        end_date = datetime.strptime(end_date, "%Y%m%d%H")
+        rnday = end_date - start_date
+
+        hgrid = Hgrid.open(
+            'https://raw.githubusercontent.com/geomesh/test-data/main/NWM/hgrid.ll',
+            crs=4326
+        )
+
+        config = ModelConfig(
+            hgrid=hgrid,
+            flags=[[5, 5, 0, 0]],
+            constituents = 'major',
+            database = 'tpxo',
+        )
+
+        coldstart = config.coldstart(
+            start_date=start_date,
+            end_date=start_date + rnday,
+            timestep=dt,
+            nspool=timedelta(hours=1),
+            elev=True,
+            dahv=True,
+        )
+
+        with tempfile.TemporaryDirectory() as dn:
+            tmpdir = pathlib.Path(dn)
+            coldstart.write(tmpdir, overwrite=True)
+
+    def test_case3_setup(self):
         # Inputs: Modify As Needed!
         start_date = "2018091000" #in YYYYMMDDHH format
         end_date = "2018091600" #in YYYYMMDDHH format
 
-        tide_source = "tpxo"
+        dt = timedelta(seconds=150.)
+        nspool = timedelta(minutes=20.)
+
+        start_date = datetime.strptime(start_date, "%Y%m%d%H")
+        end_date = datetime.strptime(end_date, "%Y%m%d%H")
+        rnday = end_date - start_date
+
+        hgrid = Hgrid.open(
+            'https://raw.githubusercontent.com/geomesh/test-data/main/NWM/hgrid.ll',
+            crs=4326
+        )
+
+        config = ModelConfig(
+            hgrid=hgrid,
+            flags=[[5, 5, 2, 2]], #constant temperature and salinity
+            constituents = 'major',
+            database = 'tpxo',
+            tthconst = [10.0],   #T value
+            sthconst = [0.0],    #S value
+            tobc = [0.1],   #nuding factor for T
+            sobc = [0.1],   #nuding factor for S
+        )
+
+        coldstart = config.coldstart(
+            start_date=start_date,
+            end_date=start_date + rnday,
+            timestep=dt,
+            nspool=timedelta(hours=1),
+            elev=True,
+            dahv=True,
+        )
+
+        with tempfile.TemporaryDirectory() as dn:
+            tmpdir = pathlib.Path(dn)
+            coldstart.write(tmpdir, overwrite=True)
+
+    def test_case4_setup(self):
+        # Inputs: Modify As Needed!
+        start_date = "2018091000" #in YYYYMMDDHH format
+        end_date = "2018091600" #in YYYYMMDDHH format
 
         dt = timedelta(seconds=150.)
         nspool = timedelta(minutes=20.)
@@ -131,8 +243,11 @@ def test_simple_tidal_setup(self):
 
         config = ModelConfig(
             hgrid=hgrid,
-            iettype=iettype.Iettype3(database=tide_source),
-            ifltype=ifltype.Ifltype3(database=tide_source),
+            flags=[[5, 5, 4, 4]], #constant temperature and salinity
+            constituents = 'major',
+            database = 'tpxo',
+            tobc = [0.1],   #nuding factor for T
+            sobc = [0.1],   #nuding factor for S
         )
 
         coldstart = config.coldstart(
diff --git a/tests/test_icogs.py b/tests/test_icogs.py
deleted file mode 100644
index f541f029..00000000
--- a/tests/test_icogs.py
+++ /dev/null
@@ -1,24 +0,0 @@
-#! /usr/bin/env python
-import unittest
-import os
-
-
-class OutputsApiTestCase(unittest.TestCase):
-
-    def test_outputs_api(self):
-        test_data = f'{os.getenv("ICOGS3D_HGRID")}'
-        # test_data = "/sciclone/pscr/lcui01/ICOGS3D_dev/outputs_run3c"
-        outputs = OutputCollection(test_data)
-        # anim = outputs.elev.animation(vmin=-3, vmax=3, fps=6)
-        from time import time
-        start_total = time()
-        print('Begin time aggregation.')
-        for dt in outputs.elev.flattened_timevector:
-            start_local = time()
-            outputs.elev.aggregate(dt)
-            print(f'took {time()-start_local} to aggregate elev for time {dt}.')
-        print(f'Total aggregation time for variable elev took {time()-start_total}')
-
-
-if __name__ == '__main__':
-    unittest.main()
diff --git a/tests/test_outputs_api.py b/tests/test_outputs_api.py
deleted file mode 100755
index 97c00312..00000000
--- a/tests/test_outputs_api.py
+++ /dev/null
@@ -1,38 +0,0 @@
-#! /usr/bin/env python
-import argparse
-# from datetime import datetime
-import logging
-import os
-import pathlib
-# import shutil
-import tarfile
-import tempfile
-import unittest
-import urllib.request
-
-from pyschism.outputs.outputs import OutputCollection
-
-
-# logging.basicConfig(level=logging.INFO, force=True)
-
-
-class OutputsApiTestCase(unittest.TestCase):
-
-    def test_outputs_api(self):
-        test_data = f"{os.getenv('HOME')}/pyschism/forecast/forecasts/2021-04-02T00:00:00+00:00/outputs/"
-        # test_data = "/sciclone/pscr/lcui01/ICOGS3D_dev/outputs_run3c"
-        outputs = OutputCollection(test_data)
-        # anim = outputs.elev.animation(vmin=-3, vmax=3, fps=6)
-        from time import time
-        start_total = time()
-        print('Begin time aggregation.')
-        for dt in outputs.elev.flattened_timevector:
-            start_local = time()
-            outputs.elev.aggregate(dt)
-            print(f'took {time()-start_local} to aggregate elev for time {dt}.')
-        print(f'Total aggregation time for variable elev took {time()-start_total}')
-
-
-
-if __name__ == '__main__':
-    unittest.main()
diff --git a/tests/test_outputs_api_2.py b/tests/test_outputs_api_2.py
deleted file mode 100755
index abe8e253..00000000
--- a/tests/test_outputs_api_2.py
+++ /dev/null
@@ -1,38 +0,0 @@
-#! /usr/bin/env python
-# import argparse
-# from datetime import datetime
-# import logging
-import os
-# import pathlib
-# import shutil
-# import tarfile
-# import tempfile
-import unittest
-# import urllib.request
-
-from pyschism.outputs.outputs import OutputCollection
-
-
-# logging.basicConfig(level=logging.INFO, force=True)
-
-class OutputsApiTestCase(unittest.TestCase):
-
-    def test_outputs_api(self):
-        test_data = f"{os.getenv('HOME')}/pyschism/forecast/forecasts/2021-04-02T00:00:00+00:00/outputs/"
-        # test_data = "/sciclone/pscr/lcui01/ICOGS3D_dev/outputs_run3c"
-        print("init outputs")
-        from time import time
-        start_total = time()
-        outputs = OutputCollection(test_data)
-        print(f'init took {time()-start_total}')
-        outputs.elev.aggregate_parallel(nprocs=16)
-        # anim = outputs.elev.animation(vmin=-3, vmax=3, fps=6)
-        # for dt in outputs.elev.flattened_timevector:
-        #     start_local = time()
-        #     outputs.elev.aggregate(dt)
-        #     print(f'took {time()-start_local} to aggregate elev for time {dt}.')
-        # print(f'Total aggregation time for variable elev took {time()-start_total}')
-
-
-if __name__ == '__main__':
-    unittest.main()
diff --git a/tests/test_param_ncor.py b/tests/test_param_ncor.py
index 198d5643..c2605695 100644
--- a/tests/test_param_ncor.py
+++ b/tests/test_param_ncor.py
@@ -26,6 +26,7 @@ def test_ncor_zero(self):
 
         config = ModelConfig(
             hgrid,
+            flags = [[3, 3, 0, 0]],
         )
         
         # create reference dates
@@ -41,8 +42,6 @@ def test_ncor_zero(self):
             dramp_ss=spinup_time,
             drampwind=spinup_time,
             nspool=timedelta(hours=1),
-            elev=True,
-            dahv=True,
         )
 
         with tempfile.TemporaryDirectory() as tempdir:
@@ -58,6 +57,7 @@ def test_ncor_one(self):
 
         config = ModelConfig(
             hgrid,
+            flags = [[3, 3, 0, 0]],
         )
         
         # create reference dates
@@ -73,8 +73,6 @@ def test_ncor_one(self):
             dramp_ss=spinup_time,
             drampwind=spinup_time,
             nspool=timedelta(hours=1),
-            elev=True,
-            dahv=True,
         )
 
         with tempfile.TemporaryDirectory() as tempdir:
diff --git a/tests/test_paramwind.py b/tests/test_paramwind.py
index c868f0c7..e3d66ee1 100644
--- a/tests/test_paramwind.py
+++ b/tests/test_paramwind.py
@@ -1,6 +1,7 @@
 import tempfile
 from datetime import timedelta, datetime
 from pathlib import Path
+import unittest
 
 from stormevents.nhc import VortexTrack
 
@@ -8,54 +9,20 @@
 from pyschism.driver import ModelConfig
 from pyschism.forcing.nws import BestTrackForcing
 
+class ModelConfigurationTestCase(unittest.TestCase):
+    def test_paramwind_from_stormname(self):
 
-def test_paramwind_from_stormname():
-
-    hgrid = Hgrid.open(
-        'https://raw.githubusercontent.com/geomesh/test-data/main/NWM/hgrid.ll'
-    )
-
-    meteo = BestTrackForcing(storm='Florence2018')
-
-    config = ModelConfig(
-        hgrid=hgrid,
-        vgrid=None,
-        fgrid=None,
-        iettype=None,
-        ifltype=None,
-        nws=meteo,
-    )
-
-    driver = config.coldstart(
-        start_date=datetime(2018, 9, 8),
-        end_date=datetime(2018, 9, 18),
-        timestep=timedelta(seconds=150),
-        nspool=24,
-    )
-
-    with tempfile.TemporaryDirectory() as dn:
-        tmpdir = Path(dn)
-        driver.write(tmpdir / 'paramwind', overwrite=True)
-
-def test_paramwind_from_file():
-    hgrid = Hgrid.open(
-        'https://raw.githubusercontent.com/geomesh/test-data/main/NWM/hgrid.ll'
-    )
-
-    track = VortexTrack.from_storm_name('Florence', 2018)
-
+        hgrid = Hgrid.open(
+            'https://raw.githubusercontent.com/geomesh/test-data/main/NWM/hgrid.ll'
+        )
 
-    with tempfile.TemporaryDirectory() as dn:
-        tmpdir = Path(dn)
-        track.to_file(tmpdir / 'track.dat')
+        meteo = BestTrackForcing(storm='Florence2018')
 
-        meteo = BestTrackForcing.from_nhc_bdeck(nhc_bdeck=tmpdir / 'track.dat')
         config = ModelConfig(
             hgrid=hgrid,
             vgrid=None,
             fgrid=None,
-            iettype=None,
-            ifltype=None,
+            flags=[[3, 3, 0, 0]],
             nws=meteo,
         )
 
@@ -66,6 +33,39 @@ def test_paramwind_from_file():
             nspool=24,
         )
 
-        driver.write(tmpdir / 'paramwind', overwrite=True)
+        with tempfile.TemporaryDirectory() as dn:
+            tmpdir = Path(dn)
+            driver.write(tmpdir / 'paramwind', overwrite=True)
+
+    def test_paramwind_from_file(self):
+        hgrid = Hgrid.open(
+            'https://raw.githubusercontent.com/geomesh/test-data/main/NWM/hgrid.ll'
+        )
+
+        track = VortexTrack.from_storm_name('Florence', 2018)
+
+
+        with tempfile.TemporaryDirectory() as dn:
+            tmpdir = Path(dn)
+            track.to_file(tmpdir / 'track.dat')
+
+            meteo = BestTrackForcing.from_nhc_bdeck(nhc_bdeck=tmpdir / 'track.dat')
+            config = ModelConfig(
+                hgrid=hgrid,
+                vgrid=None,
+                fgrid=None,
+                flags=[[3, 3, 0, 0]],
+                nws=meteo,
+            )
+
+            driver = config.coldstart(
+                start_date=datetime(2018, 9, 8),
+                end_date=datetime(2018, 9, 18),
+                timestep=timedelta(seconds=150),
+                nspool=24,
+            )
 
+            driver.write(tmpdir / 'paramwind', overwrite=True)
 
+if __name__ == '__main__':
+    unittest.main()

From 7c40ab49560458f073dbabe034566ab185ccb835 Mon Sep 17 00:00:00 2001
From: cuill <cui86ll@gmail.com>
Date: Fri, 29 Sep 2023 17:31:51 -0400
Subject: [PATCH 09/15] delete unused scripts and update example and test
 scripts

---
 examples/Bctides/gen_bctides.py        | 117 ++++----------------
 examples/Bctides/gen_bctides_prompt.py | 144 +++++++++++++++++++++++++
 pyschism/driver.py                     |  50 ++++++---
 pyschism/forcing/bctides/__init__.py   |   7 +-
 pyschism/forcing/bctides/bctides.py    |   8 +-
 5 files changed, 204 insertions(+), 122 deletions(-)
 create mode 100644 examples/Bctides/gen_bctides_prompt.py

diff --git a/examples/Bctides/gen_bctides.py b/examples/Bctides/gen_bctides.py
index 7f45a87b..123c1499 100644
--- a/examples/Bctides/gen_bctides.py
+++ b/examples/Bctides/gen_bctides.py
@@ -1,9 +1,7 @@
 from time import time
 import os
-import argparse
 from datetime import datetime, timedelta
 import logging
-import json
 
 import numpy as np
 
@@ -16,9 +14,6 @@
 )
 logging.getLogger("pyschism").setLevel(logging.DEBUG)
 
-def list_of_strings(arg):
-    return arg.split(',')
-
 if __name__ == "__main__":
 
     '''
@@ -30,96 +25,24 @@ def list_of_strings(arg):
     database = 'tpxo'
         ~/.local/share/tpxo/
     '''
-    parser = argparse.ArgumentParser(description="Create bctides.in for SCHISM with command-line arguments! e.g. python test_bctides.py hgrid.gr3 2014-12-01 397 '[[5,5,4,4],[5,5,4,4],[0,1,1,2]]' major fes2014")
-    
-    #Add arguments
-    parser.add_argument('hgrid', type=str, help='hgrid (lon/lat) file')
-    parser.add_argument('start_date', type=datetime.fromisoformat, help='model startdate')
-    parser.add_argument('rnday', type=float, help='model rnday')
-    parser.add_argument('bctypes', type=str, help="JSON format for Flags for each open boundary, '[5,5,4,4],[5,5,4,4],[0,1,1,2]'")
-    parser.add_argument('constituents', type=list_of_strings, help="Choose tidal constituents to be included, major, minor, or list of constituents ['K1', 'O1', 'M2']")
-    parser.add_argument('database', type=str, help='Tidal datbase: tpxo or fes2014')
-
-    #Parse the command-line arguments
-    args = parser.parse_args()
-    hgrid_filename = args.hgrid
-    start_date = args.start_date
-    rnday = args.rnday
-    bctypes = args.bctypes
-    constituents = args.constituents
-    database = args.database
-
-    # Parse the JSON string into a Python data structure
-    try:
-        flags = json.loads(bctypes)
-        print("Parsed bctype list:", flags)
-    except json.JSONDecodeError:
-        raise TypeError("Invalid JSON format for bctype list.")
-
-    #earth tidal potential
-    add_earth_tidal_potential = input("Would you like to add earth tidal potential? Y/N: ")
-    if add_earth_tidal_potential == "Y":
-        earth_tidal_potential = True
-    else:
-        earth_tidal_potential = False
-
-    #Check if constant values needed
-    ethconst = []
-    vthconst = []
-    tthconst = []
-    sthconst = []
-    tobc = []
-    sobc = []
-    relax = []
-
-    for ibnd, flag in enumerate(flags):
-        iettype, ifltype, itetype, isatype = [i for i in flag]
-        if iettype == 2:
-            val = input(f"Elevation value at boundary {ibnd+1}: ")
-            ethconst.append(float(val))
-        else:
-            ethconst.append(np.nan)
-
-        if ifltype == 2:
-            val = input(f"Discharge value (negative for inflow) at boundary {ibnd+1}: ")
-            vthconst.append(float(val))
-        elif ifltype == -4:
-            val = input(f"Relaxation constants (between 0 and 1) for inflow at boundary {ibnd+1}: ")
-            relax.append(float(val))
-            val = input(f"Relaxation constants (between 0 and 1) for outflow at boundary {ibnd+1}: ")
-            relax.append(float(val))
-            
-        else:
-            vthconst.append(np.nan)
-
-        if itetype == 2:
-            val = input(f"Temperature value at boundary {ibnd+1}: ")
-            tthconst.append(float(val))
-            val = input(f"Nuding factor (between 0 and 1) for temperature at boundary {ibnd+1}: ")
-            tobc.append(float(val))
-        elif itetype == 1 or itetype == 3 or itetype == 4:
-            tthconst.append(np.nan)
-            val = input(f"Nuding factor (between 0 and 1) for temperature at boundary {ibnd+1}: ")
-            tobc.append(float(val))
-        else:
-            tthconst.append(np.nan)
-            tobc.append(np.nan)
-
-        if isatype == 2:
-            val = input(f"Salinity value at boundary {ibnd+1}: ")
-            sthconst.append(float(val))
-            val = input(f"Nuding factor (between 0 and 1) for salinity at boundary {ibnd+1}: ")
-            sobc.append(float(val))
-        elif isatype == 1 or isatype == 3 or isatype == 4:
-            sthconst.append(np.nan)
-            val = input(f"Nuding factor (between 0 and 1) for salinity at boundary {ibnd+1}: ")
-            sobc.append(float(val))
-        else:
-            sthconst.append(np.nan)
-            sobc.append(np.nan)
 
     outdir = './'
-    hgrid = Hgrid.open(hgrid_filename, crs="epsg:4326")
+    hgrid = Hgrid.open('./hgrid.gr3', crs="epsg:4326")
+    start_date = datetime(2014, 12, 1)
+    rnday = 397
+    
+    #bctypes for each open boundary, here shows [ocean bnd1, ocean bnd2, river bnd]
+    flags = [[5, 5, 4, 4], [5, 5, 4, 4], [0, 1, 1, 2]]
+    constituents = 'major'
+    database = 'fes2014'
+    earth_tidal_potential = True
+    #ethconst = [0.5, 0.5, 0.5]
+    #vthconst = [-100, -100, -100]
+    #tthconst = [10.0, 10.0, 10.0]
+    sthconst = [np.nan, np.nan, 0.0]
+    tobc = [0.5, 0.5, 1]
+    sobc = [0.5, 0.5, 1]
+    #relax = [0.5, 0.5]
 
     bctides=Bctides(
         hgrid = hgrid, 
@@ -127,13 +50,13 @@ def list_of_strings(arg):
         constituents = constituents, 
         database = database,
         add_earth_tidal = earth_tidal_potential,
-        ethconst = ethconst,
-        vthconst = vthconst,
-        tthconst = tthconst,
+        #ethconst = ethconst,
+        #vthconst = vthconst,
+        #tthconst = tthconst,
         sthconst = sthconst,
         tobc = tobc,
         sobc = sobc,
-        relax = relax,
+        #relax = relax,
     )
 
     bctides.write(
diff --git a/examples/Bctides/gen_bctides_prompt.py b/examples/Bctides/gen_bctides_prompt.py
new file mode 100644
index 00000000..66970035
--- /dev/null
+++ b/examples/Bctides/gen_bctides_prompt.py
@@ -0,0 +1,144 @@
+from time import time
+import os
+import argparse
+from datetime import datetime, timedelta
+import logging
+import json
+
+import numpy as np
+
+from pyschism.mesh import Hgrid
+from pyschism.forcing.bctides import Bctides
+
+logging.basicConfig(
+    format = "[%(asctime)s] %(name)s %(levelname)s: %(message)s",
+    force=True,
+)
+logging.getLogger("pyschism").setLevel(logging.DEBUG)
+
+def list_of_strings(arg):
+    return arg.split(',')
+
+if __name__ == "__main__":
+
+    '''
+    Assume files are already located in:
+    database='fes2014'
+        ~/.local/share/fes2014/eastward_velocity/
+        ~/.local/share/fes2014/northward_velocity/  
+        ~/.local/share/fes2014/ocean_tide_extrapolated/
+    database = 'tpxo'
+        ~/.local/share/tpxo/
+    '''
+    parser = argparse.ArgumentParser(description="Create bctides.in for SCHISM with command-line arguments! e.g. python test_bctides.py hgrid.gr3 2014-12-01 397 '[[5,5,4,4],[5,5,4,4],[0,1,1,2]]' major fes2014")
+    
+    #Add arguments
+    parser.add_argument('hgrid', type=str, help='hgrid (lon/lat) file')
+    parser.add_argument('start_date', type=datetime.fromisoformat, help='model startdate')
+    parser.add_argument('rnday', type=float, help='model rnday')
+    parser.add_argument('bctypes', type=str, help="JSON format for Flags for each open boundary, '[[5,5,4,4],[5,5,4,4],[0,1,1,2]]'")
+    parser.add_argument('constituents', type=list_of_strings, help="Choose tidal constituents to be included, major, minor, or list of constituents ['K1', 'O1', 'M2']")
+    parser.add_argument('database', type=str, help='Tidal datbase: tpxo or fes2014')
+
+    #Parse the command-line arguments
+    args = parser.parse_args()
+    hgrid_filename = args.hgrid
+    start_date = args.start_date
+    rnday = args.rnday
+    bctypes = args.bctypes
+    constituents = args.constituents
+    database = args.database
+
+    # Parse the JSON string into a Python data structure
+    try:
+        flags = json.loads(bctypes)
+        print("Parsed bctype list:", flags)
+    except json.JSONDecodeError:
+        raise TypeError("Invalid JSON format for bctype list.")
+
+    #earth tidal potential
+    add_earth_tidal_potential = input("Would you like to add earth tidal potential? Y/N: ")
+    if add_earth_tidal_potential == "Y":
+        earth_tidal_potential = True
+    else:
+        earth_tidal_potential = False
+
+    #Check if constant values needed
+    ethconst = []
+    vthconst = []
+    tthconst = []
+    sthconst = []
+    tobc = []
+    sobc = []
+    relax = []
+
+    for ibnd, flag in enumerate(flags):
+        iettype, ifltype, itetype, isatype = [i for i in flag]
+        if iettype == 2:
+            val = input(f"Elevation value at boundary {ibnd+1}: ")
+            ethconst.append(float(val))
+        else:
+            ethconst.append(np.nan)
+
+        if ifltype == 2:
+            val = input(f"Discharge value (negative for inflow) at boundary {ibnd+1}: ")
+            vthconst.append(float(val))
+        elif ifltype == -4:
+            val = input(f"Relaxation constants (between 0 and 1) for inflow at boundary {ibnd+1}: ")
+            relax.append(float(val))
+            val = input(f"Relaxation constants (between 0 and 1) for outflow at boundary {ibnd+1}: ")
+            relax.append(float(val))
+            
+        else:
+            vthconst.append(np.nan)
+
+        if itetype == 2:
+            val = input(f"Temperature value at boundary {ibnd+1}: ")
+            tthconst.append(float(val))
+            val = input(f"Nuding factor (between 0 and 1) for temperature at boundary {ibnd+1}: ")
+            tobc.append(float(val))
+        elif itetype == 1 or itetype == 3 or itetype == 4:
+            tthconst.append(np.nan)
+            val = input(f"Nuding factor (between 0 and 1) for temperature at boundary {ibnd+1}: ")
+            tobc.append(float(val))
+        else:
+            tthconst.append(np.nan)
+            tobc.append(np.nan)
+
+        if isatype == 2:
+            val = input(f"Salinity value at boundary {ibnd+1}: ")
+            sthconst.append(float(val))
+            val = input(f"Nuding factor (between 0 and 1) for salinity at boundary {ibnd+1}: ")
+            sobc.append(float(val))
+        elif isatype == 1 or isatype == 3 or isatype == 4:
+            sthconst.append(np.nan)
+            val = input(f"Nuding factor (between 0 and 1) for salinity at boundary {ibnd+1}: ")
+            sobc.append(float(val))
+        else:
+            sthconst.append(np.nan)
+            sobc.append(np.nan)
+
+    outdir = './'
+    hgrid = Hgrid.open(hgrid_filename, crs="epsg:4326")
+
+    bctides=Bctides(
+        hgrid = hgrid, 
+        flags = flags,
+        constituents = constituents, 
+        database = database,
+        add_earth_tidal = earth_tidal_potential,
+        ethconst = ethconst,
+        vthconst = vthconst,
+        tthconst = tthconst,
+        sthconst = sthconst,
+        tobc = tobc,
+        sobc = sobc,
+        relax = relax,
+    )
+
+    bctides.write(
+        outdir, 
+        start_date=start_date, 
+        rnday=rnday, 
+        overwrite=True,
+    )
diff --git a/pyschism/driver.py b/pyschism/driver.py
index 5a4ed469..fa16612a 100644
--- a/pyschism/driver.py
+++ b/pyschism/driver.py
@@ -19,7 +19,7 @@
 from pyschism.forcing.nws.best_track import BestTrackForcing
 
 # from pyschism.forcing.baroclinic import BaroclinicForcing
-from pyschism.forcing.bctides import Bctides, iettype, ifltype, isatype, itetype
+from pyschism.forcing.bctides import Bctides
 from pyschism.makefile import MakefileDriver
 from pyschism.mesh import Hgrid, Vgrid, Fgrid, gridgr3, prop
 from pyschism.mesh.fgrid import ManningsN, DragCoefficient
@@ -52,7 +52,7 @@ def write(
             self.bctides.write(
                 output_directory,
                 start_date=driver.param.opt.start_date,
-                end_date=driver.param.core.rnday,
+                rnday=driver.param.core.rnday,
                 overwrite=True,
             )
 
@@ -687,10 +687,17 @@ def __init__(
         hgrid: Hgrid,
         vgrid: Vgrid = None,
         fgrid: Fgrid = None,
-        iettype: iettype.Iettype = None,
-        ifltype: ifltype.Ifltype = None,
-        itetype: itetype.Itetype = None,
-        isatype: isatype.Isatype = None,
+        flags = None,
+        constituents = 'major',
+        database = 'tpxo',
+        add_earth_tidal = True,
+        ethconst = None,
+        vthconst = None,
+        tthconst = None,
+        sthconst = None,
+        tobc = None,
+        sobc = None,
+        relax = None, 
         # itrtype: itrtype.Itrtype = None,
         nws: NWS = None,
         source_sink: Union[List[SourceSink], SourceSink] = None,
@@ -708,10 +715,16 @@ def __init__(
         self.hgrid = hgrid
         self.vgrid = vgrid
         self.fgrid = fgrid
-        self.iettype = iettype
-        self.ifltype = ifltype
-        self.itetype = itetype
-        self.isatype = isatype
+        self.flags = flags
+        self.constituents = constituents
+        self.add_earth_tidal = add_earth_tidal
+        self.ethconst = ethconst
+        self.vthconst = vthconst
+        self.tthconst = tthconst
+        self.sthconst = sthconst
+        self.tobc = tobc
+        self.sobc = sobc
+        self.relax = relax
         # self.itrtype = itrtype
         self.nws = nws
         self.source_sink = source_sink
@@ -941,12 +954,17 @@ def forcings(self):
     def bctides(self):
         if not hasattr(self, "_bctides"):
             self._bctides = Bctides(
-                self.hgrid,
-                vgrid=self.vgrid,
-                iettype=self.iettype,
-                ifltype=self.ifltype,
-                isatype=self.isatype,
-                itetype=self.itetype,
+                hgrid = self.hgrid,
+                flags = self.flags,
+                constituents = self.constituents,
+                add_earth_tidal = self.add_earth_tidal,
+                ethconst = self.ethconst,
+                vthconst = self.vthconst,
+                tthconst = self.tthconst,
+                sthconst = self.sthconst,
+                tobc = self.tobc,
+                sobc = self.sobc,
+                relax = self.relax,
                 # cutoff_depth=self.cutoff_depth,
             )
         return self._bctides
diff --git a/pyschism/forcing/bctides/__init__.py b/pyschism/forcing/bctides/__init__.py
index af5abc81..ff13d70b 100644
--- a/pyschism/forcing/bctides/__init__.py
+++ b/pyschism/forcing/bctides/__init__.py
@@ -1,9 +1,4 @@
 from pyschism.forcing.bctides.tides import Tides
-from pyschism.forcing.bctides.iettype import Iettype
-from pyschism.forcing.bctides.ifltype import Ifltype
-from pyschism.forcing.bctides.isatype import Isatype
-from pyschism.forcing.bctides.itetype import Itetype
-from pyschism.forcing.bctides.itrtype import Itrtype
 from pyschism.forcing.bctides.bctides import Bctides
 
-__all__ = ["Bctides", "Tides", "Iettype", "Ifltype", "Isatype", "Itetype", "Itrtype"]
+__all__ = ["Bctides", "Tides"]
diff --git a/pyschism/forcing/bctides/bctides.py b/pyschism/forcing/bctides/bctides.py
index 81bc5b75..9a6e16b9 100644
--- a/pyschism/forcing/bctides/bctides.py
+++ b/pyschism/forcing/bctides/bctides.py
@@ -1,5 +1,5 @@
 from datetime import datetime, timedelta
-from functools import cached_property, lru_cache 
+from functools import cached_property
 import pathlib
 import logging
 from typing import Union
@@ -9,13 +9,15 @@
 logger = logging.getLogger(__name__)
 
 class Bctides:
+    """
+    """
 
     def __init__(
         self,
         hgrid,
-        flags,
+        flags: list = None,
         constituents: Union[str, list] = 'major',
-        database = 'tpxo',
+        database: str = 'tpxo',
         add_earth_tidal: bool = True,
         cutoff_depth: float = 50.0,
         ethconst: list = None,

From 7de0c1833e3f8b34767a2cc3257784043e67ac68 Mon Sep 17 00:00:00 2001
From: cuill <cui86ll@gmail.com>
Date: Sun, 1 Oct 2023 22:16:36 -0400
Subject: [PATCH 10/15] changed to use u/v associated lon/lat instead of
 elevation's in tpxo

---
 examples/Bctides/gen_bctides.py        | 117 ++++++++++++++++----
 examples/Bctides/gen_bctides_prompt.py | 144 -------------------------
 pyschism/forcing/bctides/tpxo.py       |  31 +++++-
 3 files changed, 125 insertions(+), 167 deletions(-)
 delete mode 100644 examples/Bctides/gen_bctides_prompt.py

diff --git a/examples/Bctides/gen_bctides.py b/examples/Bctides/gen_bctides.py
index 123c1499..65f2c204 100644
--- a/examples/Bctides/gen_bctides.py
+++ b/examples/Bctides/gen_bctides.py
@@ -1,7 +1,9 @@
 from time import time
 import os
+import argparse
 from datetime import datetime, timedelta
 import logging
+import json
 
 import numpy as np
 
@@ -14,6 +16,9 @@
 )
 logging.getLogger("pyschism").setLevel(logging.DEBUG)
 
+def list_of_strings(arg):
+    return arg.split(',')
+
 if __name__ == "__main__":
 
     '''
@@ -25,24 +30,96 @@
     database = 'tpxo'
         ~/.local/share/tpxo/
     '''
+    parser = argparse.ArgumentParser(description="Create bctides.in for SCHISM with command-line arguments! e.g. python test_bctides.py hgrid.ll 2014-12-01 397 '[[5,5,4,4],[5,5,4,4],[0,1,1,2]]' major fes2014")
+    
+    #Add arguments
+    parser.add_argument('hgrid', type=str, help='hgrid.ll (lon/lat) file')
+    parser.add_argument('start_date', type=datetime.fromisoformat, help='model startdate')
+    parser.add_argument('rnday', type=float, help='model rnday')
+    parser.add_argument('bctypes', type=str, help="JSON format for Flags for each open boundary, '[[5,5,4,4],[5,5,4,4],[0,1,1,2]]'")
+    parser.add_argument('constituents', type=list_of_strings, help="Choose tidal constituents to be included, major, minor, or list of constituents ['K1', 'O1', 'M2']")
+    parser.add_argument('database', type=str, help='Tidal datbase: tpxo or fes2014')
+
+    #Parse the command-line arguments
+    args = parser.parse_args()
+    hgrid_filename = args.hgrid
+    start_date = args.start_date
+    rnday = args.rnday
+    bctypes = args.bctypes
+    constituents = args.constituents
+    database = args.database
+
+    # Parse the JSON string into a Python data structure
+    try:
+        flags = json.loads(bctypes)
+        print("Parsed bctype list:", flags)
+    except json.JSONDecodeError:
+        raise TypeError("Invalid JSON format for bctype list.")
+
+    #earth tidal potential
+    add_earth_tidal_potential = input("Would you like to add earth tidal potential? Y/N: ")
+    if add_earth_tidal_potential.lower() == "y":
+        earth_tidal_potential = True
+    else:
+        earth_tidal_potential = False
+
+    #Check if constant values needed
+    ethconst = []
+    vthconst = []
+    tthconst = []
+    sthconst = []
+    tobc = []
+    sobc = []
+    relax = []
+
+    for ibnd, flag in enumerate(flags):
+        iettype, ifltype, itetype, isatype = [i for i in flag]
+        if iettype == 2:
+            val = input(f"Elevation value at boundary {ibnd+1}: ")
+            ethconst.append(float(val))
+        else:
+            ethconst.append(np.nan)
+
+        if ifltype == 2:
+            val = input(f"Discharge value (negative for inflow) at boundary {ibnd+1}: ")
+            vthconst.append(float(val))
+        elif ifltype == -4:
+            val = input(f"Relaxation constants (between 0 and 1) for inflow at boundary {ibnd+1}: ")
+            relax.append(float(val))
+            val = input(f"Relaxation constants (between 0 and 1) for outflow at boundary {ibnd+1}: ")
+            relax.append(float(val))
+            
+        else:
+            vthconst.append(np.nan)
+
+        if itetype == 2:
+            val = input(f"Temperature value at boundary {ibnd+1}: ")
+            tthconst.append(float(val))
+            val = input(f"Nuding factor (between 0 and 1) for temperature at boundary {ibnd+1}: ")
+            tobc.append(float(val))
+        elif itetype == 1 or itetype == 3 or itetype == 4:
+            tthconst.append(np.nan)
+            val = input(f"Nuding factor (between 0 and 1) for temperature at boundary {ibnd+1}: ")
+            tobc.append(float(val))
+        else:
+            tthconst.append(np.nan)
+            tobc.append(np.nan)
+
+        if isatype == 2:
+            val = input(f"Salinity value at boundary {ibnd+1}: ")
+            sthconst.append(float(val))
+            val = input(f"Nuding factor (between 0 and 1) for salinity at boundary {ibnd+1}: ")
+            sobc.append(float(val))
+        elif isatype == 1 or isatype == 3 or isatype == 4:
+            sthconst.append(np.nan)
+            val = input(f"Nuding factor (between 0 and 1) for salinity at boundary {ibnd+1}: ")
+            sobc.append(float(val))
+        else:
+            sthconst.append(np.nan)
+            sobc.append(np.nan)
 
     outdir = './'
-    hgrid = Hgrid.open('./hgrid.gr3', crs="epsg:4326")
-    start_date = datetime(2014, 12, 1)
-    rnday = 397
-    
-    #bctypes for each open boundary, here shows [ocean bnd1, ocean bnd2, river bnd]
-    flags = [[5, 5, 4, 4], [5, 5, 4, 4], [0, 1, 1, 2]]
-    constituents = 'major'
-    database = 'fes2014'
-    earth_tidal_potential = True
-    #ethconst = [0.5, 0.5, 0.5]
-    #vthconst = [-100, -100, -100]
-    #tthconst = [10.0, 10.0, 10.0]
-    sthconst = [np.nan, np.nan, 0.0]
-    tobc = [0.5, 0.5, 1]
-    sobc = [0.5, 0.5, 1]
-    #relax = [0.5, 0.5]
+    hgrid = Hgrid.open(hgrid_filename, crs="epsg:4326")
 
     bctides=Bctides(
         hgrid = hgrid, 
@@ -50,13 +127,13 @@
         constituents = constituents, 
         database = database,
         add_earth_tidal = earth_tidal_potential,
-        #ethconst = ethconst,
-        #vthconst = vthconst,
-        #tthconst = tthconst,
+        ethconst = ethconst,
+        vthconst = vthconst,
+        tthconst = tthconst,
         sthconst = sthconst,
         tobc = tobc,
         sobc = sobc,
-        #relax = relax,
+        relax = relax,
     )
 
     bctides.write(
diff --git a/examples/Bctides/gen_bctides_prompt.py b/examples/Bctides/gen_bctides_prompt.py
deleted file mode 100644
index 66970035..00000000
--- a/examples/Bctides/gen_bctides_prompt.py
+++ /dev/null
@@ -1,144 +0,0 @@
-from time import time
-import os
-import argparse
-from datetime import datetime, timedelta
-import logging
-import json
-
-import numpy as np
-
-from pyschism.mesh import Hgrid
-from pyschism.forcing.bctides import Bctides
-
-logging.basicConfig(
-    format = "[%(asctime)s] %(name)s %(levelname)s: %(message)s",
-    force=True,
-)
-logging.getLogger("pyschism").setLevel(logging.DEBUG)
-
-def list_of_strings(arg):
-    return arg.split(',')
-
-if __name__ == "__main__":
-
-    '''
-    Assume files are already located in:
-    database='fes2014'
-        ~/.local/share/fes2014/eastward_velocity/
-        ~/.local/share/fes2014/northward_velocity/  
-        ~/.local/share/fes2014/ocean_tide_extrapolated/
-    database = 'tpxo'
-        ~/.local/share/tpxo/
-    '''
-    parser = argparse.ArgumentParser(description="Create bctides.in for SCHISM with command-line arguments! e.g. python test_bctides.py hgrid.gr3 2014-12-01 397 '[[5,5,4,4],[5,5,4,4],[0,1,1,2]]' major fes2014")
-    
-    #Add arguments
-    parser.add_argument('hgrid', type=str, help='hgrid (lon/lat) file')
-    parser.add_argument('start_date', type=datetime.fromisoformat, help='model startdate')
-    parser.add_argument('rnday', type=float, help='model rnday')
-    parser.add_argument('bctypes', type=str, help="JSON format for Flags for each open boundary, '[[5,5,4,4],[5,5,4,4],[0,1,1,2]]'")
-    parser.add_argument('constituents', type=list_of_strings, help="Choose tidal constituents to be included, major, minor, or list of constituents ['K1', 'O1', 'M2']")
-    parser.add_argument('database', type=str, help='Tidal datbase: tpxo or fes2014')
-
-    #Parse the command-line arguments
-    args = parser.parse_args()
-    hgrid_filename = args.hgrid
-    start_date = args.start_date
-    rnday = args.rnday
-    bctypes = args.bctypes
-    constituents = args.constituents
-    database = args.database
-
-    # Parse the JSON string into a Python data structure
-    try:
-        flags = json.loads(bctypes)
-        print("Parsed bctype list:", flags)
-    except json.JSONDecodeError:
-        raise TypeError("Invalid JSON format for bctype list.")
-
-    #earth tidal potential
-    add_earth_tidal_potential = input("Would you like to add earth tidal potential? Y/N: ")
-    if add_earth_tidal_potential == "Y":
-        earth_tidal_potential = True
-    else:
-        earth_tidal_potential = False
-
-    #Check if constant values needed
-    ethconst = []
-    vthconst = []
-    tthconst = []
-    sthconst = []
-    tobc = []
-    sobc = []
-    relax = []
-
-    for ibnd, flag in enumerate(flags):
-        iettype, ifltype, itetype, isatype = [i for i in flag]
-        if iettype == 2:
-            val = input(f"Elevation value at boundary {ibnd+1}: ")
-            ethconst.append(float(val))
-        else:
-            ethconst.append(np.nan)
-
-        if ifltype == 2:
-            val = input(f"Discharge value (negative for inflow) at boundary {ibnd+1}: ")
-            vthconst.append(float(val))
-        elif ifltype == -4:
-            val = input(f"Relaxation constants (between 0 and 1) for inflow at boundary {ibnd+1}: ")
-            relax.append(float(val))
-            val = input(f"Relaxation constants (between 0 and 1) for outflow at boundary {ibnd+1}: ")
-            relax.append(float(val))
-            
-        else:
-            vthconst.append(np.nan)
-
-        if itetype == 2:
-            val = input(f"Temperature value at boundary {ibnd+1}: ")
-            tthconst.append(float(val))
-            val = input(f"Nuding factor (between 0 and 1) for temperature at boundary {ibnd+1}: ")
-            tobc.append(float(val))
-        elif itetype == 1 or itetype == 3 or itetype == 4:
-            tthconst.append(np.nan)
-            val = input(f"Nuding factor (between 0 and 1) for temperature at boundary {ibnd+1}: ")
-            tobc.append(float(val))
-        else:
-            tthconst.append(np.nan)
-            tobc.append(np.nan)
-
-        if isatype == 2:
-            val = input(f"Salinity value at boundary {ibnd+1}: ")
-            sthconst.append(float(val))
-            val = input(f"Nuding factor (between 0 and 1) for salinity at boundary {ibnd+1}: ")
-            sobc.append(float(val))
-        elif isatype == 1 or isatype == 3 or isatype == 4:
-            sthconst.append(np.nan)
-            val = input(f"Nuding factor (between 0 and 1) for salinity at boundary {ibnd+1}: ")
-            sobc.append(float(val))
-        else:
-            sthconst.append(np.nan)
-            sobc.append(np.nan)
-
-    outdir = './'
-    hgrid = Hgrid.open(hgrid_filename, crs="epsg:4326")
-
-    bctides=Bctides(
-        hgrid = hgrid, 
-        flags = flags,
-        constituents = constituents, 
-        database = database,
-        add_earth_tidal = earth_tidal_potential,
-        ethconst = ethconst,
-        vthconst = vthconst,
-        tthconst = tthconst,
-        sthconst = sthconst,
-        tobc = tobc,
-        sobc = sobc,
-        relax = relax,
-    )
-
-    bctides.write(
-        outdir, 
-        start_date=start_date, 
-        rnday=rnday, 
-        overwrite=True,
-    )
diff --git a/pyschism/forcing/bctides/tpxo.py b/pyschism/forcing/bctides/tpxo.py
index 4834e677..76ad699f 100644
--- a/pyschism/forcing/bctides/tpxo.py
+++ b/pyschism/forcing/bctides/tpxo.py
@@ -70,13 +70,29 @@ def constituents(self):
         return self._constituents
 
     @property
-    def x(self) -> np.ndarray:
+    def lon_z(self) -> np.ndarray:
         return self.h['lon_z'][:, 0].data
 
     @property
-    def y(self) -> np.ndarray:
+    def lat_z(self) -> np.ndarray:
         return self.h['lat_z'][0, :].data
 
+    @property
+    def lon_u(self) -> np.ndarray:
+        return self.uv['lon_u'][:, 0].data
+
+    @property
+    def lat_u(self) -> np.ndarray:
+        return self.uv['lat_u'][0, :].data
+
+    @property
+    def lon_v(self) -> np.ndarray:
+        return self.uv['lon_v'][:, 0].data
+
+    @property
+    def lat_v(self) -> np.ndarray:
+        return self.uv['lat_v'][0, :].data
+
     @property
     def h(self):
         if not hasattr(self, '_h'):
@@ -107,8 +123,17 @@ def _get_interpolation(self, phys_var, ncvar, constituent, vertices):
         lower_c = [c.lower() for c in self.constituents]
         if phys_var == 'elevation':
             ncarray = self.h
-        elif phys_var == 'velocity':
+            self.x = self.lon_z
+            self.y = self.lat_z
+        elif ncvar == 'ua' or ncvar == 'up':
             ncarray = self.uv
+            self.x = self.lon_u
+            self.y = self.lat_u
+        elif ncvar == 'va' or ncvar == 'vp':
+            ncarray = self.uv
+            self.x = self.lon_v
+            self.y = self.lat_v
+
         zi = ncarray[ncvar][
             lower_c.index(constituent.lower()), :, :]
         xo = np.asarray(

From 6d5a33c1c6da4d023f7f79772455470f4e739bed Mon Sep 17 00:00:00 2001
From: cuill <cui86ll@gmail.com>
Date: Tue, 3 Oct 2023 11:29:06 -0400
Subject: [PATCH 11/15] add comments in bctides

---
 pyschism/forcing/bctides/bctides.py | 30 ++++++++++++++++++++++++++++-
 1 file changed, 29 insertions(+), 1 deletion(-)

diff --git a/pyschism/forcing/bctides/bctides.py b/pyschism/forcing/bctides/bctides.py
index 9a6e16b9..3a41bcfa 100644
--- a/pyschism/forcing/bctides/bctides.py
+++ b/pyschism/forcing/bctides/bctides.py
@@ -9,7 +9,9 @@
 logger = logging.getLogger(__name__)
 
 class Bctides:
-    """
+    """Bctides class
+    This is the bctides class to generate bctides.in file.
+    
     """
 
     def __init__(
@@ -28,6 +30,23 @@ def __init__(
         sobc: list = None,
         relax: list = None,
     ):
+    """Initialize Bctides ojbect
+    Parameters
+    --------
+    hgrid: Hgrid object
+    flags: nested list of bctypes
+    constituents: str or list
+    dtabase: str ('tpxo' or 'fes2014')
+    add_earth_tidal: bool
+    cutoff_depth: float
+    ethconst: list (constant elevation value for each open boundary)
+    vthconst: list (constant discharge value for each open boundary)
+    tthconst: list (constant temperature value for each open boundary)
+    sthconst: list (constant salinity value for each open boundary)
+    tobc: list (nuding factor of temperature for each open boundary)
+    sobc: list (nuding factor of salinity for each open boundary)
+    realx: list (relaxation constants for inflow and outflow)
+    """
 
         self.hgrid = hgrid
         self.add_earth_tidal = add_earth_tidal
@@ -217,6 +236,15 @@ def write(
         constituents = 'major',
         overwrite: bool = True, 
     ):
+    """
+    parameters
+    --------
+    output_directory: str
+    start_date: datetime.datetime
+    rnday: int, float or datetime.timedelta
+    constituents: str or list
+    overwrite: bool 
+    """
 
         if start_date is not None:
             self.start_date = start_date

From 25bb6bd3d727e2a86230cbb25e6f6bdc6772c3df Mon Sep 17 00:00:00 2001
From: cuill <cui86ll@gmail.com>
Date: Tue, 3 Oct 2023 16:45:13 -0400
Subject: [PATCH 12/15] fixed IndentationError

---
 pyschism/forcing/bctides/bctides.py | 52 ++++++++++++++---------------
 1 file changed, 26 insertions(+), 26 deletions(-)

diff --git a/pyschism/forcing/bctides/bctides.py b/pyschism/forcing/bctides/bctides.py
index 3a41bcfa..9d8e6075 100644
--- a/pyschism/forcing/bctides/bctides.py
+++ b/pyschism/forcing/bctides/bctides.py
@@ -30,23 +30,23 @@ def __init__(
         sobc: list = None,
         relax: list = None,
     ):
-    """Initialize Bctides ojbect
-    Parameters
-    --------
-    hgrid: Hgrid object
-    flags: nested list of bctypes
-    constituents: str or list
-    dtabase: str ('tpxo' or 'fes2014')
-    add_earth_tidal: bool
-    cutoff_depth: float
-    ethconst: list (constant elevation value for each open boundary)
-    vthconst: list (constant discharge value for each open boundary)
-    tthconst: list (constant temperature value for each open boundary)
-    sthconst: list (constant salinity value for each open boundary)
-    tobc: list (nuding factor of temperature for each open boundary)
-    sobc: list (nuding factor of salinity for each open boundary)
-    realx: list (relaxation constants for inflow and outflow)
-    """
+        """Initialize Bctides ojbect
+        Parameters
+        --------
+        hgrid: Hgrid object
+        flags: nested list of bctypes
+        constituents: str or list
+        database: str ('tpxo' or 'fes2014')
+        add_earth_tidal: bool
+        cutoff_depth: float
+        ethconst: list (constant elevation value for each open boundary)
+        vthconst: list (constant discharge value for each open boundary)
+        tthconst: list (constant temperature value for each open boundary)
+        sthconst: list (constant salinity value for each open boundary)
+        tobc: list (nuding factor of temperature for each open boundary)
+        sobc: list (nuding factor of salinity for each open boundary)
+        realx: list (relaxation constants for inflow and outflow)
+        """
 
         self.hgrid = hgrid
         self.add_earth_tidal = add_earth_tidal
@@ -236,15 +236,15 @@ def write(
         constituents = 'major',
         overwrite: bool = True, 
     ):
-    """
-    parameters
-    --------
-    output_directory: str
-    start_date: datetime.datetime
-    rnday: int, float or datetime.timedelta
-    constituents: str or list
-    overwrite: bool 
-    """
+        """
+        parameters
+        --------
+        output_directory: str
+        start_date: datetime.datetime
+        rnday: int, float or datetime.timedelta
+        constituents: str or list
+        overwrite: bool 
+        """
 
         if start_date is not None:
             self.start_date = start_date

From 5c0b3cf2ff2ec8ceb71fcc6bf28024269d8e8a57 Mon Sep 17 00:00:00 2001
From: cuill <cui86ll@gmail.com>
Date: Thu, 5 Oct 2023 10:39:22 -0400
Subject: [PATCH 13/15] made edits on hrrr2

---
 examples/Sflux/gen_sflux_hrrr2.py  |  31 ++++
 pyschism/forcing/nws/nws2/hrrr2.py | 251 +++++++++++++++++++++--------
 2 files changed, 216 insertions(+), 66 deletions(-)
 create mode 100644 examples/Sflux/gen_sflux_hrrr2.py

diff --git a/examples/Sflux/gen_sflux_hrrr2.py b/examples/Sflux/gen_sflux_hrrr2.py
new file mode 100644
index 00000000..20a62323
--- /dev/null
+++ b/examples/Sflux/gen_sflux_hrrr2.py
@@ -0,0 +1,31 @@
+from time import time
+from datetime import datetime
+import pathlib
+
+from pyschism.mesh.hgrid import Hgrid
+from pyschism.forcing.nws.nws2.hrrr2 import AWSZarrInventory
+
+from pyschism.dates import nearest_cycle
+
+'''
+Please use it with caution, some variables may have missing data, see:
+https://mesowest.utah.edu/html/hrrr/zarr_documentation/html/zarr_variables.html
+'''
+
+if __name__ == "__main__":
+    t0 = time()
+    startdate = datetime(2023, 10, 3)
+
+    #bbox from hgrid
+    #hgrid = Hgrid.open('./hgrid.gr3', crs='epsg:4326')
+    #bbox = hgrid.bbox
+    #or
+    ##bbox from list [xmin, ymin, xmax, ymax]
+    bbox = [-98, 7.8, -52.986, 42.09]
+
+    outdir = './'
+    
+    hrrr = AWSZarrInventory(bbox=bbox)
+    hrrr.write(outdir=outdir, start_date=startdate, rnday=2, air=True, prc=True, rad=False)
+    print(f'It took {(time() - t0)/60} mins')
+
diff --git a/pyschism/forcing/nws/nws2/hrrr2.py b/pyschism/forcing/nws/nws2/hrrr2.py
index d56f8fdd..4a613648 100644
--- a/pyschism/forcing/nws/nws2/hrrr2.py
+++ b/pyschism/forcing/nws/nws2/hrrr2.py
@@ -2,30 +2,54 @@
 import pathlib
 import os
 from typing import Union
+import logging
 
 import numpy as np
 import pandas as pd
 import xarray as xr
+import boto3
+from botocore import UNSIGNED
+from botocore.config import Config
 import fsspec
 import zarr
+from matplotlib.transforms import Bbox
 
 from pyschism.dates import nearest_cycle
 
+logger = logging.getLogger(__name__)
+
 class AWSZarrInventory:
-    '''
-    12/21/2021, L. Cui
-    There is an issue with DSWRF in this dataset. After communicating with 
-    MesoWest group, it is confirmed that they are unable to successfully 
-    generate zarr-formatted output for this vaiable for the zarr-formatted
-    forecast files. At this point, I'll keep this script in the repo. 
+    '''HRRR AWS Zarr bucket
+    Create sflux from HRRr Zarr files.
+
+    Notes
+    -------
+    DSWRF has missing data:
+    https://mesowest.utah.edu/html/hrrr/zarr_documentation/html/zarr_variables.html
+    If only need air/prc, this method is faster than generating sflux from grib2 files (hrrr3).
     '''
 
-    def __init__(self, bbox=None):
-        self.bbox = bbox
+    def __init__(self, bbox: Union[list, Bbox]=None):
+        '''
+        Parameters
+        -------
+        bbox: matplotlib.transforms.Bbox object or a list [xmin, ymin, xmax, ymax]
+        '''
 
-    def gen_sflux(self, outdir: Union[str, os.PathLike], start_date=None, air=True, prc=True, rad=True):
+        if bbox is None:
+            raise ValueError('Bbox is needed!')
 
-        start_date = nearest_cycle() if start_date is None else start_date
+        if not isinstance(bbox, Bbox) and not isinstance(bbox, list):
+            raise TypeError('bbox is a list or matplotlib.transforms.Bbox object!')
+
+        if isinstance(bbox, list):
+            xmin, ymin, xmax, ymax = [v for v in bbox]
+            self.bbox =  Bbox.from_extents(xmin, ymin, xmax, ymax)
+        else:
+            self.bbox = bbox
+
+    def gen_sflux(self, outdir: Union[str, os.PathLike], start_date=None, stack=None, air=True, prc=True, rad=True): 
+        self.start_date = nearest_cycle() if start_date is None else start_date
         airVars = {'stmp': ['TMP', '2m_above_ground'],
                    'spfh': ['SPFH', '2m_above_ground'],
                    'uwind': ['UGRD', '10m_above_ground'],
@@ -35,66 +59,73 @@ def gen_sflux(self, outdir: Union[str, os.PathLike], start_date=None, air=True,
         radVars = {'dlwrf': ['DLWRF', 'surface'],
                    'dswrf': ['DSWRF', 'surface']}
 
-        lon, lat, idx_ymin, idx_ymax, idx_xmin, idx_xmax = self.modified_latlon()
-
+        lon, lat, idx_ymin, idx_ymax, idx_xmin, idx_xmax = self.modified_latlon() 
         path = pathlib.Path(outdir)
         path.mkdir(parents=True, exist_ok=True)
 
-        #for cycle in range(0, 24, 6):
-
-        cycle = start_date.hour
-        base_url = f's3://hrrrzarr/sfc/{start_date.strftime("%Y%m%d")}/' \
-            + f'{start_date.strftime("%Y%m%d")}_{cycle:02d}z_fcst.zarr'
-
-        xdat = []
-        lon2 = xr.DataArray(
-            data=lon, 
-            dims=("ny_grid", "nx_grid",), 
-            name="lon",
-            attrs=dict(
-                long_name = "Longitude",
-                standard_name = "longitude",
-                units = "degrees_east"
-            )
-        )
-        xdat.append(lon2)
-        lat2 = xr.DataArray(
-            data=lat, 
-            dims=("ny_grid", "nx_grid",), 
-            name="lat",
-            attrs=dict(
-                long_name = "Latitude",
-                standard_name = "latitude",
-                units = "degrees_north"
-            )
-        )
-        xdat.append(lat2)
+        self.cycle = start_date.hour
+        base_url = f's3://hrrrzarr/sfc/{self.start_date.strftime("%Y%m%d")}/' \
+            + f'{self.start_date.strftime("%Y%m%d")}_{self.cycle:02d}z_fcst.zarr'
 
         if air:
+            xdat = []
+            xdat.append(self.lon_DataArray(lon1D=lon))
+            xdat.append(self.lat_DataArray(lat1D=lat))
             for key, value in airVars.items():
                 url = f'{base_url}/{value[1]}/{value[0]}/{value[1]}/'
-                print(url)
-                ds = xr.open_zarr(fsspec.get_mapper(url, anon=True))
+                logger.info(url)
+                ds = xr.open_zarr(fsspec.get_mapper(url, anon=True), consolidated=False)
                 ds = ds.rename({'projection_y_coordinate': 'ny_grid'})
                 ds = ds.rename({'projection_x_coordinate': 'nx_grid'})
                 val = ds[value[0]][:, idx_ymin:idx_ymax+1, idx_xmin:idx_xmax+1].astype('float32')
                 xdat.append(xr.DataArray(val, dims=("time", "ny_grid", "nx_grid",), name=key))
 
+                #add variable time
+                times = ds['time'].values.astype('float32')
+               
+                #zarr data already discard t00z, which has all junk values
+                times = (times+1)/24
+                xdat.append(self.time_DataArray(time1D=times))
+
+            fout = xr.merge(xdat)
+            fout = fout.rename_vars({'time2': 'time'})
+            #fout.to_netcdf(path / f'hrrr_{start_date.strftime("%Y%m%d")}{cycle:02d}.nc', 'w', unlimited_dims='time')
+            fout.to_netcdf(path / f'sflux_air_2.{stack:04d}.nc', 'w', unlimited_dims='time')
+            fout.close()
+
         if prc:
+            xdat = []
+            xdat.append(self.lon_DataArray(lon1D=lon))
+            xdat.append(self.lat_DataArray(lat1D=lat))
             for key, value in prcVars.items():
                 url = f'{base_url}/{value[1]}/{value[0]}/{value[1]}/'
-                print(url)
-                ds = xr.open_zarr(fsspec.get_mapper(url, anon=True))
+                logger.info(url)
+                ds = xr.open_zarr(fsspec.get_mapper(url, anon=True), consolidated=False)
                 ds = ds.rename({'projection_y_coordinate': 'ny_grid'})
                 ds = ds.rename({'projection_x_coordinate': 'nx_grid'})
                 val = ds[value[0]][:, idx_ymin:idx_ymax+1, idx_xmin:idx_xmax+1].astype('float32')
                 xdat.append(xr.DataArray(val, dims=("time", "ny_grid", "nx_grid",), name=key))
 
+                #add variable time
+                times = ds['time'].values.astype('float32')
+                times = (times+1)/24
+                #bdate = pd.to_datetime(start_date + timedelta(hours=cycle)).strftime('%Y %m %d %H').split(' ')
+                #bdate = [int(q) for q in bdate[:4]] + [0]
+                xdat.append(self.time_DataArray(time1D=times))
+            fout = xr.merge(xdat)
+            fout = fout.rename_vars({'time2': 'time'})
+            #fout.to_netcdf(path / f'hrrr_{start_date.strftime("%Y%m%d")}{cycle:02d}.nc', 'w', unlimited_dims='time')
+            fout.to_netcdf(path / f'sflux_prc_2.{stack:04d}.nc', 'w', unlimited_dims='time')
+            fout.close()
+
         if rad:
+            xdat = []
+            xdat.append(self.lon_DataArray(lon1D=lon))
+            xdat.append(self.lat_DataArray(lat1D=lat))
             for key, value in radVars.items():
                 url = f'{base_url}/{value[1]}/{value[0]}/{value[1]}/'
-                print(url)
-                ds = xr.open_zarr(fsspec.get_mapper(url, anon=True))
+                logger.info(url)
+                ds = xr.open_zarr(fsspec.get_mapper(url, anon=True), consolidated=False)
                 ds = ds.rename({'projection_y_coordinate': 'ny_grid'})
                 ds = ds.rename({'projection_x_coordinate': 'nx_grid'})
                 val = ds[value[0]][:, idx_ymin:idx_ymax+1, idx_xmin:idx_xmax+1].astype('float32')
@@ -102,28 +133,53 @@ def gen_sflux(self, outdir: Union[str, os.PathLike], start_date=None, air=True,
                 #if np.sum(np.isnan(val)) > 0:
                 #    val[:] = 0.0
                 xdat.append(xr.DataArray(val, dims=("time", "ny_grid", "nx_grid",), name=key))
+
+                #add variable time
+                times = ds['time'].values.astype('float32')
+                times = (times+1)/24
+                #bdate = pd.to_datetime(start_date + timedelta(hours=cycle)).strftime('%Y %m %d %H').split(' ')
+                #bdate = [int(q) for q in bdate[:4]] + [0]
+                xdat.append(self.time_DataArray(time1D=times))
+            fout = xr.merge(xdat)
+            fout = fout.rename_vars({'time2': 'time'})
+            #fout.to_netcdf(path / f'hrrr_{start_date.strftime("%Y%m%d")}{cycle:02d}.nc', 'w', unlimited_dims='time')
+            fout.to_netcdf(path / f'sflux_rad_2.{stack:04d}.nc', 'w', unlimited_dims='time')
+            fout.close()
         
-        time = ds['time'].values.astype('float32')
-        time = (time+1)/24
-        bdate = pd.to_datetime(start_date + timedelta(hours=cycle)).strftime('%Y %m %d %H').split(' ')
-        bdate = [int(q) for q in bdate[:4]] + [0]
-        xdat.append(
-            xr.DataArray(
-                time, dims=('time'), 
-                name='time2',
-                attrs=dict(
-                    long_name="Time",
-                    standard_name = "time",
-                    base_date = bdate,
-                    units = f"days since {start_date.year}-{start_date.month}-{start_date.day} {cycle:02d}:00 UTC"
-                )))
-        fout = xr.merge(xdat)
-        fout=fout.rename_vars({'time2': 'time'})
-        fout.to_netcdf(path / f'hrrr_{start_date.strftime("%Y%m%d")}{cycle:02d}.nc', 'w', unlimited_dims='time')
+        #fout = xr.merge(xdat)
+        #fout=fout.rename_vars({'time2': 'time'})
+        #fout.to_netcdf(path / f'hrrr_{start_date.strftime("%Y%m%d")}{cycle:02d}.nc', 'w', unlimited_dims='time')
+
+    def write(
+        self,
+        outdir,
+        start_date: datetime = None,
+        rnday: Union[float, timedelta] = 2,
+        air: bool = True,
+        prc: bool = False,
+        rad: bool = False,
+    ):
+
+        start_date = nearest_cycle() if start_date is None else start_date
+
+        if not isinstance(rnday, timedelta):
+            try:
+                rnday = timedelta(days=rnday)
+            except:
+                raise TypeError("rnday must be either float or datetime.timedelta!")
+
+        timevector = np.arange(
+            start_date,
+            rnday,
+            np.timedelta64(1, "D"),
+            dtype="datetime64",
+        ).astype(datetime)
+
+        for i, date in enumerate(timevector):
+            self.gen_sflux(outdir, date, i+1, air=air, prc=prc, rad=rad)
 
     def modified_latlon(self):
-        latlon_h5_file = "HRRR_latlon.h5"
-        ds = xr.open_dataset(latlon_h5_file)
+        ds = xr.open_dataset(self.grid_file)
         lon = ds['longitude'].astype("float32")
         lat = ds['latitude'].astype("float32")
         xmin = self.bbox.xmin
@@ -140,3 +196,66 @@ def modified_latlon(self):
         idx_xmax = np.max(idxs[:,1])
         #print(f'idx_ymin is {idx_ymin}, idx_ymax is {idx_ymax}, idx_xmin is {idx_xmin}, idx_xmax is {idx_xmax}')
         return lon[idx_ymin:idx_ymax+1, idx_xmin:idx_xmax+1], lat[idx_ymin:idx_ymax+1, idx_xmin:idx_xmax+1], idx_ymin, idx_ymax, idx_xmin, idx_xmax
+
+
+    def lon_DataArray(self, lon1D):
+        return xr.DataArray(
+            data=lon1D, 
+            dims=("ny_grid", "nx_grid",), 
+            name="lon",
+            attrs=dict(
+                long_name = "Longitude",
+                standard_name = "longitude",
+                units = "degrees_east"
+            )
+        )
+
+    def lat_DataArray(self, lat1D):
+        return xr.DataArray(
+            data=lat1D, 
+            dims=("ny_grid", "nx_grid",), 
+            name="lat",
+            attrs=dict(
+                long_name = "Latitude",
+                standard_name = "latitude",
+                units = "degrees_north"
+            )
+        )
+
+    def time_DataArray(self, time1D):
+        #bdate = pd.to_datetime(self.start_date + timedelta(hours=self.cycle)).strftime('%Y %m %d %H').split(' ')
+        bdate = pd.to_datetime(self.start_date).strftime('%Y %m %d %H').split(' ')
+        bdate = [int(q) for q in bdate[:4]] + [0]
+        
+        return xr.DataArray(
+            time1D, 
+            dims=('time'), 
+            name='time2',
+            attrs=dict(
+                long_name="Time",
+                standard_name = "time",
+                base_date = bdate,
+                units = f"days since {self.start_date.year}-{self.start_date.month}-{self.start_date.day} {self.cycle:02d}:00 UTC"
+            )
+        )
+
+    @property
+    def s3(self):
+        try:
+            return self._s3
+        except AttributeError:
+            self._s3 = boto3.client(
+                "s3", config = Config(signature_version=UNSIGNED))
+            return self._s3
+
+    @property
+    def bucket(self):
+        return "hrrrzarr"
+
+    @property
+    def grid_file(self):
+        #if not hasattr(self, "_grid_file"):
+        self._grid_file = "HRRR_latlon.h5"
+        with open(self._grid_file, "wb") as f:
+            self.s3.download_fileobj(f"{self.bucket}", "grid/HRRR_latlon.h5", f)
+        return self._grid_file

From 9638daa1980bcf485f41409d332214f28a3b2f3d Mon Sep 17 00:00:00 2001
From: cuill <cui86ll@gmail.com>
Date: Wed, 11 Oct 2023 16:24:05 -0400
Subject: [PATCH 14/15] limit CPUs in gfs2 and hrrr3, and comment out
 non-working cli

---
 pyschism/cmd/bctides.py             |    6 +-
 pyschism/cmd/bootstrap/bootstrap.py |    2 +-
 pyschism/cmd/common.py              | 1106 +++++++++++++--------------
 pyschism/cmd/forecast.py            |    4 +-
 pyschism/forcing/nws/nws2/gfs2.py   |   17 +-
 pyschism/forcing/nws/nws2/hrrr3.py  |   18 +-
 6 files changed, 580 insertions(+), 573 deletions(-)

diff --git a/pyschism/cmd/bctides.py b/pyschism/cmd/bctides.py
index 2154bd68..b1379acc 100644
--- a/pyschism/cmd/bctides.py
+++ b/pyschism/cmd/bctides.py
@@ -71,10 +71,10 @@ def add_bctides_options_to_parser(parser):
     common.add_vgrid_to_parser(parser)
     common.add_tidal_constituents_to_parser(parser)
     common.add_tidal_database_to_parser(parser)
-    common.add_baroclinic_database_to_parser(parser)
+    #common.add_baroclinic_database_to_parser(parser)
     common.add_bctides_options_to_parser(parser)
-    common.add_ibctype_to_parser(parser)
-    parser.add_argument("--parallel-download", action='store_true')
+    #common.add_ibctype_to_parser(parser)
+    #parser.add_argument("--parallel-download", action='store_true')
 
 
 #     add_nudge_to_parser(parser)
diff --git a/pyschism/cmd/bootstrap/bootstrap.py b/pyschism/cmd/bootstrap/bootstrap.py
index 3f95d97e..8ab541d2 100644
--- a/pyschism/cmd/bootstrap/bootstrap.py
+++ b/pyschism/cmd/bootstrap/bootstrap.py
@@ -65,4 +65,4 @@ def add_bootstrap_options_to_parser(parser):
     )
     parser.add_argument('--Z0', type=float)
     parser.add_argument("--cutoff-depth", type=float, default=50.0)
-    common.add_ibctype_to_parser(parser)
+    #common.add_ibctype_to_parser(parser)
diff --git a/pyschism/cmd/common.py b/pyschism/cmd/common.py
index 26f211f9..2f84ffb4 100644
--- a/pyschism/cmd/common.py
+++ b/pyschism/cmd/common.py
@@ -259,27 +259,27 @@ def add_tidal_database_to_parser(parser):
     )
 
 
-def add_baroclinic_database_to_parser(parser):
-    class BaroclinicDatabases(Enum):
-        RTOFS = hycom.RTOFS
-        GOFS = hycom.GOFS
-
-    class BaroclinicDatabaseAction(argparse.Action):
-        def __call__(self, parser, namespace, values, option_string=None):
-
-            setattr(namespace, self.dest, values)
-
-    hycom_group = parser.add_argument_group("HYCOM options")
-    # *************************************************************************************************
-    hycom_group.add_argument(
-        "--hycom",
-        # "--baroclinic-database",
-        choices=[x.name.lower() for x in BaroclinicDatabases],
-        dest="baroclinic_database",
-        default="gofs",
-        action=BaroclinicDatabaseAction,
-        help="Options for obtaining HYCOM data. Defaults to GOFS.",
-    )
+# def add_baroclinic_database_to_parser(parser):
+#     class BaroclinicDatabases(Enum):
+#         RTOFS = hycom.RTOFS
+#         GOFS = hycom.GOFS
+# 
+#     class BaroclinicDatabaseAction(argparse.Action):
+#         def __call__(self, parser, namespace, values, option_string=None):
+# 
+#             setattr(namespace, self.dest, values)
+# 
+#     hycom_group = parser.add_argument_group("HYCOM options")
+#     # *************************************************************************************************
+#     hycom_group.add_argument(
+#         "--hycom",
+#         # "--baroclinic-database",
+#         choices=[x.name.lower() for x in BaroclinicDatabases],
+#         dest="baroclinic_database",
+#         default="gofs",
+#         action=BaroclinicDatabaseAction,
+#         help="Options for obtaining HYCOM data. Defaults to GOFS.",
+#     )
 
 
 # class CustomBoundaryAction(argparse.Action):
@@ -312,536 +312,536 @@ def get_per_boundary_help_message(varname, ibctype):
     )
 
 
-def add_iettype_to_parser(parser):
-    class IettypeAction(argparse.Action):
-        def __call__(self, parser, namespace, values, option_string=None):
-            values = json.loads(values)
-            for key, val in values.items():
-                values.update({key: self.iettype(int(key), val)})
-            setattr(namespace, self.dest, values)
-
-        def iettype(self, key, value):
-            values = [1, 2, 3, 4, 5, -1]
-            assert value in values, (
-                "per-boundary iettype values must be " f"one of {values}, not {value}."
-            )
-            raise NotImplementedError(
-                "Per-boundary specification not yet " "supported."
-            )
-
-    class Iettype3Action(argparse.Action):
-        def __call__(self, parser, namespace, values, option_string=None):
-            tmp_parser = argparse.ArgumentParser(add_help=False)
-            add_tidal_constituents_to_parser(tmp_parser)
-            add_tidal_database_to_parser(tmp_parser)
-            tmp_args = tmp_parser.parse_known_args()[0]
-            setattr(
-                namespace,
-                self.dest,
-                self.const(
-                    constituents=tmp_args.constituents, database=tmp_args.tidal_database
-                ),
-            )
-
-    class Iettype4Action(argparse.Action):
-        def __call__(self, parser, namespace, values, option_string=None):
-            if not bool(set(sys.argv).intersection(["-h", "--help"])):
-                if namespace.vgrid.is2D() is True:
-                    raise NotImplementedError("--iettype-4 not available for 2D model.")
-                tmp_parser = argparse.ArgumentParser(add_help=False)
-                add_baroclinic_database_to_parser(tmp_parser)
-                tmp_args = tmp_parser.parse_known_args()[0]
-                setattr(namespace, self.dest, self.const(tmp_args.baroclinic_database))
-
-    class Iettype5Action(argparse.Action):
-        def __call__(self, parser, namespace, values, option_string=None):
-            if not bool(set(sys.argv).intersection(["-h", "--help"])):
-                if namespace.vgrid.is2D() is True:
-                    raise NotImplementedError("--iettype-5 not available for 2D model.")
-                tmp_parser = argparse.ArgumentParser(add_help=False)
-                add_tidal_constituents_to_parser(tmp_parser)
-                add_tidal_database_to_parser(tmp_parser)
-                add_baroclinic_database_to_parser(tmp_parser)
-                tmp_args = tmp_parser.parse_known_args()[0]
-                iettype3 = bctides.iettype.Iettype3(
-                    constituents=tmp_args.constituents, database=tmp_args.tidal_database
-                )
-                iettype4 = bctides.iettype.Iettype4(
-                    data_source=tmp_args.baroclinic_database
-                )
-                setattr(namespace, self.dest, self.const(iettype3, iettype4))
-
-    iettype_group = parser.add_argument_group(
-        "Elevation boundary condition options"
-    ).add_mutually_exclusive_group()
-    iettype_group.add_argument(
-        "--iettype",
-        "--per-boundary-elevation",
-        dest="iettype",
-        metavar="JSON_STYLE_STRING",
-        action=IettypeAction,
-        # const=bctides.iettype.Iettype,
-        help=get_per_boundary_help_message("elevation", "iettype"),
-    )
-    iettype_group.add_argument(
-        "--iettype-1",
-        "--elevation-time-history",
-        dest="iettype",
-        metavar="ELEVATION_TIME_HISTORY",
-        type=bctides.iettype.Iettype1,
-        help="Global elevation options for time history.  (Not implemented)",
-    )
-    iettype_group.add_argument(
-        "--iettype-2",
-        "--constant-elevation-value",
-        dest="iettype",
-        metavar="VALUE",
-        type=bctides.iettype.Iettype2,
-        help="Global constant elevation option.",
-    )
-
-    iettype_group.add_argument(
-        "--iettype-3",
-        "--harmonic-tidal-elevation",
-        dest="iettype",
-        nargs=0,
-        const=bctides.iettype.Iettype3,
-        action=Iettype3Action,
-        help="Enables tidal elevation harmonic forcing.",
-    )
-
-    iettype_group.add_argument(
-        "--iettype-4",
-        "--subtidal-elevation",
-        "--elev-2d",
-        dest="iettype",
-        nargs=0,
-        const=bctides.iettype.Iettype4,
-        action=Iettype4Action,
-        help="Enables subtidal elevation forcing.",
-    )
-
-    iettype_group.add_argument(
-        "--iettype-5",
-        "--subtidal-and-tidal-elevation",
-        dest="iettype",
-        nargs=0,
-        const=bctides.iettype.Iettype5,
-        action=Iettype5Action,
-        help="Enables the combined elevation tidal harmonic forcing and "
-        "subtidal elevation forcing.",
-    )
-
-    iettype_group.add_argument(
-        "--iettype--1",
-        "--zero-elevation",
-        dest="iettype",
-        const=bctides.iettype.Iettype_1,
-        action="store_const",
-        help="Zero elevation at boundaries. (Not implemented)",
-    )
-
-
-def add_ifltype_to_parser(parser):
-    class IfltypeAction(argparse.Action):
-        def __call__(self, parser, namespace, values, option_string=None):
-            values = json.loads(values)
-            for key, val in values.items():
-                values.update({key: self.iettype(int(key), val)})
-            setattr(namespace, self.dest, values)
-
-        def ifltype(self, key, value):
-            values = [1, 2, 3, 4, 5, -1]
-            assert value in values, (
-                "per-boundary ifltype values must be " f"one of {values}, not {value}."
-            )
-            raise NotImplementedError(
-                "Per-boundary specification not yet " "supported."
-            )
-
-    class Ifltype3Action(argparse.Action):
-        def __call__(self, parser, namespace, values, option_string=None):
-            tmp_parser = argparse.ArgumentParser(add_help=False)
-            add_tidal_constituents_to_parser(tmp_parser)
-            add_tidal_database_to_parser(tmp_parser)
-            tmp_args = tmp_parser.parse_known_args()[0]
-            setattr(
-                namespace,
-                self.dest,
-                self.const(
-                    constituents=tmp_args.constituents,
-                    database=tmp_args.tidal_database,
-                ),
-            )
-
-    class Ifltype4Action(argparse.Action):
-        def __call__(self, parser, namespace, values, option_string=None):
-            if not bool(set(sys.argv).intersection(["-h", "--help"])):
-                if namespace.vgrid.is2D() is True:
-                    raise NotImplementedError("--ifltype-4 not available for 2D model.")
-                tmp_parser = argparse.ArgumentParser(add_help=False)
-                add_baroclinic_database_to_parser(tmp_parser)
-                tmp_args = tmp_parser.parse_known_args()[0]
-                setattr(namespace, self.dest, self.const(tmp_args.baroclinic_database))
-
-    class Ifltype5Action(argparse.Action):
-        def __call__(self, parser, namespace, values, option_string=None):
-            if not bool(set(sys.argv).intersection(["-h", "--help"])):
-                if namespace.vgrid.is2D() is True:
-                    raise NotImplementedError("--ifltype-5 not available for 2D model.")
-                tmp_parser = argparse.ArgumentParser(add_help=False)
-                add_tidal_constituents_to_parser(tmp_parser)
-                add_tidal_database_to_parser(tmp_parser)
-                add_baroclinic_database_to_parser(tmp_parser)
-                tmp_args = tmp_parser.parse_known_args()[0]
-                ifltype3 = bctides.ifltype.Ifltype3(
-                    constituents=tmp_args.constituents, database=tmp_args.tidal_database
-                )
-                ifltype4 = bctides.ifltype.Ifltype4(
-                    data_source=tmp_args.baroclinic_database
-                )
-                setattr(namespace, self.dest, self.const(ifltype3, ifltype4))
-
-    ifltype_group = parser.add_argument_group(
-        "Velocity boundary condition options"
-    ).add_mutually_exclusive_group()
-    # ifltype_group.add_argument(
-    #     "--ifltype",
-    #     '--per-boundary-velocity',
-    #     dest="ifltype",
-    #     metavar='JSON_STYLE_STRING',
-    #     action=IfltypeAction,
-    #     help=get_per_boundary_help_message('velocity', 'ifltype')
-    # )
-    ifltype_group.add_argument(
-        "--ifltype-1",
-        "--flux-time-history",
-        metavar="VELOCITY_TIME_HISTORY",
-        dest="ifltype",
-        # type=ifltype.Ifltype1
-        help="Global boundary velocity time history file. (Not implemented)",
-    )
-    ifltype_group.add_argument(
-        "--ifltype-2",
-        "--constant-flux-value",
-        dest="ifltype",
-        metavar="VALUE",
-        # type=ifltype.Ifltype2,
-        help="Global constant flow option.",
-    )
-    ifltype_group.add_argument(
-        "--ifltype-3",
-        "--harmonic-tidal-velocity",
-        dest="ifltype",
-        nargs=0,
-        const=bctides.ifltype.Ifltype3,
-        action=Ifltype3Action,
-        help="Enables tidal velocity harmonic forcing.",
-    )
-
-    ifltype_group.add_argument(
-        "--ifltype-4",
-        "--subtidal-velocity",
-        dest="ifltype",
-        nargs=0,
-        const=bctides.ifltype.Ifltype4,
-        action=Ifltype4Action,
-        help="Enables subtidal velocity forcing.",
-    )
-
-    ifltype_group.add_argument(
-        "--ifltype-5",
-        "--subtidal-and-tidal-velocity",
-        action=Ifltype5Action,
-        nargs=0,
-        dest="ifltype",
-        const=bctides.ifltype.Ifltype5,
-        help="Enables the combined velocity tidal harmonic forcing and "
-        "subtidal velocity forcing.",
-    )
-
-    ifltype_group.add_argument(
-        "--ifltype--1",
-        "--uv-zero",
-        "--flather",
-        # action='store_const',
-        dest="ifltype",
-        # const=ifltype.Ifltype_1,
-    )
-
-
-def add_itetype_to_parser(parser):
-    def get_nudge_bnds(namespace):
-        tmp_parser = argparse.ArgumentParser(add_help=False)
-        add_baroclinic_database_to_parser(tmp_parser)
-        add_nudge_to_parser(tmp_parser)
-        tmp_args = tmp_parser.parse_known_args()[0]
-        bnd_ids = namespace.hgrid.boundaries.open["id"].tolist()
-        if tmp_args.nudge_temp is None:
-            if namespace.vgrid.is2D():
-                tmp_args.nudge_temp = {bnd_id: False for bnd_id in bnd_ids}
-            else:
-                tmp_args.nudge_temp = {bnd_id: True for bnd_id in bnd_ids}
-        elif tmp_args.nudge_temp is True:
-            tmp_args.nudge_temp = {bnd_id: True for bnd_id in bnd_ids}
-        elif tmp_args.nudge_temp is False:
-            tmp_args.nudge_temp = {bnd_id: False for bnd_id in bnd_ids}
-        else:
-            remaining_bounds = list(
-                set(list(tmp_args.nudge_temp.keys())) - set(bnd_ids)
-            )
-            if np.all(list(tmp_args.nudge_temp.values())) is True:
-                for bnd_id in remaining_bounds:
-                    tmp_args.nudge_temp.setdefault(bnd_id, False)
-            else:
-                for bnd_id in remaining_bounds:
-                    tmp_args.nudge_temp.setdefault(bnd_id, True)
-        remaining_bounds = list(set(tmp_args.nudge_temp.keys()) - set(bnd_ids))
-        if len(remaining_bounds) > 0:
-            f = ["No nudging boundaries found with ids:"]
-            for bnd in remaining_bounds:
-                f.append(bnd)
-            raise ValueError(" ".join(f))
-        return tmp_args.nudge_temp
-
-    class Itetype4Action(argparse.Action):
-        def __call__(self, parser, namespace, values, option_string=None):
-            tmp_parser = argparse.ArgumentParser(add_help=False)
-            add_nudge_to_parser(tmp_parser)
-            tmp_args = tmp_parser.parse_known_args()[0]
-            if namespace.hgrid is not None:
-                nudge_bounds = get_nudge_bnds(namespace)
-                itetype4 = bctides.itetype.Itetype4(
-                    nudge=nudge_bounds["1"],
-                    rlmax=tmp_args.rlmax_temp,
-                    rnu_day=tmp_args.rnu_day_temp,
-                )
-                setattr(namespace, self.dest, itetype4)
-
-    itetype_group = parser.add_argument_group(
-        "Temperature boundary condition options"
-    ).add_mutually_exclusive_group()
-    # itetype_group.add_argument(
-    #     "--itetype",
-    #     dest="itetype",
-    #     # nargs=1,  # 0 or more
-    #     # const={},
-    #     # action=CustomBoundaryAction,
-    # )
-    itetype_group.add_argument(
-        "--itetype-1",
-        "--temperature-time-history",
-        dest="itetype",
-        # type=itetype.Itetype1
-    )
-    itetype_group.add_argument(
-        "--itetype-2",
-        "--constant-temperature-value",
-        dest="itetype",
-        # type=itetype.Itetype2
-    )
-
-    itetype_group.add_argument(
-        "--itetype-3",
-        "--temperature-initial-condition",
-        dest="itetype",
-        # type=itetype.Itetype3,
-    )
-
-    itetype_group.add_argument(
-        "--itetype-4",
-        "--temp-3d",
-        action=Itetype4Action,
-        dest="itetype",
-        nargs=0,
-        const=bctides.itetype.Itetype4,
-    )
-
-
-def add_isatype_to_parser(parser):
-    def get_nudge_bnds(namespace):
-        tmp_parser = argparse.ArgumentParser(add_help=False)
-        add_baroclinic_database_to_parser(tmp_parser)
-        add_nudge_to_parser(tmp_parser)
-        tmp_args = tmp_parser.parse_known_args()[0]
-        bnd_ids = namespace.hgrid.boundaries.open["id"].tolist()
-        if tmp_args.nudge_salt is None:
-            if namespace.vgrid.is2D():
-                tmp_args.nudge_salt = {bnd_id: False for bnd_id in bnd_ids}
-            else:
-                tmp_args.nudge_salt = {bnd_id: True for bnd_id in bnd_ids}
-        elif tmp_args.nudge_salt is True:
-            tmp_args.nudge_salt = {bnd_id: True for bnd_id in bnd_ids}
-        elif tmp_args.nudge_salt is False:
-            tmp_args.nudge_salt = {bnd_id: False for bnd_id in bnd_ids}
-        else:
-            remaining_bounds = list(
-                set(list(tmp_args.nudge_salt.keys())) - set(bnd_ids)
-            )
-            if np.all(list(tmp_args.nudge_salt.values())) is True:
-                for bnd_id in remaining_bounds:
-                    tmp_args.nudge_salt.setdefault(bnd_id, False)
-            else:
-                for bnd_id in remaining_bounds:
-                    tmp_args.nudge_salt.setdefault(bnd_id, True)
-        remaining_bounds = list(set(tmp_args.nudge_salt.keys()) - set(bnd_ids))
-        if len(remaining_bounds) > 0:
-            f = ["No nudging boundaries found with ids:"]
-            for bnd in remaining_bounds:
-                f.append(bnd)
-            raise ValueError(" ".join(f))
-        return tmp_args.nudge_salt
-
-    class Isatype4Action(argparse.Action):
-        def __call__(self, parser, namespace, values, option_string=None):
-            tmp_parser = argparse.ArgumentParser(add_help=False)
-            add_nudge_to_parser(tmp_parser)
-            tmp_args = tmp_parser.parse_known_args()[0]
-            if namespace.hgrid is not None:
-                nudge_bounds = get_nudge_bnds(namespace)
-                isatype4 = bctides.isatype.Isatype4(
-                    nudge=nudge_bounds["1"],
-                    rlmax=tmp_args.rlmax_salt,
-                    rnu_day=tmp_args.rnu_day_salt,
-                )
-                setattr(namespace, self.dest, isatype4)
-
-    isatype_group = parser.add_argument_group(
-        "Salinity boundary condition options"
-    ).add_mutually_exclusive_group()
-    # isatype_group.add_argument(
-    #     "--isatype",
-    #     dest="isatype",
-    #     # nargs=1,  # 0 or more
-    #     # const={},
-    #     # action=CustomBoundaryAction,
-    # )
-    isatype_group.add_argument(
-        "--isatype-1",
-        "--salt-th",
-        dest="isatype",
-        # type=isatype.Isatype1
-    )
-    isatype_group.add_argument(
-        "--isatype-2",
-        "--salt-val",
-        dest="isatype",
-        # type=isatype.Isatype2
-    )
-
-    isatype_group.add_argument(
-        "--isatype-3",
-        # "--salt-ic",
-        dest="isatype",
-        # type=isatype.Isatype3,
-    )
-
-    isatype_group.add_argument(
-        "--isatype-4",
-        "--salt-3d",
-        action=Isatype4Action,
-        dest="isatype",
-        nargs=0,
-        const=bctides.isatype.Isatype4,
-    )
-
-
-class NudgeAction(argparse.Action):
-    def __call__(self, parser, namespace, values, option_string=None):
-        nudge = {}
-        if "disable" in option_string:
-            for disable_bnd in values:
-                nudge.setdefault(disable_bnd, False)
-            if len(nudge) == 0:
-                nudge = False
-        else:
-            for enable_bnd in values:
-                nudge.setdefault(enable_bnd, True)
-            if len(nudge) == 0:
-                nudge = True
-        setattr(namespace, self.dest, nudge)
-
-
-def add_temperature_nudge_to_parser(parser):
-    temp_group = parser.add_argument_group("Nudge options for temperature")
-    nudge_parser = temp_group.add_mutually_exclusive_group()
-    nudge_parser.add_argument(
-        "--nudge-temp",
-        "--nudge-temperature",
-        "--nudge-t",
-        dest="nudge_temp",
-        nargs="*",
-        action=NudgeAction,
-        help="Enables temperature nudging. If no arguments are given, the nudging "
-        "is enabled for all the boundaries. If a per-boundary nudging is desired, "
-        "you can pass the boundary id's of the boundaries where nudging should "
-        "be enabled, for example passing --nudge-temp 1 2 4 will enable nudging "
-        "for boundaries with id's 1, 2 and 4, and will disable the nudging "
-        "for the remaining boundaries. Note: This option is mutually exclusive "
-        "to --disable-nudge-temperature.",
-    )
-    nudge_parser.add_argument(
-        "--disable-nudge-temp",
-        "--disable-nudge-temperature",
-        "--disable-nudge-t",
-        dest="nudge_temp",
-        nargs="*",
-        action=NudgeAction,
-        help="Disables temperature nudging. If no arguments are given, the nudging "
-        "is disabled for all the boundaries. If per-boundary nudging disabling is desired, "
-        "you can pass the boundary id's of the boundaries where nudging should "
-        "be disabled, for example passing --disable-nudge-temp 1 2 4 will disable nudging "
-        "for boundaries with id's 1, 2 and 4, and will enable the nudging "
-        "for the remaining boundaries. Note: This option is mutually exclusive "
-        "to --nudge-temperature.",
-    )
-    temp_group.set_defaults(nudge_temp=None)
-    temp_group.add_argument("--rlmax-temp", default=1.5, type=float)
-    temp_group.add_argument("--rnu_day-temp", default=0.25, type=float)
-
-
-def add_salinity_nudge_to_parser(parser):
-    salt_group = parser.add_argument_group("Nudge options for salinity")
-    salt_group.add_argument(
-        "--nudge-salt",
-        "--nudge-salinity",
-        "--nudge-s",
-        dest="nudge_salt",
-        nargs="?",
-        action=NudgeAction,
-    )
-    salt_group.add_argument(
-        "--disable-nudge-salt",
-        "--disable-nudge-salinity",
-        "--disable-nudge-s",
-        dest="nudge_salt",
-        nargs="?",
-        action=NudgeAction,
-    )
-    salt_group.set_defaults(nudge_salt=None)
-    salt_group.add_argument("--rlmax-salt", default=1.5, type=float)
-    salt_group.add_argument("--rnu_day-salt", default=0.25, type=float)
-
-
-def add_nudge_to_parser(parser):
-    add_temperature_nudge_to_parser(parser)
-    add_salinity_nudge_to_parser(parser)
-
-
-def add_ibctype_to_parser(parser):
-    add_iettype_to_parser(parser)
-    add_ifltype_to_parser(parser)
-    add_itetype_to_parser(parser)
-    add_isatype_to_parser(parser)
-    # add_itrtype_to_parser(parser)
-    add_nudge_to_parser(parser)
-
-
-def add_windrot_to_parser(parser):
-    parser.add_argument(
-        "--windrot",
-        # action=WindrotAction,
-    )
+# def add_iettype_to_parser(parser):
+#     class IettypeAction(argparse.Action):
+#         def __call__(self, parser, namespace, values, option_string=None):
+#             values = json.loads(values)
+#             for key, val in values.items():
+#                 values.update({key: self.iettype(int(key), val)})
+#             setattr(namespace, self.dest, values)
+# 
+#         def iettype(self, key, value):
+#             values = [1, 2, 3, 4, 5, -1]
+#             assert value in values, (
+#                 "per-boundary iettype values must be " f"one of {values}, not {value}."
+#             )
+#             raise NotImplementedError(
+#                 "Per-boundary specification not yet " "supported."
+#             )
+# 
+#     class Iettype3Action(argparse.Action):
+#         def __call__(self, parser, namespace, values, option_string=None):
+#             tmp_parser = argparse.ArgumentParser(add_help=False)
+#             add_tidal_constituents_to_parser(tmp_parser)
+#             add_tidal_database_to_parser(tmp_parser)
+#             tmp_args = tmp_parser.parse_known_args()[0]
+#             setattr(
+#                 namespace,
+#                 self.dest,
+#                 self.const(
+#                     constituents=tmp_args.constituents, database=tmp_args.tidal_database
+#                 ),
+#             )
+# 
+#     class Iettype4Action(argparse.Action):
+#         def __call__(self, parser, namespace, values, option_string=None):
+#             if not bool(set(sys.argv).intersection(["-h", "--help"])):
+#                 if namespace.vgrid.is2D() is True:
+#                     raise NotImplementedError("--iettype-4 not available for 2D model.")
+#                 tmp_parser = argparse.ArgumentParser(add_help=False)
+#                 add_baroclinic_database_to_parser(tmp_parser)
+#                 tmp_args = tmp_parser.parse_known_args()[0]
+#                 setattr(namespace, self.dest, self.const(tmp_args.baroclinic_database))
+# 
+#     class Iettype5Action(argparse.Action):
+#         def __call__(self, parser, namespace, values, option_string=None):
+#             if not bool(set(sys.argv).intersection(["-h", "--help"])):
+#                 if namespace.vgrid.is2D() is True:
+#                     raise NotImplementedError("--iettype-5 not available for 2D model.")
+#                 tmp_parser = argparse.ArgumentParser(add_help=False)
+#                 add_tidal_constituents_to_parser(tmp_parser)
+#                 add_tidal_database_to_parser(tmp_parser)
+#                 add_baroclinic_database_to_parser(tmp_parser)
+#                 tmp_args = tmp_parser.parse_known_args()[0]
+#                 iettype3 = bctides.iettype.Iettype3(
+#                     constituents=tmp_args.constituents, database=tmp_args.tidal_database
+#                 )
+#                 iettype4 = bctides.iettype.Iettype4(
+#                     data_source=tmp_args.baroclinic_database
+#                 )
+#                 setattr(namespace, self.dest, self.const(iettype3, iettype4))
+# 
+#     iettype_group = parser.add_argument_group(
+#         "Elevation boundary condition options"
+#     ).add_mutually_exclusive_group()
+#     iettype_group.add_argument(
+#         "--iettype",
+#         "--per-boundary-elevation",
+#         dest="iettype",
+#         metavar="JSON_STYLE_STRING",
+#         action=IettypeAction,
+#         # const=bctides.iettype.Iettype,
+#         help=get_per_boundary_help_message("elevation", "iettype"),
+#     )
+#     iettype_group.add_argument(
+#         "--iettype-1",
+#         "--elevation-time-history",
+#         dest="iettype",
+#         metavar="ELEVATION_TIME_HISTORY",
+#         type=bctides.iettype.Iettype1,
+#         help="Global elevation options for time history.  (Not implemented)",
+#     )
+#     iettype_group.add_argument(
+#         "--iettype-2",
+#         "--constant-elevation-value",
+#         dest="iettype",
+#         metavar="VALUE",
+#         type=bctides.iettype.Iettype2,
+#         help="Global constant elevation option.",
+#     )
+# 
+#     iettype_group.add_argument(
+#         "--iettype-3",
+#         "--harmonic-tidal-elevation",
+#         dest="iettype",
+#         nargs=0,
+#         const=bctides.iettype.Iettype3,
+#         action=Iettype3Action,
+#         help="Enables tidal elevation harmonic forcing.",
+#     )
+# 
+#     iettype_group.add_argument(
+#         "--iettype-4",
+#         "--subtidal-elevation",
+#         "--elev-2d",
+#         dest="iettype",
+#         nargs=0,
+#         const=bctides.iettype.Iettype4,
+#         action=Iettype4Action,
+#         help="Enables subtidal elevation forcing.",
+#     )
+# 
+#     iettype_group.add_argument(
+#         "--iettype-5",
+#         "--subtidal-and-tidal-elevation",
+#         dest="iettype",
+#         nargs=0,
+#         const=bctides.iettype.Iettype5,
+#         action=Iettype5Action,
+#         help="Enables the combined elevation tidal harmonic forcing and "
+#         "subtidal elevation forcing.",
+#     )
+# 
+#     iettype_group.add_argument(
+#         "--iettype--1",
+#         "--zero-elevation",
+#         dest="iettype",
+#         const=bctides.iettype.Iettype_1,
+#         action="store_const",
+#         help="Zero elevation at boundaries. (Not implemented)",
+#     )
+# 
+# 
+# def add_ifltype_to_parser(parser):
+#     class IfltypeAction(argparse.Action):
+#         def __call__(self, parser, namespace, values, option_string=None):
+#             values = json.loads(values)
+#             for key, val in values.items():
+#                 values.update({key: self.iettype(int(key), val)})
+#             setattr(namespace, self.dest, values)
+# 
+#         def ifltype(self, key, value):
+#             values = [1, 2, 3, 4, 5, -1]
+#             assert value in values, (
+#                 "per-boundary ifltype values must be " f"one of {values}, not {value}."
+#             )
+#             raise NotImplementedError(
+#                 "Per-boundary specification not yet " "supported."
+#             )
+# 
+#     class Ifltype3Action(argparse.Action):
+#         def __call__(self, parser, namespace, values, option_string=None):
+#             tmp_parser = argparse.ArgumentParser(add_help=False)
+#             add_tidal_constituents_to_parser(tmp_parser)
+#             add_tidal_database_to_parser(tmp_parser)
+#             tmp_args = tmp_parser.parse_known_args()[0]
+#             setattr(
+#                 namespace,
+#                 self.dest,
+#                 self.const(
+#                     constituents=tmp_args.constituents,
+#                     database=tmp_args.tidal_database,
+#                 ),
+#             )
+# 
+#     class Ifltype4Action(argparse.Action):
+#         def __call__(self, parser, namespace, values, option_string=None):
+#             if not bool(set(sys.argv).intersection(["-h", "--help"])):
+#                 if namespace.vgrid.is2D() is True:
+#                     raise NotImplementedError("--ifltype-4 not available for 2D model.")
+#                 tmp_parser = argparse.ArgumentParser(add_help=False)
+#                 add_baroclinic_database_to_parser(tmp_parser)
+#                 tmp_args = tmp_parser.parse_known_args()[0]
+#                 setattr(namespace, self.dest, self.const(tmp_args.baroclinic_database))
+# 
+#     class Ifltype5Action(argparse.Action):
+#         def __call__(self, parser, namespace, values, option_string=None):
+#             if not bool(set(sys.argv).intersection(["-h", "--help"])):
+#                 if namespace.vgrid.is2D() is True:
+#                     raise NotImplementedError("--ifltype-5 not available for 2D model.")
+#                 tmp_parser = argparse.ArgumentParser(add_help=False)
+#                 add_tidal_constituents_to_parser(tmp_parser)
+#                 add_tidal_database_to_parser(tmp_parser)
+#                 add_baroclinic_database_to_parser(tmp_parser)
+#                 tmp_args = tmp_parser.parse_known_args()[0]
+#                 ifltype3 = bctides.ifltype.Ifltype3(
+#                     constituents=tmp_args.constituents, database=tmp_args.tidal_database
+#                 )
+#                 ifltype4 = bctides.ifltype.Ifltype4(
+#                     data_source=tmp_args.baroclinic_database
+#                 )
+#                 setattr(namespace, self.dest, self.const(ifltype3, ifltype4))
+# 
+#     ifltype_group = parser.add_argument_group(
+#         "Velocity boundary condition options"
+#     ).add_mutually_exclusive_group()
+#     # ifltype_group.add_argument(
+#     #     "--ifltype",
+#     #     '--per-boundary-velocity',
+#     #     dest="ifltype",
+#     #     metavar='JSON_STYLE_STRING',
+#     #     action=IfltypeAction,
+#     #     help=get_per_boundary_help_message('velocity', 'ifltype')
+#     # )
+#     ifltype_group.add_argument(
+#         "--ifltype-1",
+#         "--flux-time-history",
+#         metavar="VELOCITY_TIME_HISTORY",
+#         dest="ifltype",
+#         # type=ifltype.Ifltype1
+#         help="Global boundary velocity time history file. (Not implemented)",
+#     )
+#     ifltype_group.add_argument(
+#         "--ifltype-2",
+#         "--constant-flux-value",
+#         dest="ifltype",
+#         metavar="VALUE",
+#         # type=ifltype.Ifltype2,
+#         help="Global constant flow option.",
+#     )
+#     ifltype_group.add_argument(
+#         "--ifltype-3",
+#         "--harmonic-tidal-velocity",
+#         dest="ifltype",
+#         nargs=0,
+#         const=bctides.ifltype.Ifltype3,
+#         action=Ifltype3Action,
+#         help="Enables tidal velocity harmonic forcing.",
+#     )
+# 
+#     ifltype_group.add_argument(
+#         "--ifltype-4",
+#         "--subtidal-velocity",
+#         dest="ifltype",
+#         nargs=0,
+#         const=bctides.ifltype.Ifltype4,
+#         action=Ifltype4Action,
+#         help="Enables subtidal velocity forcing.",
+#     )
+# 
+#     ifltype_group.add_argument(
+#         "--ifltype-5",
+#         "--subtidal-and-tidal-velocity",
+#         action=Ifltype5Action,
+#         nargs=0,
+#         dest="ifltype",
+#         const=bctides.ifltype.Ifltype5,
+#         help="Enables the combined velocity tidal harmonic forcing and "
+#         "subtidal velocity forcing.",
+#     )
+# 
+#     ifltype_group.add_argument(
+#         "--ifltype--1",
+#         "--uv-zero",
+#         "--flather",
+#         # action='store_const',
+#         dest="ifltype",
+#         # const=ifltype.Ifltype_1,
+#     )
+# 
+# 
+# def add_itetype_to_parser(parser):
+#     def get_nudge_bnds(namespace):
+#         tmp_parser = argparse.ArgumentParser(add_help=False)
+#         add_baroclinic_database_to_parser(tmp_parser)
+#         add_nudge_to_parser(tmp_parser)
+#         tmp_args = tmp_parser.parse_known_args()[0]
+#         bnd_ids = namespace.hgrid.boundaries.open["id"].tolist()
+#         if tmp_args.nudge_temp is None:
+#             if namespace.vgrid.is2D():
+#                 tmp_args.nudge_temp = {bnd_id: False for bnd_id in bnd_ids}
+#             else:
+#                 tmp_args.nudge_temp = {bnd_id: True for bnd_id in bnd_ids}
+#         elif tmp_args.nudge_temp is True:
+#             tmp_args.nudge_temp = {bnd_id: True for bnd_id in bnd_ids}
+#         elif tmp_args.nudge_temp is False:
+#             tmp_args.nudge_temp = {bnd_id: False for bnd_id in bnd_ids}
+#         else:
+#             remaining_bounds = list(
+#                 set(list(tmp_args.nudge_temp.keys())) - set(bnd_ids)
+#             )
+#             if np.all(list(tmp_args.nudge_temp.values())) is True:
+#                 for bnd_id in remaining_bounds:
+#                     tmp_args.nudge_temp.setdefault(bnd_id, False)
+#             else:
+#                 for bnd_id in remaining_bounds:
+#                     tmp_args.nudge_temp.setdefault(bnd_id, True)
+#         remaining_bounds = list(set(tmp_args.nudge_temp.keys()) - set(bnd_ids))
+#         if len(remaining_bounds) > 0:
+#             f = ["No nudging boundaries found with ids:"]
+#             for bnd in remaining_bounds:
+#                 f.append(bnd)
+#             raise ValueError(" ".join(f))
+#         return tmp_args.nudge_temp
+# 
+#     class Itetype4Action(argparse.Action):
+#         def __call__(self, parser, namespace, values, option_string=None):
+#             tmp_parser = argparse.ArgumentParser(add_help=False)
+#             add_nudge_to_parser(tmp_parser)
+#             tmp_args = tmp_parser.parse_known_args()[0]
+#             if namespace.hgrid is not None:
+#                 nudge_bounds = get_nudge_bnds(namespace)
+#                 itetype4 = bctides.itetype.Itetype4(
+#                     nudge=nudge_bounds["1"],
+#                     rlmax=tmp_args.rlmax_temp,
+#                     rnu_day=tmp_args.rnu_day_temp,
+#                 )
+#                 setattr(namespace, self.dest, itetype4)
+# 
+#     itetype_group = parser.add_argument_group(
+#         "Temperature boundary condition options"
+#     ).add_mutually_exclusive_group()
+#     # itetype_group.add_argument(
+#     #     "--itetype",
+#     #     dest="itetype",
+#     #     # nargs=1,  # 0 or more
+#     #     # const={},
+#     #     # action=CustomBoundaryAction,
+#     # )
+#     itetype_group.add_argument(
+#         "--itetype-1",
+#         "--temperature-time-history",
+#         dest="itetype",
+#         # type=itetype.Itetype1
+#     )
+#     itetype_group.add_argument(
+#         "--itetype-2",
+#         "--constant-temperature-value",
+#         dest="itetype",
+#         # type=itetype.Itetype2
+#     )
+# 
+#     itetype_group.add_argument(
+#         "--itetype-3",
+#         "--temperature-initial-condition",
+#         dest="itetype",
+#         # type=itetype.Itetype3,
+#     )
+# 
+#     itetype_group.add_argument(
+#         "--itetype-4",
+#         "--temp-3d",
+#         action=Itetype4Action,
+#         dest="itetype",
+#         nargs=0,
+#         const=bctides.itetype.Itetype4,
+#     )
+# 
+# 
+# def add_isatype_to_parser(parser):
+#     def get_nudge_bnds(namespace):
+#         tmp_parser = argparse.ArgumentParser(add_help=False)
+#         add_baroclinic_database_to_parser(tmp_parser)
+#         add_nudge_to_parser(tmp_parser)
+#         tmp_args = tmp_parser.parse_known_args()[0]
+#         bnd_ids = namespace.hgrid.boundaries.open["id"].tolist()
+#         if tmp_args.nudge_salt is None:
+#             if namespace.vgrid.is2D():
+#                 tmp_args.nudge_salt = {bnd_id: False for bnd_id in bnd_ids}
+#             else:
+#                 tmp_args.nudge_salt = {bnd_id: True for bnd_id in bnd_ids}
+#         elif tmp_args.nudge_salt is True:
+#             tmp_args.nudge_salt = {bnd_id: True for bnd_id in bnd_ids}
+#         elif tmp_args.nudge_salt is False:
+#             tmp_args.nudge_salt = {bnd_id: False for bnd_id in bnd_ids}
+#         else:
+#             remaining_bounds = list(
+#                 set(list(tmp_args.nudge_salt.keys())) - set(bnd_ids)
+#             )
+#             if np.all(list(tmp_args.nudge_salt.values())) is True:
+#                 for bnd_id in remaining_bounds:
+#                     tmp_args.nudge_salt.setdefault(bnd_id, False)
+#             else:
+#                 for bnd_id in remaining_bounds:
+#                     tmp_args.nudge_salt.setdefault(bnd_id, True)
+#         remaining_bounds = list(set(tmp_args.nudge_salt.keys()) - set(bnd_ids))
+#         if len(remaining_bounds) > 0:
+#             f = ["No nudging boundaries found with ids:"]
+#             for bnd in remaining_bounds:
+#                 f.append(bnd)
+#             raise ValueError(" ".join(f))
+#         return tmp_args.nudge_salt
+# 
+#     class Isatype4Action(argparse.Action):
+#         def __call__(self, parser, namespace, values, option_string=None):
+#             tmp_parser = argparse.ArgumentParser(add_help=False)
+#             add_nudge_to_parser(tmp_parser)
+#             tmp_args = tmp_parser.parse_known_args()[0]
+#             if namespace.hgrid is not None:
+#                 nudge_bounds = get_nudge_bnds(namespace)
+#                 isatype4 = bctides.isatype.Isatype4(
+#                     nudge=nudge_bounds["1"],
+#                     rlmax=tmp_args.rlmax_salt,
+#                     rnu_day=tmp_args.rnu_day_salt,
+#                 )
+#                 setattr(namespace, self.dest, isatype4)
+# 
+#     isatype_group = parser.add_argument_group(
+#         "Salinity boundary condition options"
+#     ).add_mutually_exclusive_group()
+#     # isatype_group.add_argument(
+#     #     "--isatype",
+#     #     dest="isatype",
+#     #     # nargs=1,  # 0 or more
+#     #     # const={},
+#     #     # action=CustomBoundaryAction,
+#     # )
+#     isatype_group.add_argument(
+#         "--isatype-1",
+#         "--salt-th",
+#         dest="isatype",
+#         # type=isatype.Isatype1
+#     )
+#     isatype_group.add_argument(
+#         "--isatype-2",
+#         "--salt-val",
+#         dest="isatype",
+#         # type=isatype.Isatype2
+#     )
+# 
+#     isatype_group.add_argument(
+#         "--isatype-3",
+#         # "--salt-ic",
+#         dest="isatype",
+#         # type=isatype.Isatype3,
+#     )
+# 
+#     isatype_group.add_argument(
+#         "--isatype-4",
+#         "--salt-3d",
+#         action=Isatype4Action,
+#         dest="isatype",
+#         nargs=0,
+#         const=bctides.isatype.Isatype4,
+#     )
+# 
+# 
+# class NudgeAction(argparse.Action):
+#     def __call__(self, parser, namespace, values, option_string=None):
+#         nudge = {}
+#         if "disable" in option_string:
+#             for disable_bnd in values:
+#                 nudge.setdefault(disable_bnd, False)
+#             if len(nudge) == 0:
+#                 nudge = False
+#         else:
+#             for enable_bnd in values:
+#                 nudge.setdefault(enable_bnd, True)
+#             if len(nudge) == 0:
+#                 nudge = True
+#         setattr(namespace, self.dest, nudge)
+# 
+# 
+# def add_temperature_nudge_to_parser(parser):
+#     temp_group = parser.add_argument_group("Nudge options for temperature")
+#     nudge_parser = temp_group.add_mutually_exclusive_group()
+#     nudge_parser.add_argument(
+#         "--nudge-temp",
+#         "--nudge-temperature",
+#         "--nudge-t",
+#         dest="nudge_temp",
+#         nargs="*",
+#         action=NudgeAction,
+#         help="Enables temperature nudging. If no arguments are given, the nudging "
+#         "is enabled for all the boundaries. If a per-boundary nudging is desired, "
+#         "you can pass the boundary id's of the boundaries where nudging should "
+#         "be enabled, for example passing --nudge-temp 1 2 4 will enable nudging "
+#         "for boundaries with id's 1, 2 and 4, and will disable the nudging "
+#         "for the remaining boundaries. Note: This option is mutually exclusive "
+#         "to --disable-nudge-temperature.",
+#     )
+#     nudge_parser.add_argument(
+#         "--disable-nudge-temp",
+#         "--disable-nudge-temperature",
+#         "--disable-nudge-t",
+#         dest="nudge_temp",
+#         nargs="*",
+#         action=NudgeAction,
+#         help="Disables temperature nudging. If no arguments are given, the nudging "
+#         "is disabled for all the boundaries. If per-boundary nudging disabling is desired, "
+#         "you can pass the boundary id's of the boundaries where nudging should "
+#         "be disabled, for example passing --disable-nudge-temp 1 2 4 will disable nudging "
+#         "for boundaries with id's 1, 2 and 4, and will enable the nudging "
+#         "for the remaining boundaries. Note: This option is mutually exclusive "
+#         "to --nudge-temperature.",
+#     )
+#     temp_group.set_defaults(nudge_temp=None)
+#     temp_group.add_argument("--rlmax-temp", default=1.5, type=float)
+#     temp_group.add_argument("--rnu_day-temp", default=0.25, type=float)
+# 
+# 
+# def add_salinity_nudge_to_parser(parser):
+#     salt_group = parser.add_argument_group("Nudge options for salinity")
+#     salt_group.add_argument(
+#         "--nudge-salt",
+#         "--nudge-salinity",
+#         "--nudge-s",
+#         dest="nudge_salt",
+#         nargs="?",
+#         action=NudgeAction,
+#     )
+#     salt_group.add_argument(
+#         "--disable-nudge-salt",
+#         "--disable-nudge-salinity",
+#         "--disable-nudge-s",
+#         dest="nudge_salt",
+#         nargs="?",
+#         action=NudgeAction,
+#     )
+#     salt_group.set_defaults(nudge_salt=None)
+#     salt_group.add_argument("--rlmax-salt", default=1.5, type=float)
+#     salt_group.add_argument("--rnu_day-salt", default=0.25, type=float)
+# 
+# 
+# def add_nudge_to_parser(parser):
+#     add_temperature_nudge_to_parser(parser)
+#     add_salinity_nudge_to_parser(parser)
+# 
+# 
+# def add_ibctype_to_parser(parser):
+#     add_iettype_to_parser(parser)
+#     add_ifltype_to_parser(parser)
+#     add_itetype_to_parser(parser)
+#     add_isatype_to_parser(parser)
+#     # add_itrtype_to_parser(parser)
+#     add_nudge_to_parser(parser)
+# 
+# 
+# def add_windrot_to_parser(parser):
+#     parser.add_argument(
+#         "--windrot",
+#         # action=WindrotAction,
+#     )
 
 
 def add_sflux_options_to_parser(parser):
@@ -892,7 +892,7 @@ def add_sflux_to_parser(parser):
         metavar="sflux_level",
     )
     add_sflux_options_to_parser(parser)
-    add_windrot_to_parser(parser)
+    #add_windrot_to_parser(parser)
 
 
 def add_nws_to_parser(parser):
@@ -912,7 +912,7 @@ def add_nws_to_parser(parser):
         metavar="sflux_level",
     )
     add_sflux_options_to_parser(parser)
-    add_windrot_to_parser(parser)
+    #add_windrot_to_parser(parser)
     nws_parser.add_argument(
         "--nws-3",
         dest="nws",
diff --git a/pyschism/cmd/forecast.py b/pyschism/cmd/forecast.py
index 4dccdcc4..67cc2622 100644
--- a/pyschism/cmd/forecast.py
+++ b/pyschism/cmd/forecast.py
@@ -302,10 +302,10 @@ def add_forecast_init_to_parser(parser):
     common.add_gridgr3_to_parser(parser)
     common.add_ic_to_parser(parser)
     common.add_prop_to_parser(parser)
-    common.add_ibctype_to_parser(parser)
+    #common.add_ibctype_to_parser(parser)
     common.add_tidal_constituents_to_parser(parser)
     common.add_tidal_database_to_parser(parser)
-    common.add_baroclinic_database_to_parser(parser)
+    #common.add_baroclinic_database_to_parser(parser)
     common.add_bctides_options_to_parser(parser)
     common.add_nws_to_parser(parser)
     common.add_source_sink_to_parser(parser)
diff --git a/pyschism/forcing/nws/nws2/gfs2.py b/pyschism/forcing/nws/nws2/gfs2.py
index a57a5d86..5267fcd1 100644
--- a/pyschism/forcing/nws/nws2/gfs2.py
+++ b/pyschism/forcing/nws/nws2/gfs2.py
@@ -87,21 +87,26 @@ def files(self):
         return grbfiles
 
 class GFS:
-    def __init__(self, start_date=None, rnday=None, pscr=None, record=1, bbox=None):
-
+    def __init__(self, start_date=None, rnday=None, pscr=None, record=1, bbox=None, outdir=None): 
         start_date = nearest_cycle() if start_date is None else start_date 
         logger.info(f'start_date is {start_date}')
         self.bbox = bbox
 
+        if outdir is None:
+            self.path = pathlib.Path(start_date.strftime("%Y%m%d"))
+            self.path.mkdir(parents=True, exist_ok=True)
+        else:
+            self.path = outdir
+
         end_date = start_date + timedelta(days=rnday)
         
         datevector = np.arange(start_date, end_date, np.timedelta64(1, 'D'),
                 dtype='datetime64')
         datevector = pd.to_datetime(datevector)
 
-        npool = len(datevector) if len(datevector) < mp.cpu_count() else mp.cpu_count()
+        npool = len(datevector) if len(datevector) < mp.cpu_count()/2 else mp.cpu_count()/2
         logger.info(f'npool is {npool}')
-        pool = mp.Pool(npool)
+        pool = mp.Pool(int(npool))
 
         pool.starmap(self.gen_sflux, [(date, record, pscr) for date in datevector])
 
@@ -113,8 +118,6 @@ def gen_sflux(self, date, record, pscr):
         grbfiles = inventory.files
         cycle = date.hour
 
-        path = pathlib.Path(date.strftime("%Y%m%d"))
-        path.mkdir(parents=True, exist_ok=True)
         
         prate = []
         dlwrf = []
@@ -258,7 +261,7 @@ def gen_sflux(self, date, record, pscr):
             'long_name': 'Downward long-wave radiation flux'
         }
                          
-        fout.to_netcdf(path / f'gfs_{date.strftime("%Y%m%d")}{cycle:02d}.nc','w', 'NETCDF3_CLASSIC', unlimited_dims='time')
+        fout.to_netcdf(self.path / f'gfs_{date.strftime("%Y%m%d")}{cycle:02d}.nc','w', 'NETCDF3_CLASSIC', unlimited_dims='time')
 
 
     def modified_latlon(self, grbfile):
diff --git a/pyschism/forcing/nws/nws2/hrrr3.py b/pyschism/forcing/nws/nws2/hrrr3.py
index 28b64f5c..cd0004af 100644
--- a/pyschism/forcing/nws/nws2/hrrr3.py
+++ b/pyschism/forcing/nws/nws2/hrrr3.py
@@ -95,11 +95,17 @@ def files(self):
 
 class HRRR:
 
-    def __init__(self, start_date=None, rnday=None, pscr=None, record=2, bbox=None):
+    def __init__(self, start_date=None, rnday=None, pscr=None, record=2, bbox=None, outdir=None):
 
         start_date = nearest_cycle() if start_date is None else start_date 
         self.bbox = bbox
 
+        if outdir is None:
+            self.path = pathlib.Path(start_date.strftime("%Y%m%d"))
+            self.path.mkdir(parents=True, exist_ok=True)
+        else:
+            self.path = outdir
+
         end_date = start_date + timedelta(days=rnday)
         logger.info(f'start_date is {start_date}, end_date is {end_date}')
         
@@ -107,9 +113,10 @@ def __init__(self, start_date=None, rnday=None, pscr=None, record=2, bbox=None):
                 dtype='datetime64')
         datevector = pd.to_datetime(datevector)
 
-        npool = len(datevector) if len(datevector) < mp.cpu_count() else mp.cpu_count()
+        #Use all CPUs may cause memory issue
+        npool = len(datevector) if len(datevector) < mp.cpu_count()/2 else mp.cpu_count()/2
         logger.info(f'npool is {npool}')
-        pool = mp.Pool(npool)
+        pool = mp.Pool(int(npool))
 
         pool.starmap(self.gen_sflux, [(date, record, pscr) for date in datevector])
 
@@ -121,9 +128,6 @@ def gen_sflux(self, date, record, pscr):
         grbfiles = inventory.files
         cycle = date.hour
 
-        path = pathlib.Path(date.strftime("%Y%m%d"))
-        path.mkdir(parents=True, exist_ok=True)
-
         stmp = [] 
         spfh = []
         uwind = []
@@ -251,7 +255,7 @@ def gen_sflux(self, date, record, pscr):
             'long_name': 'Downward long-wave radiation flux'
         }
                          
-        fout.to_netcdf(path / f'hrrr_{date.strftime("%Y%m%d")}{cycle:02d}.nc','w', 'NETCDF3_CLASSIC', unlimited_dims='time')
+        fout.to_netcdf(self.path / f'hrrr_{date.strftime("%Y%m%d")}{cycle:02d}.nc','w', 'NETCDF3_CLASSIC', unlimited_dims='time')
 
     def modified_latlon(self, grbfile):
 

From f5fa0c26b518632414e033fe3d78d46bb75d891d Mon Sep 17 00:00:00 2001
From: cuill <cui86ll@gmail.com>
Date: Thu, 12 Oct 2023 11:09:52 -0400
Subject: [PATCH 15/15] minor edit on class names

---
 pyschism/forcing/source_sink/sflux2source/SourceSinkIn.py | 4 ++--
 pyschism/forcing/source_sink/sflux2source/TimeHistory.py  | 2 +-
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/pyschism/forcing/source_sink/sflux2source/SourceSinkIn.py b/pyschism/forcing/source_sink/sflux2source/SourceSinkIn.py
index cf11331a..869b7100 100644
--- a/pyschism/forcing/source_sink/sflux2source/SourceSinkIn.py
+++ b/pyschism/forcing/source_sink/sflux2source/SourceSinkIn.py
@@ -24,7 +24,7 @@ def BinA(A, B):
     return [np.logical_not(np.ma.getmask(result)), np.ma.getmask(result)]
 
 
-class SourceSinkIn():
+class SourceSinkIn:
     """ class for *.prop or other similar formats"""
     def __init__(self, filename=None, number_of_groups=2, ele_groups=[]):
         self.n_group = number_of_groups  # 0: source; 1: sink
@@ -72,7 +72,7 @@ def get_ele(self, pval):
         """return ele id with specified prop value"""
 
 
-class source_sink():
+class source_sink:
     """class for handling all source/sink inputs"""
     def __init__(self, source_dir=None, start_time_str='2000-01-01 00:00:00', timedeltas=[0.0, 86400.0*365],
                  source_eles=[], sink_eles=[], vsource_data=None, vsink_data=None):
diff --git a/pyschism/forcing/source_sink/sflux2source/TimeHistory.py b/pyschism/forcing/source_sink/sflux2source/TimeHistory.py
index aa9f4f93..6fe7e729 100644
--- a/pyschism/forcing/source_sink/sflux2source/TimeHistory.py
+++ b/pyschism/forcing/source_sink/sflux2source/TimeHistory.py
@@ -23,7 +23,7 @@ def running_mean(X, N):
     return (cumsum[N:] - cumsum[:-N]) / float(N)
 
 
-class TimeHistory():
+class TimeHistory:
 
 
     """Class for manipulating *.th"""