Repository: schism-dev/pyschism (language: python)

Issue:
CVE-2007-4559 Patch
# Patching CVE-2007-4559

Hi, we are security researchers from the Advanced Research Center at [Trellix](https://www.trellix.com). We have began a campaign to patch a widespread bug named CVE-2007-4559. CVE-2007-4559 is a 15 year old bug in the Python tarfile package. By using extract() or extractall() on a tarfile object without sanitizing input, a maliciously crafted .tar file could perform a directory path traversal attack. We found at least one unsantized extractall() in your codebase and are providing a patch for you via pull request. The patch essentially checks to see if all tarfile members will be extracted safely and throws an exception otherwise. We encourage you to use this patch or your own solution to secure against CVE-2007-4559. Further technical information about the vulnerability can be found in this [blog](https://www.trellix.com/en-us/about/newsroom/stories/research/tarfile-exploiting-the-world.html).

If you have further questions you may contact us through this projects lead researcher [Kasimir Schulz](mailto:kasimir.schulz@trellix.com).



Discussion / hints:

Fix patch:
diff --git a/examples/Bctides/gen_bctides.py b/examples/Bctides/gen_bctides.py
new file mode 100644
index 00000000..39097a3c
--- /dev/null
+++ b/examples/Bctides/gen_bctides.py
@@ -0,0 +1,70 @@
+from time import time
+import os
+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
+
+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/
+    '''
+    #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 = 396
+    end_date = start_date + timedelta(days=61)
+    outdir = './'
+    hgrid = Hgrid.open("./hgrid.gr3", crs="epsg:4326")
+
+    #elevation
+    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)
+
+    #temperature
+    ite4 = itetype.Itetype4()
+    ite2 = itetype.Itetype2(0, 1)
+
+    bctides=Bctides(
+        hgrid, 
+        iettype={'1': iet5, '2': iet5}, 
+        ifltype={'1': ifl5, '2': ifl5, '3': ifl1}, 
+        isatype={'1': isa4, '2':isa4, '3': isa2}, 
+        itetype={'1': ite4, '2':ite4, '3': ite2}
+    )
+
+    bctides.write(
+        outdir, 
+        start_date=start_date, 
+        end_date=end_date, 
+        bctides=True, 
+        overwrite=True
+    )
diff --git a/examples/Gr3/Drag/GoME+0.001.reg b/examples/Gr3/Drag/GoME+0.001.reg
new file mode 100644
index 00000000..01133ab2
--- /dev/null
+++ b/examples/Gr3/Drag/GoME+0.001.reg
@@ -0,0 +1,23 @@
+Region written by ACE/gredit
+1
+20 1
+-70.711855 41.847215
+-70.588878 41.984546
+-70.567176 42.172472
+-70.473134 42.729025
+-70.205477 43.220526
+-69.474846 43.632519
+-68.563365 43.950549
+-67.196145 44.427594
+-66.819978 44.832360
+-66.805510 44.933551
+-67.268484 45.092566
+-68.353580 45.186529
+-69.735269 45.034742
+-71.804185 43.632519
+-72.079076 43.090423
+-71.898226 42.541098
+-71.138659 42.201384
+-70.899938 42.078509
+-70.827599 42.042369
+-70.762493 41.883354
diff --git a/examples/Gr3/Drag/Lake_Charles_0.reg b/examples/Gr3/Drag/Lake_Charles_0.reg
new file mode 100644
index 00000000..607ea60e
--- /dev/null
+++ b/examples/Gr3/Drag/Lake_Charles_0.reg
@@ -0,0 +1,20 @@
+Region written by ACE/gredit
+1
+17 1
+-93.248568 30.287045
+-93.189017 30.265876
+-93.200469 30.222967
+-93.231390 30.185778
+-93.270327 30.117695
+-93.296094 30.045607
+-93.316708 30.037597
+-93.324152 29.986106
+-93.323579 29.943769
+-93.317853 29.807030
+-93.323579 29.773275
+-93.363089 29.773275
+-93.363089 29.885412
+-93.348201 29.988395
+-93.341330 30.082795
+-93.350492 30.142869
+-93.321862 30.214957
diff --git a/examples/examples_gridgr3/01_gen_drag.py b/examples/Gr3/Drag/gen_drag.py
similarity index 64%
rename from examples/examples_gridgr3/01_gen_drag.py
rename to examples/Gr3/Drag/gen_drag.py
index fdb14016..de9a63e9 100644
--- a/examples/examples_gridgr3/01_gen_drag.py
+++ b/examples/Gr3/Drag/gen_drag.py
@@ -1,16 +1,18 @@
 from pyschism.mesh.hgrid import Hgrid
 from pyschism.mesh.fgrid import DragCoefficient
 
-hgrid=Hgrid.open('/sciclone/schism10/whuang07/NWM/Case1/RUN07b_ZG/hgrid.gr3', crs='epsg:4326')
+hgrid=Hgrid.open('./hgrid.gr3', crs='epsg:4326')
 depth1=-1.0
 depth2=-3.0
 bfric_river=0.0025
 bfric_land=0.025
 fgrid=DragCoefficient.linear_with_depth(hgrid, depth1, depth2, bfric_river, bfric_land)
+
+#modify values in regions
 regions=['GoME+0.001.reg', 'Lake_Charles_0.reg']
 values=[0.001, 0.0]
-flags=[1, 0]
+flags=[1, 0] # 0: reset to given value, 1: add given value to the current value
 for reg, value, flag in zip(regions, values, flags):
-    fgrid.modify_by_region(hgrid, f'Drag_Python/{reg}', 0.001,depth1, 1)
+    fgrid.modify_by_region(hgrid, f'./{reg}', value, depth1, flag)
 
 fgrid.write('drag.gr3', overwrite=True)
diff --git a/examples/Gr3/Elev_IC/gen_elev.py b/examples/Gr3/Elev_IC/gen_elev.py
new file mode 100644
index 00000000..40f96915
--- /dev/null
+++ b/examples/Gr3/Elev_IC/gen_elev.py
@@ -0,0 +1,9 @@
+from pyschism.mesh.base import Gr3
+from pyschism.mesh.gridgr3 import ElevIc
+
+if __name__ == '__main__':
+
+    hgrid = Gr3.open('./hgrid.gr3', crs='epsg:4326') 
+
+    elevic = ElevIc.default(hgrid=hgrid, offset=0.1)
+    elevic.write('elev.ic', overwrite=True)
diff --git a/examples/Gr3/Shapiro/coastal_0.2.cpp.reg b/examples/Gr3/Shapiro/coastal_0.2.cpp.reg
new file mode 100644
index 00000000..a5b996a1
--- /dev/null
+++ b/examples/Gr3/Shapiro/coastal_0.2.cpp.reg
@@ -0,0 +1,100 @@
+
+1
+97 1
+-2087974.19892804 2735461.69114659
+-2046267.92519898 2858522.40115863
+-2035841.34196965 2913215.91823516
+-2050438.86055683 3011208.43213511
+-2027500.33426503 3093248.84537155
+-1967025.86885042 3141105.86693440
+-1848162.52744231 3200357.41374636
+-1754323.30236814 3264166.71804448
+-1664654.96156731 3280119.07046258
+-1583328.03311239 3245935.65378028
+-1479062.77871217 3216310.14785139
+-1443612.61986214 3191242.29593169
+-1395650.62825981 3182126.70366696
+-1328920.94616267 3223146.84159508
+-1295556.11362641 3248214.71207795
+-1245508.97358513 3227704.65888461
+-1220485.30601345 3236820.15783981
+-1189205.77591216 3286955.88916706
+-1157926.34976052 3327976.01260565
+-1095367.30002346 3327976.01664844
+-1034893.64740146 3341649.43003119
+-1009869.99155960 3359880.53370987
+-957737.49618650 3364438.33868968
+-916031.54254958 3357601.68995285
+-870154.91374796 3316581.57039991
+-863899.05146884 3289234.86073454
+-813851.84959525 3257330.23677857
+-767975.22361845 3259609.19461440
+-724183.98731015 3225425.72426741
+-659539.67575409 3143385.59543170
+-605321.90278038 3029440.84194120
+-544848.26862581 2865360.36119174
+-536507.01806169 2774204.62507747
+-567786.52933550 2737742.30829268
+-513568.75721901 2728626.70139270
+-413474.36214879 2749136.76160622
+-327977.07926365 2790156.88185876
+-313380.06901297 2835734.69426225
+-311294.68014425 2906380.48159590
+-309209.39297320 2999815.17512957
+-323806.50493764 3100086.51825603
+-353000.72887166 3207194.62229566
+-394706.67592320 3284677.05921303
+-419730.32554551 3398621.70152670
+-388450.81439184 3503450.95949357
+-286271.13232808 3606001.14873164
+-213285.67427848 3665252.48296850
+-123617.81714054 3713109.25278223
+-73570.61978679 3724503.70588242
+-42291.10873009 3715388.09887429
+-38120.53437502 3767802.67219048
+-25608.70961326 3792870.53580273
+20267.81168508 3817938.39941500
+49462.03556304 3820217.24550671
+76570.97226339 3820217.24550672
+78656.25944087 3854400.71612716
+143300.46703233 3904536.44335163
+176665.26526503 3929604.19564317
+166238.82937764 3995692.17947167
+128703.45678998 4095963.52259992
+168324.11655511 4189398.21613237
+205860.12695341 4280554.58522318
+235053.71302052 4330689.67947690
+272589.08560818 4364873.15009736
+312209.74537334 4430961.13392587
+322636.28295726 4494770.16034198
+360171.65554500 4510722.41694602
+504057.28436660 4549463.57974983
+593725.14149262 4560858.03284954
+718843.08408122 4565415.83635253
+760549.13276806 4576810.40077240
+752207.88234603 4640619.42718910
+731354.90884917 4679360.58999345
+687563.57299579 4738611.92422672
+702160.68494185 4800141.99323051
+773060.85582705 4845719.91695078
+873155.25071970 4889018.88325953
+987846.65615984 4934596.80698683
+1067088.07797614 4996126.98733013
+1213058.99619721 5041704.91115431
+1206803.13460117 5025752.65453882
+1081685.19017228 4945991.26010490
+1073343.93974961 4900413.44769795
+1104623.45139774 4829767.66036284
+1146329.39920358 4804699.79675655
+1225570.82246726 4836604.31001547
+1315238.68271874 4909528.94360159
+1375712.42082002 4934596.80740049
+1525853.92409735 4977895.88606557
+1659313.13598129 5025752.65830226
+1782345.81095815 5110071.75004494
+1402821.25358723 5155649.66541496
+356001.08119003 4957385.82450028
+-1476977.37949447 4209908.26326187
+-2146364.75561788 3555864.56100148
+-2219352.08710189 3104642.21781792
+-2198496.24906293 2544033.50007754
diff --git a/examples/Gr3/Shapiro/coastal_0.5_1.cpp.reg b/examples/Gr3/Shapiro/coastal_0.5_1.cpp.reg
new file mode 100644
index 00000000..9afe01cb
--- /dev/null
+++ b/examples/Gr3/Shapiro/coastal_0.5_1.cpp.reg
@@ -0,0 +1,91 @@
+
+1
+88 1
+-2110399.18520885 2759473.72383486
+-2081889.68589075 2785073.56786075
+-2063382.21061246 2828875.00377477
+-2045376.06793056 2860259.28829079
+-2041631.15575467 2889859.31897550
+-2050273.45831639 2944864.84574453
+-2066099.21736074 2999260.29785134
+-2059761.59013415 3050835.65094407
+-2045204.47150868 3096399.91130284
+-2008077.08435444 3136493.72825179
+-1978150.61985263 3155831.32141113
+-1921370.89089629 3186179.04642157
+-1863396.84021269 3212362.60420253
+-1814267.94187892 3257106.74178290
+-1768008.39598740 3283788.44727055
+-1716889.84195376 3301314.99811018
+-1678640.52151209 3309575.78697443
+-1637781.41386630 3311793.30771635
+-1571558.69822787 3290215.93267294
+-1540943.09169962 3287693.05928811
+-1528630.61329063 3293718.58695410
+-1518398.13790829 3303665.52616034
+-1490368.31685577 3299219.50771532
+-1474697.10182308 3283786.10215843
+-1441844.31707609 3251616.99261998
+-1401353.14921681 3241526.64277997
+-1373217.96164349 3250890.79753969
+-1350856.85789183 3254498.93996371
+-1329013.18087201 3245960.57331146
+-1314354.47629823 3268355.32456693
+-1299741.32237654 3275148.86019249
+-1256432.42564340 3244295.82893895
+-1234975.04701564 3238344.11763361
+-1224678.33033147 3244528.12009131
+-1240677.70001922 3265343.46807355
+-1268144.46963296 3282599.09052713
+-1273006.93803143 3291224.62897666
+-1252962.16745369 3313543.91366011
+-1247674.62477571 3324124.55235404
+-1253138.00288865 3342626.38854693
+-1267743.33668832 3345798.48187961
+-1280910.61764761 3338486.47622982
+-1293391.19487411 3353238.48602371
+-1260905.06140930 3355870.47426355
+-1244134.21608518 3369999.57121295
+-1213750.08526171 3377935.87525777
+-1170399.47392774 3373085.21418655
+-1121329.95529754 3376528.37073100
+-1114086.61240031 3402970.49031818
+-1101245.90130761 3378039.99759169
+-1086464.84718525 3367543.66329230
+-1030836.07351491 3373442.46364716
+-978429.58011614 3382330.04899383
+-920653.90342542 3370067.21322368
+-864055.36757060 3345945.81026967
+-836220.23518445 3306187.48786777
+-808123.82362033 3302131.02595357
+-744705.32699062 3323022.98375028
+-705832.26634758 3334614.91226253
+-663671.58049660 3315495.13165063
+-630645.89563965 3266171.57453546
+-584275.61149532 3234802.93564763
+-572941.16499378 3188870.08573578
+-586700.64540248 3128440.66283929
+-584973.56490607 3091235.73558275
+-560978.75302966 3026675.09217801
+-527341.88613622 2953960.51839202
+-514483.51927664 2933130.78254744
+-496156.57358228 2931248.33890320
+-472456.63392595 2868545.87998493
+-446799.36500404 2863355.00193600
+-422835.81610921 2817054.59683932
+-408438.40169585 2786797.55494828
+-361811.92086997 2793852.43468602
+-328531.34296750 2811030.89660020
+-318210.00837721 2843101.23329873
+-306038.64532001 2976010.86865396
+-316986.67249013 3027394.31296513
+-352790.65305549 3133808.23023285
+-352440.33273801 3173393.35272890
+-400717.80498957 3255997.68513505
+-428030.67124798 3327950.82494407
+-446629.96550044 3411786.06908229
+-436912.32113368 3464794.02532570
+-403951.72282368 3541328.20505098
+-348470.58474651 3597329.59237861
+-390979.85246274 4257886.12245173
+-3112804.85869375 3441796.78992251
diff --git a/examples/Gr3/Shapiro/coastal_0.5_2.cpp.reg b/examples/Gr3/Shapiro/coastal_0.5_2.cpp.reg
new file mode 100644
index 00000000..94464e2e
--- /dev/null
+++ b/examples/Gr3/Shapiro/coastal_0.5_2.cpp.reg
@@ -0,0 +1,49 @@
+
+1
+46 1
+-347729.54378964 3592328.90283183
+-287049.50749085 3633932.26631487
+-238940.60829766 3671166.00903764
+-220603.01641726 3680502.99873341
+-213230.99489887 3690855.45476839
+-207335.79283072 3717710.72795215
+-184848.75778042 3743707.79989715
+-161089.81352332 3760345.03390648
+-123576.63259316 3770143.90830229
+-96616.52158290 3766577.74763832
+-82464.16923896 3768389.45051966
+-85053.83281319 3775526.16922139
+-70136.31755149 3805436.46811014
+-42335.13910624 3832741.72637059
+-5402.84839885 3851514.81537110
+37485.50162728 3858735.17772310
+47975.08821008 3851685.63446105
+58366.43239770 3849668.35407018
+65331.74771613 3862728.97462225
+102467.88175704 3893406.54560510
+143119.78886254 3916779.50296324
+163160.45588122 3921204.80916108
+169367.14931525 3968781.87274773
+154625.28775288 4002578.92647125
+125597.78017109 4071551.94429403
+114612.51529213 4104596.88269993
+123386.81656259 4128439.25114752
+146962.08637912 4163893.71614900
+176439.25891736 4223803.51916991
+202178.90905483 4264183.69741783
+249443.33387059 4359889.60984369
+293867.24060720 4401399.22129333
+314166.55820454 4456659.46172878
+319013.84200704 4496205.91718178
+325248.06109014 4509818.55177440
+362470.87894428 4516554.51944404
+451208.41200986 4537963.22321620
+531323.70942565 4568968.81645830
+676680.41556462 4599593.57973840
+725380.40754349 4634044.56243357
+735184.97725482 4660183.65325188
+715142.18969791 4687771.05611060
+668300.75587200 4755359.15678493
+158276.73532632 4703201.81839697
+-386028.18928844 4529800.38140454
+-576207.03719000 3910267.46264677
diff --git a/examples/examples_gridgr3/03_gen_shapiro2.py b/examples/Gr3/Shapiro/gen_shapiro.py
similarity index 55%
rename from examples/examples_gridgr3/03_gen_shapiro2.py
rename to examples/Gr3/Shapiro/gen_shapiro.py
index ac016311..22a4185b 100644
--- a/examples/examples_gridgr3/03_gen_shapiro2.py
+++ b/examples/Gr3/Shapiro/gen_shapiro.py
@@ -2,22 +2,21 @@
 from pyschism.mesh.gridgr3 import Shapiro
 
 if __name__ == '__main__':
-    hgrid = Hgrid.open('/sciclone/schism10/whuang07/NWM/Case1/RUN07b_ZG/hgrid.gr3', crs='EPSG:4326')
+    hgrid = Hgrid.open('./hgrid.gr3', crs='EPSG:4326')
 
     shapiro_max = 0.5
     threshold_slope = 0.5
     depths = [-99999, 20, 50]  # tweaks in shallow waters
     shapiro_vals1 = [0.2, 0.2, 0.05]  # tweaks in shallow waters
-    regdir = '/sciclone/data10/lcui01/schism/src/Utility/Pre-Processing/NWM/Shapiro_python'
     regions = [
-        f'{regdir}/coastal_0.2.cpp.reg',
-        f'{regdir}/coastal_0.5_1.cpp.reg',
-        f'{regdir}/coastal_0.5_2.cpp.reg'
+        f'./coastal_0.2.cpp.reg',
+        f'./coastal_0.5_1.cpp.reg',
+        f'./coastal_0.5_2.cpp.reg'
     ]  # tweaks in regions, the order matters
     shapiro_vals2 = [0.2, 0.5, 0.5]  # tweaks in regions, the order matters
     i_set_add_s = [0, 0, 0]
 
     shapiro = Shapiro.slope_filter(
         hgrid, shapiro_vals1, depths, shapiro_max, threshold_slope,
-        regions, shapiro_vals2, i_set_add_s, lonc=77.07, latc=24.0)
-    shapiro.write('shapiro_pyschism.gr3', overwrite=True)
+        regions, shapiro_vals2, i_set_add_s, lonc=-77.07, latc=24.0)
+    shapiro.write('shapiro.gr3', overwrite=True)
diff --git a/examples/Gr3/gen_gr3_constant.py b/examples/Gr3/gen_gr3_constant.py
new file mode 100644
index 00000000..2ce3ad57
--- /dev/null
+++ b/examples/Gr3/gen_gr3_constant.py
@@ -0,0 +1,14 @@
+from pyschism.mesh.hgrid import Gr3
+
+if __name__ == '__main__':
+    hgrid = Gr3.open('./hgrid.gr3', crs='epsg:4326')
+    grd = hgrid.copy()
+
+    gr3_names = ['albedo', 'diffmax', 'diffmin', 'watertype', 'windrot_geo2proj']
+    values = [0.1, 1.0, 1e-6, 1.0, 0.0]
+
+    for name, value in zip(gr3_names, values):
+       grd.description = name
+       grd.nodes.values[:] = value
+       grd.write(f'{name}.gr3', overwrite=True)
+
diff --git a/examples/HYCOM/download_hycom.py b/examples/HYCOM/download_hycom.py
new file mode 100644
index 00000000..b97d242f
--- /dev/null
+++ b/examples/HYCOM/download_hycom.py
@@ -0,0 +1,37 @@
+from datetime import datetime
+import logging
+from time import time
+
+from matplotlib.transforms import Bbox
+
+from pyschism.forcing.hycom.hycom2schism import DownloadHycom
+from pyschism.mesh.hgrid import Hgrid
+
+'''
+Download hycom data for Fortran scripts.
+Default is to download data for generating initial condition (use hgrid 
+    as parameter to cover the whole domain).
+Optionally, download data for generating th.nc (bnd=True) and nu.nc (nudge=True) 
+    (use bbox as parameter to cut small domain).
+'''
+logging.basicConfig(
+    format="[%(asctime)s] %(name)s %(levelname)s: %(message)s",
+    force=True,
+)
+logger = logging.getLogger('pyschism')
+logger.setLevel(logging.INFO)
+
+if __name__ == '__main__':
+    hgrid = Hgrid.open('./hgrid.gr3', crs='epsg:4326')
+    date=datetime(2018, 1, 1)
+
+    #xmin = -65
+    #xmax = -50
+    #ymin = 0
+    #ymax = 53
+    #bbox = Bbox.from_extents(xmin, ymin, xmax, ymax)
+    hycom = DownloadHycom(hgrid=hgrid)
+    t0 = time()
+    hycom.fetch_data(date, rnday=1, bnd=False, nudge=False, sub_sample=3, outdir='./')
+    print(f'It took {(time()-t0)/60} mins to download')
+    
diff --git a/examples/HYCOM/gen_bnd.py b/examples/HYCOM/gen_bnd.py
new file mode 100644
index 00000000..681eb6c1
--- /dev/null
+++ b/examples/HYCOM/gen_bnd.py
@@ -0,0 +1,33 @@
+from datetime import datetime
+import logging
+
+from pyschism.mesh.hgrid import Hgrid
+from pyschism.forcing.hycom.hycom2schism import OpenBoundaryInventory
+
+'''
+outputs:
+    elev2D.th.nc (elev=True)
+    SAL_3D.th.nc (TS=True)
+    TEM_3D.th.nc (TS=True)
+    uv3D.th.nc   (UV=True)
+'''
+
+logging.basicConfig(
+    format="[%(asctime)s] %(name)s %(levelname)s: %(message)s",
+    force=True,
+)
+logger = logging.getLogger('pyschism')
+logger.setLevel(logging.INFO)
+
+if __name__ == "__main__":
+    hgrid=Hgrid.open('./hgrid.gr3', crs='epsg:4326')
+    vgrid='./vgrid.in'
+    outdir='./'
+    start_date=datetime(2018, 8, 1)
+    rnday=20
+
+    #boundary
+    bnd=OpenBoundaryInventory(hgrid, vgrid)
+
+    #ocean_bnd_ids - segment indices, starting from zero.
+    bnd.fetch_data(outdir, start_date, rnday, elev2D=True, TS=True, UV=True, ocean_bnd_ids=[0,1,2])
diff --git a/examples/HYCOM/gen_nudge.py b/examples/HYCOM/gen_nudge.py
new file mode 100644
index 00000000..a430ea42
--- /dev/null
+++ b/examples/HYCOM/gen_nudge.py
@@ -0,0 +1,37 @@
+from datetime import datetime, timedelta
+import logging
+
+from pyschism.mesh import Hgrid, Vgrid
+from pyschism.forcing.hycom.hycom2schism import Nudge
+
+'''
+outputs:
+    SAL_nudge.gr3
+    SAL_nu.nc
+    TEM_nudge.gr3
+    TEM_nu.nc
+'''
+
+logging.basicConfig(
+    format="[%(asctime)s] %(name)s %(levelname)s: %(message)s",
+    force=True,
+)
+logger = logging.getLogger('pyschism')
+logger.setLevel(logging.INFO)
+
+if __name__ == '__main__':
+
+
+    hgrid=Hgrid.open('./hgrid.gr3', crs='epsg:4326')
+    vgrid='./vgrid.in'
+    outdir='./'
+    start_date = datetime(2018, 1, 1)
+    rnday=30
+
+    #ocean_bnd_ids - segment indices, starting from zero.
+    nudge=Nudge(hgrid=hgrid, ocean_bnd_ids=[0, 1, 9])
+
+    #rlmax - max relax distance in m or degree
+    #rnu_day - max relax strength in days 
+    #restart = True will append to the existing nc file, works when first try doesn't break.
+    nudge.fetch_data(outdir, vgrid, start_date, rnday, restart=False, rlmax=1, rnu_day=1)
diff --git a/examples/NWM/gen_sourcesink.py b/examples/NWM/gen_sourcesink.py
new file mode 100644
index 00000000..dd83a4dc
--- /dev/null
+++ b/examples/NWM/gen_sourcesink.py
@@ -0,0 +1,51 @@
+from datetime import datetime
+from time import time
+import pathlib
+import logging
+
+from pyschism.forcing.source_sink.nwm import NationalWaterModel, NWMElementPairings
+from pyschism.mesh import Hgrid
+
+logging.basicConfig(
+    format="[%(asctime)s] %(name)s %(levelname)s: %(message)s",
+    force=True,
+)
+logging.captureWarnings(True)
+
+log_level = logging.DEBUG
+logging.getLogger('pyschism').setLevel(log_level)
+
+if __name__ == '__main__':
+
+    startdate = datetime(2017, 12, 1)
+    rnday = 396
+    hgrid = Hgrid.open("./hgrid.gr3", crs="epsg:4326")
+
+    t0 = time()
+
+    #source/sink json files, if not exist, it will call NWMElementPairings to generate.
+    sources_pairings = pathlib.Path('./sources.json')
+    sinks_pairings = pathlib.Path('./sinks.json')
+    output_directory = pathlib.Path('./')
+
+    #input directory which saves downloaded nc files
+    #cache = pathlib.Path(f'./{startdate.strftime("%Y%m%d")}')
+    #cache.mkdir(exist_ok=True, parents=True)
+    cache = pathlib.Path('./NWM_v2.1')
+
+    # check if source/sink json file exists
+    if all([sources_pairings.is_file(), sinks_pairings.is_file()]) is False:
+        pairings = NWMElementPairings(hgrid)
+        sources_pairings.parent.mkdir(exist_ok=True, parents=True)
+        pairings.save_json(sources=sources_pairings, sinks=sinks_pairings)
+    else:
+        pairings = NWMElementPairings.load_json(
+            hgrid, 
+            sources_pairings, 
+            sinks_pairings)
+
+    #check nc files, if not exist will download
+    nwm=NationalWaterModel(pairings=pairings, cache=cache)
+
+    nwm.write(output_directory, hgrid, startdate, rnday, overwrite=True)
+    print(f'It took {time()-t0} seconds to generate source/sink')
diff --git a/examples/Sflux/gen_sflux_era5.py b/examples/Sflux/gen_sflux_era5.py
new file mode 100644
index 00000000..6b1c9b9f
--- /dev/null
+++ b/examples/Sflux/gen_sflux_era5.py
@@ -0,0 +1,37 @@
+from datetime import datetime
+from time import time
+import pathlib
+
+from pyschism.forcing.nws.nws2.era5 import ERA5
+
+from pyschism.mesh.hgrid import Hgrid
+
+if __name__ == '__main__':
+
+    startdate=datetime(2018, 1, 1)
+    rnday=366
+
+    t0=time()
+    hgrid=Hgrid.open('./hgrid.gr3',crs='EPSG:4326')
+    bbox = hgrid.get_bbox('EPSG:4326', output_type='bbox')
+
+    er=ERA5()
+    outdir = pathlib.Path('./')
+    interval = 1
+    er.write(
+        outdir=outdir, 
+        start_date=startdate, 
+        rnday=rnday, 
+        air=True,   #sflux_air_1
+        rad=True,   #sflux_rad_1
+        prc=True,   #sflux_prc_1
+        bbox=bbox, 
+        output_interval=interval,  #raw data is hourly
+        overwrite=True,
+        tmpdir=outdir) #default tmpdir is system's tmp (/tmp), it is possbile /tmp has no enough space for large dataset.
+
+    #write sflux_inputs.txt
+    with open("./sflux_inputs.txt", "w") as f:
+        f.write("&sflux_inputs\n/\n")
+
+    print(f'It took {(time()-t0)/60} minutes to generate {rnday} days')
diff --git a/examples/examples_sflux/gen_sflux_gfs2.py b/examples/Sflux/gen_sflux_gfs2.py
similarity index 100%
rename from examples/examples_sflux/gen_sflux_gfs2.py
rename to examples/Sflux/gen_sflux_gfs2.py
diff --git a/examples/examples_sflux/gen_sflux_hrrr3.py b/examples/Sflux/gen_sflux_hrrr3.py
similarity index 100%
rename from examples/examples_sflux/gen_sflux_hrrr3.py
rename to examples/Sflux/gen_sflux_hrrr3.py
diff --git a/examples/examples_sflux/link_sflux.py b/examples/Sflux/link_sflux.py
similarity index 100%
rename from examples/examples_sflux/link_sflux.py
rename to examples/Sflux/link_sflux.py
diff --git a/examples/combine_sources/README b/examples/combine_sources/README
new file mode 100644
index 00000000..98af4a8f
--- /dev/null
+++ b/examples/combine_sources/README
@@ -0,0 +1,40 @@
+generate_sourcesink_nc.py
+The script will combine two sources (NWM + sflux_prc) and write vsource/msource/vsink  to one netcdf file.
+
+Required inputs:
+1. path to the run directory:
+    rundir = "/sciclone/schism10/lcui01/schism20/ICOGS/ICOGS2D/GA_RUN02"
+
+2. model start time:
+    model_start_time_str = ""2018-08-01 00:00:00"
+
+The script assumes the following files are in the rundir/:
+    hgrid.gr3
+    
+    source_sink.in   # These four 
+    vsource.th
+    vsink.th
+    msource.th 
+
+    sflux/flux_prc_1.****.nc (precipitation will be interpolated to elements and converted to volume source)
+   
+Finally, the script will write one netcdf file into rundir/Orignial+sflux_source_sink/source.nc, which looks like this:
+netcdf source {
+dimensions:
+	nsources = 1283607 ;
+	nsinks = 34 ;
+	ntracers = 2 ;
+	time_msource = 2 ;
+	time_vsource = 1489 ;
+	time_vsink = 1489 ;
+	one = 1 ;
+variables:
+	int source_elem(nsources) ;
+	float vsource(time_vsource, nsources) ;
+	float msource(time_msource, ntracers, nsources) ;
+	int sink_elem(nsinks) ;
+	float vsink(time_vsink, nsinks) ;
+	float time_step_vsource(one) ;
+	float time_step_msource(one) ;
+	float time_step_vsink(one) ;
+}
diff --git a/examples/combine_sources/generate_sourcesink_nc.py b/examples/combine_sources/generate_sourcesink_nc.py
new file mode 100644
index 00000000..0cce0dea
--- /dev/null
+++ b/examples/combine_sources/generate_sourcesink_nc.py
@@ -0,0 +1,98 @@
+from time import time
+import numpy as np
+from glob import glob 
+from scipy import spatial
+import xarray as xr
+import pandas as pd
+
+
+from pyschism.mesh import Hgrid
+from pyschism.forcing.source_sink.sflux2source.SourceSinkIn import source_sink
+
+
+def nearest_neighbour(points_a, points_b):
+    tree = spatial.cKDTree(points_b)
+    return tree.query(points_a)[1]
+
+def Sflux_var_at_points(sflux_dir, sflux_set='air_1', points=None, var_name=None, nt=None):
+    '''
+    extract variable values at specified points from SCHISM-format sflux files
+    '''
+
+    sflux_files = sorted(glob(f'{sflux_dir}/sflux_{sflux_set}.*.nc'))
+    # read the first sflux file
+    sflux = xr.open_dataset(sflux_files[0])
+    if nt is None:
+        nt = len(sflux['time'].data)
+    sflux.close()
+
+    idx = nearest_neighbour(np.c_[points[:, 0], points[:, 1]], np.c_[sflux.lon.data.ravel(), sflux.lat.data.ravel()])
+    sflux_files = sorted(glob(f'{sflux_dir}/sflux_{sflux_set}.*.nc'))
+    nt_all = nt * len(sflux_files)
+    npts = len(idx)
+
+    var_at_pts = np.empty((nt_all, npts), dtype=float) * np.nan
+    time = np.empty((nt_all), dtype='datetime64[ns]')
+
+    for i, sflux_file in enumerate(sflux_files):
+        sflux = xr.open_dataset(sflux_file)
+        var = sflux[var_name].data
+        var = var.reshape(*var.shape[:-2], -1)
+        var_at_pts[i*nt:(i+1)*nt, :] = var[:nt, idx]
+        time[i*nt:(i+1)*nt] = sflux.time.data[:nt]
+        print(f'processing {sflux_file}')
+        sflux.close()
+
+    return [time, var_at_pts]
+
+
+if __name__ == "__main__":
+    rundir = '/sciclone/schism10/lcui01/schism20/ICOGS/ICOGS2D/GA_RUN02'
+    model_start_time_str = "2018-08-01 00:00:00"
+
+    # read hgrid and get element info
+    hgrid = Hgrid.open(f'{rundir}/hgrid.gr3', crs='epsg:4326')
+    lon_ctr, lat_ctr, dp_ctr = hgrid.elements.compute_centroid()
+    #get element area (crs=26918 to project the coordinate)
+    area = hgrid.elements.get_areas(26918)
+
+    ## read the original source/sink files
+    original_source_sink = source_sink(source_dir=rundir, start_time_str=model_start_time_str)
+
+    # make another set of source/sink files based on sflux
+    # read variable values from sflux
+    var_at_ele_center = Sflux_var_at_points(f'{rundir}/sflux/', sflux_set='prc_1', points=np.c_[lon_ctr, lat_ctr], var_name='prate')
+    # read time info from sflux
+    times = pd.to_datetime(var_at_ele_center[0].astype(str))
+    start_time_str = times[0].strftime("%Y-%m-%d %H:%M:%S")
+    timedeltas = (times-times[0]).total_seconds().astype(int)
+    nt = times.shape[0]
+    # calculate new sources for all elements: prate x area
+    sources_sflux = np.zeros((nt, hgrid.elements.__len__()), dtype=float)
+    for it in range(nt):  # for each time
+        sources_sflux[it, :] = 1e-3 * var_at_ele_center[1][it] * area  # prate is in kg/m^2/s, i.e., 1e-3 m3/m2/s
+
+    ## make new source sink 
+    added_ss = source_sink(
+        source_dir=None,
+        source_eles=(np.array(range(hgrid.elements.__len__()))+1).tolist(),
+        sink_eles=[1],  # dummy
+        start_time_str=start_time_str,
+        timedeltas=timedeltas,
+        vsource_data=sources_sflux
+    )
+    added_ss.update_vars()
+
+    ## add original source and sflux source
+    print('Adding two sources:')
+    t0 = time()
+    total_ss = original_source_sink + added_ss
+    if np.isnan(total_ss.vsource.data.astype(float)).any():
+        raise Exception('nan found in sources')
+    if np.isnan(total_ss.vsink.data.astype(float)).any():
+        raise Exception('nan found in sinks')
+
+    print(f'It took {(time()-t0)/60} minutes to add two sources')
+    
+    ## write
+    total_ss.nc_writer(f'{rundir}/Orignial+sflux_source_sink/') 
diff --git a/examples/examples_gridgr3/02_gen_shapiro.py b/examples/examples_gridgr3/02_gen_shapiro.py
deleted file mode 100644
index f8d41516..00000000
--- a/examples/examples_gridgr3/02_gen_shapiro.py
+++ /dev/null
@@ -1,15 +0,0 @@
-from pyschism.mesh.hgrid import Hgrid
-from pyschism.mesh.base import Gr3
-from pyschism.mesh.gridgr3 import Shapiro
-
-hgrid=Hgrid.open('/sciclone/schism10/whuang07/NWM/Case1/RUN07b_ZG/hgrid.gr3', crs='epsg:4326')
-shapiro=Shapiro.from_binary(hgrid,lonc = -77.07, latc = 24.0)
-regions=['coastal_0.2.cpp.reg', 'coastal_0.5_1.cpp.reg', 'coastal_0.5_2.cpp.reg']
-values=[0.2, 0.5, 0.5]
-flags=[0, 0, 0]
-depth1=-99999
-#hgrid.write('hgrid_after_cpp.gr3', overwrite=True)
-
-for reg, value, flag in zip(regions, values, flags):
-    shapiro.modify_by_region(hgrid, f'Shapiro_python/{reg}', value, depth1, flag)
-shapiro.write('shapiro_region.gr3', overwrite=True)
diff --git a/examples/examples_hindcast/01_nwm.py b/examples/examples_hindcast/01_nwm.py
deleted file mode 100644
index f5b25aba..00000000
--- a/examples/examples_hindcast/01_nwm.py
+++ /dev/null
@@ -1,32 +0,0 @@
-from datetime import datetime
-from time import time
-
-import logging
-
-from pyschism.forcing import source_sink
-from pyschism.mesh import Hgrid
-
-
-logging.basicConfig(
-    format="[%(asctime)s] %(name)s %(levelname)s: %(message)s",
-    force=True,
-)
-logging.captureWarnings(True)
-
-log_level = logging.DEBUG
-source_sink.nwm.logger.setLevel(log_level)
-source_sink.base.logger.setLevel(log_level)
-
-startdate=datetime(2018, 8, 24)
-rnday=40
-hgrid = Hgrid.open("./hgrid.gr3", crs="epsg:4326")
-
-t0=time()
-nwm=source_sink.nwm.NationalWaterModel()
-nwm.write('./', hgrid, startdate, rnday, overwrite=True)
-print(f'It took {time()-t0} seconds to generate source/sink')
-nwm.pairings.save_json(sources=True, sinks=True)
-nwm.pairings.sources_gdf.to_file("sources_run07b_new.shp")
-nwm.pairings.sinks_gdf.to_file("sinks_run07b_new.shp")
-nwm.pairings.intersection.to_file("intersection_run07b_new.shp")
-nwm.pairings.reaches.to_file("reaches_run07b_new.shp")
diff --git a/examples/examples_hindcast/02_era5.py b/examples/examples_hindcast/02_era5.py
deleted file mode 100644
index 8fecf832..00000000
--- a/examples/examples_hindcast/02_era5.py
+++ /dev/null
@@ -1,18 +0,0 @@
-from datetime import datetime
-from time import time
-import pathlib
-
-from pyschism.forcing.nws.nws2.era5 import ERA5
-
-from pyschism.mesh.hgrid import Hgrid
-
-startdate=datetime(2018,7,1)
-
-t0=time()
-hgrid=Hgrid.open('hgrid.gr3',crs='EPSG:4326')
-bbox = hgrid.get_bbox('EPSG:4326', output_type='bbox')
-
-er=ERA5()
-outdir = pathlib.Path('./ERA5')
-er.write(outdir=outdir, start_date=startdate, rnday=10, air=True, rad=True, prc=True, bbox=bbox, overwrite=True)
-print(f'It took {(time()-t0)/60} minutes to generate 10 days')
diff --git a/examples/examples_sflux/gen_sflux_era5.py b/examples/examples_sflux/gen_sflux_era5.py
deleted file mode 100644
index 4f127dc9..00000000
--- a/examples/examples_sflux/gen_sflux_era5.py
+++ /dev/null
@@ -1,31 +0,0 @@
-from datetime import datetime
-from time import time
-import pathlib
-import logging
-
-from pyschism.forcing.nws.nws2.era5 import ERA5
-
-from pyschism.mesh.hgrid import Hgrid
-
-logging.basicConfig(
-    format="[%(asctime)s] %(name)s %(levelname)s: %(message)s",
-    force=True,
-)
-logging.captureWarnings(True)
-
-log_level = logging.DEBUG
-logging.getLogger('pyschism').setLevel(log_level)
-
-if __name__ == '__main__':
-
-    startdate=datetime(2017, 4, 26)
-    rnday=18
-
-    t0=time()
-    hgrid=Hgrid.open('./hgrid.gr3',crs='EPSG:4326')
-    bbox = hgrid.get_bbox('EPSG:4326', output_type='bbox')
-
-    er=ERA5()
-    outdir = pathlib.Path('./')
-    er.write(outdir=outdir, start_date=startdate, rnday=rnday, air=True, rad=True, prc=True, bbox=bbox, overwrite=True)
-    print(f'It took {(time()-t0)/60} minutes to generate {rnday} days')
diff --git a/pyschism/dates.py b/pyschism/dates.py
index 7f18bc87..9f30ecea 100644
--- a/pyschism/dates.py
+++ b/pyschism/dates.py
@@ -11,14 +11,14 @@ def __init__(self):
         self.start_date = None
 
     def __set__(self, obj, start_date: datetime):
-        start_date = nearest_cycle() if start_date is None else \
-                     localize_datetime(start_date).astimezone(pytz.utc)
-
-        if obj.end_date is not None:
-            if start_date > obj.end_date:
-                raise ValueError(
-                    'start_date is greater than end_date. '
-                    'Try clearing end_date first.')
+        if start_date is not None:
+            start_date = localize_datetime(start_date).astimezone(pytz.utc)
+
+            if obj.end_date is not None:
+                if start_date > obj.end_date:
+                    raise ValueError(
+                        'start_date is greater than end_date. '
+                        'Try clearing end_date first.')
         self.start_date = start_date
 
     def __get__(self, obj, val) -> datetime:
@@ -37,9 +37,13 @@ def __set__(self, obj, run_days: Union[timedelta, float]):
 
         if isinstance(run_days, (int, float)):
             run_days = timedelta(days=float(run_days))
+        
+        elif isinstance(run_days, datetime):
+            run_days = localize_datetime(run_days) - obj.start_date
+            
         elif not isinstance(run_days, timedelta):
             raise TypeError(
-                f'Argument run_days must be {timedelta} or float.')
+                f'Argument run_days must be {timedelta} or float, not {run_days}.')
         
         if obj.start_date is not None:
             obj.end_date = obj.start_date + run_days
@@ -60,15 +64,21 @@ def __init__(self):
 
     def __set__(self, obj, val: Union[float, timedelta, datetime]):
 
+        if isinstance(val, (int, float)):
+            val = timedelta(days=float(val))
+
         if isinstance(val, datetime):
             val = localize_datetime(val).astimezone(pytz.utc)
 
-        elif not isinstance(val, timedelta):
+        elif not isinstance(val, timedelta) and obj.start_date is not None:
             val = obj.start_date + timedelta(days=float(val))
 
-        elif isinstance(val, timedelta):
+        elif isinstance(val, timedelta) and obj.start_date is not None:
             val = obj.start_date + val
 
+        else:
+            val = None
+
         self.end_date = val
 
     def __get__(self, obj, val) -> datetime:
@@ -109,7 +119,7 @@ def nearest_zulu(input_datetime=None, method='floor'):
 
 
 def nearest_cycle(input_datetime=None, period=6, method='floor'):
-
+    # Bug: method='ceil' doesn´t work correctly for cycle 18
     if input_datetime is None:
         input_datetime = localize_datetime(datetime.utcnow())
     assert method in ['floor', 'ceil']
diff --git a/pyschism/forcing/bctides/bctides.py b/pyschism/forcing/bctides/bctides.py
index b967af5e..063e8d36 100644
--- a/pyschism/forcing/bctides/bctides.py
+++ b/pyschism/forcing/bctides/bctides.py
@@ -1,13 +1,13 @@
 from datetime import datetime, timedelta
+from functools import cached_property, lru_cache
+import pathlib
 from typing import Dict, Union
-from functools import lru_cache
 import logging
-import pathlib
 
 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.elev2d import Elev2D
@@ -29,18 +29,12 @@ def __set__(self, obj, val):
         if val is not None:
             if isinstance(val, dict):
                 for bnd_id, ibctype in val.items():
-                    if not isinstance(val, (self.bctype, type(None))):
+                    if not isinstance(ibctype, (self.bctype, type(None))):
                         raise TypeError(
-                            f"Argument {val} must be of type {self.bctype} "
+                            f"Argument {ibctype} for boundary {bnd_id} must be of type {self.bctype} "
                             f" or None, not type {type(ibctype)}."
                         )
-                # TODO
-                raise NotImplementedError("Need to find the column name")
-                idxs = obj.gdf[
-                    (obj.gdf["id"] == bnd_id & obj.gdf["id"] == bnd_id)
-                ].index.values
-                for idx in idxs:
-                    obj.gdf.at[idx, val] = obj
+                    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(
@@ -68,6 +62,7 @@ class Bctides(metaclass=BctidesMeta):
 
     start_date = dates.StartDate()
     end_date = dates.EndDate()
+    run_days = dates.RunDays()
 
     def __init__(
         self,
@@ -135,10 +130,10 @@ def write(
         start_date: datetime = None,
         end_date: Union[datetime, timedelta] = None,
         bctides: Union[bool, str] = True,
-        elev2D: Union[bool, str] = True,
-        uv3D: Union[bool, str] = True,
-        tem3D: Union[bool, str] = True,
-        sal3D: 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,
@@ -146,7 +141,7 @@ def write(
         if start_date is not None:
             self.start_date = start_date
         if end_date is not None:
-            self.end_date = end_date
+            self.run_days = end_date
         # self.tidal_database.write(path, )
         output_directory = pathlib.Path(output_directory)
         output_directory.mkdir(exist_ok=overwrite, parents=True)
@@ -155,25 +150,24 @@ 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:
-
-        #            # I admit this exec is hacky.
-        #            # pros: works well, it's simple, we don't need a return value
-        #            # cons: might be confusing to read.
-        #            # This generates all the nudges and writes the nudge files.
-        #            exec(
-        #                f"from pyschism.forcing.bctides.nudge import {tracer}_Nudge;"
-        #                f"_tracer = output_directory / f'{tracer}_nudge.gr3' if {tracer.lower()}3D is True else {tracer};"
-        #                f"_tr={tracer}_Nudge(self, data_source, rlmax=data_source.rlmax, rnu_day=data_source.rnu_day);"
-        #                f'logger.info(f"Writing {tracer} nudge to file '
-        #                + r'{_tracer}");'
-        #                "_tr.write(_tracer, overwrite=overwrite)"
-        #            )
-        #            break
+        # 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
@@ -260,18 +254,13 @@ def get_forcing_string(self, boundary, global_constituents):
 
         def get_focing_digit(bctype):
             if bctype is not None:
-                # sensitive to MRO.
-                return str(
-                    getattr(
-                        bctype, f"{bctype.__class__.__bases__[0].__name__.lower()}")
-                )
+                return bctype.forcing_digit
             return "0"
 
         line = [
             f"{len(boundary.indexes)}",
-            *[digit for digit in map(get_focing_digit, bctypes)],
+            *[str(digit) for digit in map(get_focing_digit, bctypes)],
         ]
-
         f = [" ".join(line)]
         for bctype in bctypes:
             if bctype is not None:
@@ -283,15 +272,18 @@ def get_focing_digit(bctype):
         return "\n".join(f)
 
     @property
+    def rnday(self):
+        return self.run_days
+
+    @cached_property
     def gdf(self):
-        if not hasattr(self, "_gdf"):
-            self._gdf = self.hgrid.boundaries.open.copy()
-            self._gdf["iettype"] = None
-            self._gdf["ifltype"] = None
-            self._gdf["isatype"] = None
-            self._gdf["itetype"] = None
-            self._gdf["itrtype"] = None
-        return self._gdf
+        gdf = self.hgrid.boundaries.open.copy()
+        gdf["iettype"] = None
+        gdf["ifltype"] = None
+        gdf["isatype"] = None
+        gdf["itetype"] = None
+        gdf["itrtype"] = None
+        return gdf
 
     @property
     def ntip(self):
@@ -301,77 +293,115 @@ def ntip(self):
     def nbfr(self):
         return len(self.tides.get_active_forcing_constituents())
 
-    @property
-    def rnday(self):
-        return self.end_date - self.start_date
-
-    @property
+    @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,
+            )
 
-        if not hasattr(self, "_tides"):
-
-            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,
-                        )
-                
-                @lru_cache
-                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)
-
-                @lru_cache
-                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(),
-                            ]
-                        )
-                    )
-
-            self._tides = TidalConstituentCombiner(self.gdf)
+    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 self._tides
 
+if __name__ == "__main__":
+    ad_hoc_test()
diff --git a/pyschism/forcing/bctides/bctypes.py b/pyschism/forcing/bctides/bctypes.py
index e8b86a6b..c9ecc601 100644
--- a/pyschism/forcing/bctides/bctypes.py
+++ b/pyschism/forcing/bctides/bctypes.py
@@ -1,4 +1,4 @@
-from abc import abstractmethod
+from abc import abstractmethod, abstractproperty
 
 from pyschism.forcing.base import ModelForcing
 
@@ -9,5 +9,10 @@ class BoundaryForcing(ModelForcing):
     def get_boundary_string(self, boundary) -> str:
         pass
 
+    @abstractproperty
+    def forcing_digit(self):
+        """ bctype int value """
+        raise NotImplementedError
+
 
 Bctype = BoundaryForcing
diff --git a/pyschism/forcing/bctides/fes2014.py b/pyschism/forcing/bctides/fes2014.py
new file mode 100644
index 00000000..aad2ae73
--- /dev/null
+++ b/pyschism/forcing/bctides/fes2014.py
@@ -0,0 +1,155 @@
+import logging
+import os
+import pathlib
+import glob
+
+import appdirs
+from netCDF4 import Dataset
+import numpy as np
+from scipy.interpolate import griddata
+from scipy.interpolate.fitpack2 import RectBivariateSpline
+
+from pyschism.forcing.bctides.base import TidalDataProvider
+
+
+logger = logging.getLogger(__name__)
+
+FES2014_TIDES_EXTRA = 'ocean_tide_extrapolated'
+FES2014_EASTWARD_VEL = 'eastward_velocity'
+FES2014_NORTHWARD_VEL = 'northward_velocity'
+
+
+def raise_missing_file(fpath, fname):
+    raise FileNotFoundError('\n'.join([
+        f'No FES2014 file found at "{fpath}".',
+        'New users will need to register and download a copy of '
+        f'the FES2014 NetCDF file (specifically `{fname}`) '
+        'from https://www.aviso.altimetry.fr/en/data/products/auxiliary-products/global-tide-fes.html',
+        'Once you obtain netCDF files, you can follow one of the '
+        'following options: ',
+        f'1) copy or symlink the file to "{fpath}"',
+        f'2) set the environment variable `{fname}` to point'
+        ' to the file',
+    ]))
+
+
+class FES2014(TidalDataProvider):
+
+    def __init__(self, resource=None):
+        self.resource = resource
+
+    def get_elevation(self, constituent, vertices):
+        logger.info('Querying FES2014 for elevation constituent '
+                    f'{constituent}.')
+        amp = self._get_interpolation(
+            'elevation', 'amplitude', constituent, vertices)
+        phase = self._get_interpolation(
+            'elevation', 'phase', constituent, vertices)
+        return amp, phase
+
+    def get_velocity(self, constituent, vertices):
+        logger.info('Querying FES2014 for velocity constituent '
+                    f'{constituent}.')
+        uamp = self._get_interpolation(
+            'eastward_vel', 'Ua', constituent, vertices)
+        uphase = self._get_interpolation(
+            'eastward_vel', 'Ug', constituent, vertices)
+        vamp = self._get_interpolation(
+            'northward_vel', 'Va', constituent, vertices)
+        vphase = self._get_interpolation(
+            'northward_vel', 'Vg', constituent, vertices)
+        return uamp, uphase, vamp, vphase
+
+    @property
+    def constituents(self):
+        return ['M2', 'S2', 'N2', 'K2', 'K1', 'O1', 'P1',
+                'Q1', 'Mm', 'Mf', 'M4', 'MN4', 'MS4', '2N2', 'S1']
+        if not hasattr(self, '_constituents'):
+            files = glob.glob(appdirs.user_data_dir('fes2014') + '/ocean_tide_extrapolated/*.nc')
+            self._constituents = [f.split('/')[-1].split('.')[0] for f in files]
+        return self._constituents
+
+    def _get_resource(self, variable, constituent) -> Dataset:
+        resource = self._resource[variable][constituent]
+        if resource is not None:
+            return Dataset(resource)
+        datapath = pathlib.Path(appdirs.user_data_dir('fes2014'))
+        if variable == 'elevation':
+            fname = datapath / f'{FES2014_TIDES_EXTRA}/{constituent.lower()}.nc'
+        if variable == 'eastward_vel':
+            fname = datapath / f'{FES2014_EASTWARD_VEL}/{constituent.lower()}.nc'
+        if variable == 'northward_vel':
+            fname = datapath / f'{FES2014_NORTHWARD_VEL}/{constituent.lower()}.nc'
+        return Dataset(fname)
+
+    def _get_interpolation(self, phys_var, ncvar, constituent, vertices):
+        ds = self._get_resource(phys_var, constituent) 
+        lon = ds['lon'][:]
+        lat = ds['lat'][:]
+        dxs = np.unique(np.diff(lon))
+        dys = np.unique(np.diff(lat))
+        if len(dxs) !=1 or len(dys) != 1:
+           raise ValueError(f'{phys_var}: lon, lat of {constituent}.nc not uniform! ')  
+        dx = dxs[0]
+        dy = dys[0]
+
+        #get interp index
+        xi = np.asarray(
+            [x + 360. if x < 0. else x for x in vertices[:, 0]]).flatten()
+        yi = vertices[:, 1].flatten()
+        idx = np.floor((xi - lon[0]) / dx).astype('int')
+        mask = np.nonzero((lon[idx] - xi) > 0)[0]
+        idx[mask] = idx[mask] - 1
+        idy = np.floor((yi - lat[0]) / dy).astype('int')
+        mask = np.nonzero((lat[idy] - yi) > 0)[0]
+        idy[mask] = idy[mask] - 1
+        xrat = (xi - lon[idx]) / (lon[idx+1] - lon[idx])
+        yrat = (yi - lat[idy]) / (lat[idy+1] - lat[idy])
+        if np.sum((xrat > 1) | (xrat < 0) | (yrat > 1) | (yrat < 0)) != 0:
+            raise ValueError(f'xrat or yrat > 1 or < 0')
+       
+        zi = ds[ncvar][:,:]
+        # vm = 100 junk for amplitude, vm = 370 junk for phase
+        if ncvar == 'amplitude' or ncvar == 'Ua' or ncvar == 'Va':
+            vm = 100
+            zi = zi / 100
+        if ncvar == 'phase' or ncvar == 'Ug' or ncvar == 'Vg':
+            vm = 370
+            #Put phases into [0, 360)
+            mask = zi<0
+            zi[mask] = zi[mask] + 360
+        v0 = np.c_[zi[idy, idx], zi[idy, idx+1], zi[idy+1, idx], zi[idy+1, idx+1]].T
+        vmax = v0.max(axis=0)
+        vmin = v0.min(axis=0)
+        #handle phase jump
+        if ncvar == 'phase' or ncvar == 'Ug' or ncvar == 'Vg':
+            for kk in np.nonzero(abs(vmax - vmin) > 180)[0]:
+                mask = abs(v0[:, kk] - vmax[kk]) > 180
+                v0[mask, kk] = v0[mask, kk] + 360
+        v1 = v0[0] * (1 - xrat) + v0[1] * xrat
+        v2 = v0[2] * (1 - xrat) + v0[3] * xrat
+        values = v1 * (1 - yrat) + v2 * yrat
+        mask = np.nonzero((vmax > vm) * (vmin <= vm) * (vmin >= 0))[0]
+        values[mask] = vmin[mask]
+        if sum((vmax > vm) * ((vmin > vm) | (vmin < 0))) != 0:
+            raise ValueError('All junk values for {phys_var} {constituent}')
+        
+        ds.close()
+        return values
+
+    @property
+    def resource(self):
+        return self._resource
+
+    @resource.setter
+    def resource(self, resource):
+        if resource is None:
+            resource = {'elevation': {}, 'eastward_vel': {}, 'northward_vel': {}}
+            for constituent in self.constituents:
+                for key in resource.keys():
+                    resource[key].update({constituent: None})
+        else:
+            raise NotImplementedError('Check that static files exist.')
+            for file in pathlib.Path(resource).glob('*.nc'):
+                print(file)
+        self._resource = resource
diff --git a/pyschism/forcing/bctides/iettype.py b/pyschism/forcing/bctides/iettype.py
index 1b70f071..42598ff9 100644
--- a/pyschism/forcing/bctides/iettype.py
+++ b/pyschism/forcing/bctides/iettype.py
@@ -8,11 +8,16 @@
 
 
 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):
@@ -40,15 +45,34 @@ def iettype(self) -> int:
 
 
 class TidalElevation(Iettype):
-    def __init__(self, constituents="all", database="hamtide"):
+    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="EPSG:4326")[boundary.indexes, :]
+                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}")
@@ -60,7 +84,7 @@ def get_boundary_string(self, hgrid, boundary, global_constituents: List[str] =
                     for i in range(len(boundary.indexes)):
                         f.append(f"{0:.8e} {0:.8e}")
                 else:
-                    vertices = hgrid.get_xy(crs="EPSG:4326")[boundary.indexes, :]
+                    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}")
@@ -72,6 +96,11 @@ 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
@@ -139,9 +168,25 @@ def iettype(self):
         return -1
 
 
-Iettype1 = UniformTimeHistoryElevation
-Iettype2 = ConstantElevation
-Iettype3 = TidalElevation
-Iettype4 = SpatiallyVaryingTimeHistoryElevation
-Iettype5 = TidalAndSpatiallyVaryingElevationTimeHistory
-Iettype_1 = ZeroElevation
+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
index 46e99e15..cf708039 100644
--- a/pyschism/forcing/bctides/ifltype.py
+++ b/pyschism/forcing/bctides/ifltype.py
@@ -7,11 +7,16 @@
 
 
 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):
@@ -25,13 +30,15 @@ def ifltype(self):
 
 class ConstantVelocity(Ifltype):
     def __init__(self, value):
-        raise NotImplementedError(f"{self.__class__.__name__}")
         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"):
@@ -151,10 +158,29 @@ def ifltype(self):
         return -5
 
 
-Ifltype1 = UniformTimeHistoryVelocity
-Ifltype2 = ConstantVelocity
-Ifltype3 = TidalVelocity
-Ifltype4 = SpatiallyVaryingTimeHistoryVelocity
-Ifltype5 = TidalAndSpatiallyVaryingVelocityTimeHistory
-Ifltype_1 = ZeroVelocity
-Flather = ZeroVelocity
+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
index 4c6a09ba..3ddfcb6d 100644
--- a/pyschism/forcing/bctides/isatype.py
+++ b/pyschism/forcing/bctides/isatype.py
@@ -6,11 +6,16 @@
 
 
 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):
@@ -23,14 +28,25 @@ def isatype(self) -> int:
 
 
 class ConstantSalinity(Isatype):
-    def __init__(self, value):
-        raise NotImplementedError(f"{self.__class__.__name__}")
+    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):
@@ -74,7 +90,17 @@ def isatype(self):
         return 4
 
 
-Isatype1 = UniformTimeHistorySalinity
-Isatype2 = ConstantSalinity
-Isatype3 = SalinityInitialConditions
-Isatype4 = SpatiallyVaryingTimeHistorySalinity
+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
index 227c8c73..c8add4df 100644
--- a/pyschism/forcing/bctides/itetype.py
+++ b/pyschism/forcing/bctides/itetype.py
@@ -16,6 +16,10 @@ class Itetype(Bctype):
     def itetype(self) -> int:
         pass
 
+    @property
+    def forcing_digit(self):
+        return self.itetype
+
 
 class UniformTimeHistoryTemperature(Itetype):
 
@@ -30,14 +34,26 @@ def itetype(self) -> int:
 
 class ConstantTemperature(Itetype):
 
-    # def __init__(self, value, tobc: float = 1.):
-    #     self.value = value
-    #     super().__init__(tobc)
+    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):
 
@@ -86,8 +102,21 @@ def itetype(self):
         return 4
 
 
-Itetype0 = Itetype
-Itetype1 = UniformTimeHistoryTemperature
-Itetype2 = ConstantTemperature
-Itetype3 = TemperatureInitialConditions
-Itetype4 = SpatiallyVaryingTimeHistoryTemperature
+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
index 2307fa41..f8fbdaf2 100644
--- a/pyschism/forcing/bctides/itrtype.py
+++ b/pyschism/forcing/bctides/itrtype.py
@@ -7,6 +7,10 @@ class Itrtype(Bctype):
     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 ''
 
@@ -53,7 +57,17 @@ def itrtype(self):
         return 4
 
 
-Itrtype1 = UniformTimeHistoryTracer
-Itrtype2 = ConstantTracer
-Itrtype3 = TracerInitialConditions
-Itrtype4 = SpatiallyVaryingTimeHistoryTracer
+class Itrtype1(UniformTimeHistoryTracer):
+    pass
+
+
+class Itrtype2(ConstantTracer):
+    pass
+
+
+class Itrtype3(TracerInitialConditions):
+    pass
+
+
+class Itrtype4(SpatiallyVaryingTimeHistoryTracer):
+    pass
diff --git a/pyschism/forcing/bctides/tides.py b/pyschism/forcing/bctides/tides.py
index 02a6aed7..0fa07ac6 100644
--- a/pyschism/forcing/bctides/tides.py
+++ b/pyschism/forcing/bctides/tides.py
@@ -9,6 +9,7 @@
 
 # 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
@@ -22,6 +23,7 @@
 
 class TidalDatabase(Enum):
     TPXO = TPXO
+    FES2014 = FES2014 
     HAMTIDE = HAMTIDE
 
     def _missing_(self, name):
@@ -47,6 +49,19 @@ class Tides:
     major_constituents = MAJOR_CONSTITUENTS
     minor_constituents = MINOR_CONSTITUENTS
     constituents = ALL_CONSTITUENTS
+    _nodal_factors = {}
+    _earth_equilibrium_arguments = {}
+    _tidal_species_types = {}
+    _amplitude_constants = {}
+    _amplitudes = {
+            'elevation': {},
+            'velocity': {},
+        }
+
+    _phases = {
+        'elevation': {},
+        'velocity': {},
+    }
 
     def __init__(
             self,
@@ -109,14 +124,43 @@ def __call__(self, start_date: datetime, rnday: Union[int, datetime],
                 self.get_greenwich_factor(start_date, rnday, constituent))  # FACE* # noqa:E501
 
     def get_elevation(self, constituent, vertices):
+
         if constituent.lower() == 'z0':
             return np.full((vertices.shape[0],), float(self._Z0)), \
                     np.full((vertices.shape[0],), 0.)
+
+        if constituent in self._amplitudes['elevation']:
+            amps = self._amplitudes['elevation'][constituent]
+            phases = self._phases['elevation'][constituent]
+            if isinstance(amps, (float, int)):
+                amps = np.full((vertices.shape[0],), float(amps))
+            if isinstance(phases, (float, int)):
+                phases = np.full((vertices.shape[0],), float(phases))
+            return amps, phases
+        
+        # if vertices is None:
+            # raise ValueError(
+            #     'Argument vertices must not be None for tidal constituent '
+            #     f'{constituent}.')
         return self.tidal_database.get_elevation(constituent, vertices)
 
     def get_velocity(self, constituent, vertices):
         if constituent.lower() == 'z0':
             return tuple(4*[np.full((vertices.shape[0],), 0.)])
+        
+        if constituent in self._amplitudes['velocity']:
+            amps = self._amplitudes['velocity'][constituent]
+            phases = self._phases['velocity'][constituent]
+            if isinstance(amps, (float, int)):
+                amps = np.full((vertices.shape[0],), float(amps))
+            if isinstance(phases, (float, int)):
+                phases = np.full((vertices.shape[0],), float(phases))
+            return amps, phases
+
+        # if vertices is None:
+            # raise ValueError(
+            #     'Argument vertices must not be None for tidal constituent '
+            #     f'{constituent}.')
         return self.tidal_database.get_velocity(constituent, vertices)
 
     def use_all(self, potential=True, forcing=True):
@@ -148,6 +192,35 @@ def drop_constituent(self, constituent):
             raise ValueError("Argument constituent must be one of "
                              f"{self.active_constituents}")
         self._active_constituents.pop(constituent)
+        
+    def add_constituent(
+            self,
+            name: str,
+            angular_frequency,
+            nodal_factor = 1.,
+            earth_equilibrium_argument: float = 0.,
+            amplitude_constant: float = None,
+            tidal_species_type: int = None,
+            elevation_amplitude: dict = None,
+            elevation_phase: dict = 0.,
+            velocity_amplitude: dict = None,
+            velocity_phase = 0.,
+    ):
+        # self._amplitudes['elevation'].update({name: amplitude})
+        self.orbital_frequencies.update({name: angular_frequency})
+        # self.
+        self._nodal_factors.update({name: nodal_factor})
+        self._earth_equilibrium_arguments.update({name: earth_equilibrium_argument})
+        self._tidal_species_types.update({name: tidal_species_type})
+        self._amplitude_constants.update({name: amplitude_constant})  
+        self._amplitudes['elevation'].update({name: elevation_amplitude})
+        self._amplitudes['velocity'].update({name: velocity_amplitude})
+        self._phases['elevation'].update({name: elevation_phase})
+        self._phases['velocity'].update({name: velocity_phase})
+        self._active_constituents[name] = {
+                "potential": False,
+                "forcing": True,
+            }
 
     def add_Z0(self, value):
         self._active_constituents['Z0'] = {
@@ -176,6 +249,7 @@ def get_tidal_species_type(self, constituent):
             return self.tidal_species_type[constituent]
 
     def get_orbital_frequency(self, constituent):
+        # if constituent in self.
         return self.orbital_frequencies[constituent]
 
     def get_initial_conditions(self, constituent, vertices):
@@ -191,8 +265,10 @@ def set_Z0(self, Z0):
     def _manage_dates(f: Callable):
         def decorator(self, start_date, rnday, constituent):
             val = f(self, start_date, rnday, constituent)
-            del(self.start_date_utc)
-            del(self.end_date_utc)
+            if hasattr(self, 'start_date_utc'):
+                del(self.start_date_utc)
+            if hasattr(self, 'end_date_utc'):
+                del(self.end_date_utc)
             return val
         return decorator
 
@@ -200,6 +276,8 @@ def decorator(self, start_date, rnday, constituent):
     def get_nodal_factor(self, start_date: datetime,
                          rnday: Union[float, timedelta],
                          constituent: str):
+        if constituent not in ALL_CONSTITUENTS:
+            return self._nodal_factors[constituent]
         if start_date.tzinfo is not None and \
                 start_date.tzinfo.utcoffset(start_date) is not None:
             self.start_date_utc = start_date.astimezone(timezone(timedelta(0)))
@@ -300,6 +378,8 @@ def decorator(self, start_date, rnday, constituent):
     def get_greenwich_factor(self, start_date: datetime,
                              rnday: Union[float, timedelta],
                              constituent: str):
+        if constituent in self._earth_equilibrium_arguments:
+            return self._earth_equilibrium_arguments[constituent]
         if start_date.tzinfo is not None and \
                 start_date.tzinfo.utcoffset(start_date) is not None:
             self.start_date_utc = start_date.astimezone(timezone(timedelta(0)))
@@ -495,7 +575,7 @@ def active_constituents(self):
     def all_constituents(self):
         if isinstance(self.tidal_database, HAMTIDE):
             return self.major_constituents
-        elif isinstance(self.tidal_database, TPXO):
+        elif isinstance(self.tidal_database, TPXO) or isinstance(self.tidal_database, FES2014):
             return (*self.major_constituents, *self.minor_constituents)
         else:
             raise NotImplementedError(
diff --git a/pyschism/forcing/hycom/hycom2schism.py b/pyschism/forcing/hycom/hycom2schism.py
index bceec245..fbb2eb25 100644
--- a/pyschism/forcing/hycom/hycom2schism.py
+++ b/pyschism/forcing/hycom/hycom2schism.py
@@ -11,7 +11,7 @@
 
 import numpy as np
 import scipy as sp
-from numba import jit, prange
+#from numba import jit, prange
 import netCDF4 as nc
 from netCDF4 import Dataset
 from matplotlib.transforms import Bbox
@@ -23,18 +23,35 @@
 
 logger = logging.getLogger(__name__)
 
-def convert_longitude(ds):
-    lon_name = 'lon'
-    ds['_lon_adjusted'] = xr.where(
-        ds[lon_name] > 180,
-        ds[lon_name] - 360,
-        ds[lon_name])
+def convert_longitude(ds, bbox):
+#https://stackoverflow.com/questions/53345442/about-changing-longitude-array-from-0-360-to-180-to-180-with-python-xarray
+#Light_B's solution didn't generate the correct result
+#Michael's solution works, but it takes significantly longer to write nc file (~30 mins compared with 5 mins)
+#TODO: figure out why it takes much longer with the second method
+    #lon_attr = ds.coords['lon'].attrs
+    if bbox.xmin < 0:
+        logger.info(f'Convert HYCOM longitude from [0, 360) to [-180, 180):')
+        #ds.coords['lon'] = (ds.coords['lon'] + 180) % 360 - 180
+        ds['_lon_adjusted'] = xr.where(ds['lon'] > 180, ds['lon'] - 360, ds['lon'])
+    elif bbox.xmin > 0:
+        logger.info(f'Convert HYCOM longitude from [-180, 180) to [0, 360): ')
+        #ds.coords['lon'] = (ds.coords['lon'] + 360) % 360 - 180
+        ds['_lon_adjusted'] = xr.where(ds['lon'] < 0, ds['lon'] + 360, ds['lon'])
+
+    t0 = time()
     ds = (
         ds.swap_dims({lon_name: '_lon_adjusted'})
         .sel(**{'_lon_adjusted': sorted(ds._lon_adjusted)})
         .drop(lon_name)
     )
     ds = ds.rename({'_lon_adjusted': lon_name})
+    #ds = ds.sortby(ds.lon)
+    #ds.coords['lon'].attrs = lon_attr
+    logger.info(f'swap dims took {time()-t0} seconds!')
+
+    #make sure it is clipped to the bbox
+    ds = ds.sel(lon=slice(bbox.xmin - 0.5, bbox.xmax + 0.5))
+
     return ds
 
 def get_database(date, Bbox=None):
@@ -55,13 +72,13 @@ def get_database(date, Bbox=None):
     elif date >= datetime(1994, 1, 1) and date < datetime(2016, 1, 1):
         database = f'GLBv0.08/expt_53.X/data/{date.year}'
     else:
-        logger.info('No data for {date}')
+        raise ValueError(f'No data fro {date}!') 
     return database
 
-def get_idxs(date, database, bbox):
+def get_idxs(date, database, bbox, lonc=None, latc=None):
 
-    if date.strftime("%Y-%m-%d") >= datetime.now().strftime("%Y-%m-%d"):
-        date2 = datetime.now() - timedelta(days=1)
+    if date >= datetime.utcnow():
+        date2 = datetime.utcnow() - timedelta(days=1)
         baseurl = f'https://tds.hycom.org/thredds/dodsC/{database}/FMRC/runs/GLBy0.08_930_FMRC_RUN_{date2.strftime("%Y-%m-%dT12:00:00Z")}?depth[0:1:-1],lat[0:1:-1],lon[0:1:-1],time[0:1:-1]'
     else:
         baseurl=f'https://tds.hycom.org/thredds/dodsC/{database}?lat[0:1:-1],lon[0:1:-1],time[0:1:-1],depth[0:1:-1]'
@@ -73,8 +90,22 @@ def get_idxs(date, database, bbox):
     lon=ds['lon'][:]
     lat=ds['lat'][:]
     dep=ds['depth'][:]
-    lat_idxs=np.where((lat>=bbox.ymin-2.0)&(lat<=bbox.ymax+2.0))[0]
-    lon_idxs=np.where((lon>=bbox.xmin-2.0) & (lon<=bbox.xmax+2.0))[0]
+    
+    #check if hycom's lon is the same range as schism's
+    same = True
+    if not (bbox.xmin >= lon.min() and bbox.xmax <= lon.max()):
+        same = False
+        if lon.min() >= 0:
+            logger.info(f'Convert HYCOM longitude from [0, 360) to [-180, 180):')
+            idxs = lon>=180
+            lon[idxs] = lon[idxs]-360
+        elif lon.min() <= 0:
+            logger.info(f'Convert HYCOM longitude from [-180, 180) to [0, 360):')
+            idxs = lon<=0
+            lon[idxs] = lon[idxs]+360
+
+    lat_idxs=np.where((lat>=bbox.ymin-0.5)&(lat<=bbox.ymax+0.5))[0]
+    lon_idxs=np.where((lon>=bbox.xmin-0.5) & (lon<=bbox.xmax+0.5))[0]
     lon=lon[lon_idxs]
     lat=lat[lat_idxs]
     #logger.info(lon_idxs)
@@ -86,14 +117,13 @@ def get_idxs(date, database, bbox):
     lat_idx2=lat_idxs[-1].item()
     #logger.info(f'lat_idx1 is {lat_idx1}, lat_idx2 is {lat_idx2}')
     
-    for ilon in np.arange(len(lon)):
-        if lon[ilon] > 180:
-            lon[ilon] = lon[ilon]-360.
-    #lonc=(np.max(lon)+np.min(lon))/2.0
+    if lonc is None:
+        lonc = lon.mean()
     #logger.info(f'lonc is {lonc}')
-    #latc=(np.max(lat)+np.min(lat))/2.0
+    if latc is None:
+        latc = lat.mean()
     #logger.info(f'latc is {latc}')
-    x2, y2=transform_ll_to_cpp(lon, lat)
+    x2, y2=transform_ll_to_cpp(lon, lat, lonc, latc)
 
     idxs=np.where( date == times)[0]
     #check if time_idx is empty
@@ -114,13 +144,9 @@ def get_idxs(date, database, bbox):
 
     ds.close()
 
-    return time_idx, lon_idx1, lon_idx2, lat_idx1, lat_idx2, x2, y2
+    return time_idx, lon_idx1, lon_idx2, lat_idx1, lat_idx2, x2, y2, same
 
 def transform_ll_to_cpp(lon, lat, lonc=-77.07, latc=24.0):
-    #lonc=(np.max(lon)+np.min(lon))/2.0
-    #logger.info(f'lonc is {lonc}')
-    #latc=(np.max(lat)+np.min(lat))/2.0
-    #logger.info(f'latc is {latc}')
     longitude=lon/180*np.pi
     latitude=lat/180*np.pi
     radius=6378206.4
@@ -135,6 +161,12 @@ def interp_to_points_3d(dep, y2, x2, bxyz, val):
     idxs = np.where(abs(val) > 10000)
     val[idxs] = float('nan')
 
+    if not np.all(x2[:-1] <= x2[1:]):
+        logger.info('x2 is not in stricitly ascending order! Sorting x2 and val')
+        idxs = np.argsort(x2)
+        x2 = x2[idxs]
+        val = val[:, :, idxs]
+
     val_fd = sp.interpolate.RegularGridInterpolator((dep,y2,x2),np.squeeze(val),'linear', bounds_error=False, fill_value = float('nan'))
     val_int = val_fd(bxyz)
     idxs = np.isnan(val_int)
@@ -150,6 +182,12 @@ def interp_to_points_2d(y2, x2, bxy, val):
     idxs = np.where(abs(val) > 10000)
     val[idxs] = float('nan')
 
+    if not np.all(x2[:-1] <= x2[1:]):
+        logger.info('x2 is not in stricitly ascending order! Sorting x2 and val')
+        idxs = np.argsort(x2)
+        x2 = x2[idxs]
+        val = val[:, idxs]
+
     val_fd = sp.interpolate.RegularGridInterpolator((y2,x2),np.squeeze(val),'linear', bounds_error=False, fill_value = float('nan'))
     val_int = val_fd(bxy)
     idxs = np.isnan(val_int)
@@ -177,7 +215,7 @@ def __init__(self, hgrid, vgrid=None):
         self.hgrid = hgrid
         self.vgrid = Vgrid.default() if vgrid is None else vgrid
 
-    def fetch_data(self, outdir: Union[str, os.PathLike], start_date, rnday, elev2D=True, TS=True, UV=True, adjust2D=False, lats=None, msl_shifts=None): 
+    def fetch_data(self, outdir: Union[str, os.PathLike], start_date, rnday, ocean_bnd_ids = [0], elev2D=True, TS=True, UV=True, restart=False, adjust2D=False, lats=None, msl_shifts=None): 
         outdir = pathlib.Path(outdir)
 
         self.start_date = start_date
@@ -190,8 +228,10 @@ def fetch_data(self, outdir: Union[str, os.PathLike], start_date, rnday, elev2D=
         #Get open boundary 
         gdf=self.hgrid.boundaries.open.copy()
         opbd=[]
-        for boundary in gdf.itertuples():
-            opbd.extend(list(boundary.indexes))
+        #for boundary in gdf.itertuples():
+        #    opbd.extend(list(boundary.indexes))
+        for ibnd in ocean_bnd_ids:
+            opbd.extend(list(gdf.iloc[ibnd].indexes))
         blon = self.hgrid.coords[opbd,0]
         blat = self.hgrid.coords[opbd,1]
         #logger.info(f'blon min {np.min(blon)}, max {np.max(blon)}')
@@ -226,7 +266,7 @@ def fetch_data(self, outdir: Union[str, os.PathLike], start_date, rnday, elev2D=
         one=1
         #ndt=np.zeros([ntimes])
 
-        if elev2D:
+        if elev2D and restart == False:
             #timeseries_el=np.zeros([ntimes,NOP,nComp1])
             #create netcdf 
             dst_elev = Dataset(outdir / 'elev2D.th.nc', 'w', format='NETCDF4')
@@ -246,8 +286,11 @@ def fetch_data(self, outdir: Union[str, os.PathLike], start_date, rnday, elev2D=
 
             dst_elev.createVariable('time_series', 'f', ('time', 'nOpenBndNodes', 'nLevels', 'nComponents'))
             #dst_elev['time_series'][:,:,:,:] = timeseries_el
+        elif elev2D and restart:
+            dst_elev = Dataset(outdir / 'elev2D.th.nc', 'a', format='NETCDF4')
+            time_idx_restart = dst_elev['time'][:].shape[0]
 
-        if TS:
+        if TS and restart == False:
             #timeseries_s=np.zeros([ntimes,NOP,nvrt,nComp1])
             dst_salt = Dataset(outdir / 'SAL_3D.th.nc', 'w', format='NETCDF4')
             #dimensions
@@ -268,7 +311,7 @@ def fetch_data(self, outdir: Union[str, os.PathLike], start_date, rnday, elev2D=
             #temp
             #timeseries_t=np.zeros([ntimes,NOP,nvrt,nComp1])
 
-            dst_temp =  Dataset(outdir / 'TEM_3D.th.nc', 'w', format='NETCDF4')
+            dst_temp = Dataset(outdir / 'TEM_3D.th.nc', 'w', format='NETCDF4')
             #dimensions
             dst_temp.createDimension('nOpenBndNodes', NOP)
             dst_temp.createDimension('one', one)
@@ -284,8 +327,12 @@ def fetch_data(self, outdir: Union[str, os.PathLike], start_date, rnday, elev2D=
 
             dst_temp.createVariable('time_series', 'f', ('time', 'nOpenBndNodes', 'nLevels', 'nComponents'))
             #dst_temp['time_series'][:,:,:,:] = timeseries_t
+        elif TS and restart:
+            dst_salt = Dataset(outdir / 'SAL_3D.th.nc', 'a', format='NETCDF4')
+            dst_temp = Dataset(outdir / 'TEM_3D.th.nc', 'a', format='NETCDF4')
+            time_idx_restart = dst_salt['time'][:].shape[0]
 
-        if UV:
+        if UV and restart == False:
             #timeseries_uv=np.zeros([ntimes,NOP,nvrt,nComp2])
             dst_uv = Dataset(outdir / 'uv3D.th.nc', 'w', format='NETCDF4')
             #dimensions
@@ -304,25 +351,44 @@ def fetch_data(self, outdir: Union[str, os.PathLike], start_date, rnday, elev2D=
             dst_uv.createVariable('time_series', 'f', ('time', 'nOpenBndNodes', 'nLevels', 'nComponents'))
             #dst_uv['time_series'][:,:,:,:] = timeseries_uv
 
+        elif UV and restart:
+            dst_uv = Dataset(outdir / 'uv3D.th.nc', 'a', format='NETCDF4')
+            time_idx_restart = dst_uv['time'][:].shape[0]
+
         logger.info('**** Accessing GOFS data*****')
         t0=time()
-        for it, date in enumerate(self.timevector):
 
+        if restart == False: 
+            timevector = self.timevector
+            it0 = 0
+        elif restart:
+            #restart from one day earlier
+            timevector = self.timevector[time_idx_restart-1:]
+            it0 = time_idx_restart-1
+
+        for it1, date in enumerate(timevector):
+            
+            it = it0 + it1
+        
             database=get_database(date)
             logger.info(f'Fetching data for {date} from database {database}')
 
             #loop over each open boundary
             ind1 = 0
             ind2 = 0
-            for boundary in gdf.itertuples():
-
-                opbd = list(boundary.indexes)
+            #for boundary in gdf.itertuples():
+            for ibnd in ocean_bnd_ids:
+                #opbd = list(boundary.indexes)
+                opbd = list(gdf.iloc[ibnd].indexes)
                 ind1 = ind2
                 ind2 = ind1 + len(opbd)
                 #logger.info(f'ind1 = {ind1}, ind2 = {ind2}')
                 blon = self.hgrid.coords[opbd,0]
                 blat = self.hgrid.coords[opbd,1]
-                xi,yi = transform_ll_to_cpp(blon, blat)
+                blonc = blon.mean()
+                blatc = blat.mean()
+                #logger.info(f'blonc = {blon.mean()}, blatc = {blat.mean()}')
+                xi,yi = transform_ll_to_cpp(blon, blat, blonc, blatc)
                 bxy = np.c_[yi, xi]
 
                 if TS or UV:
@@ -337,21 +403,12 @@ def fetch_data(self, outdir: Union[str, os.PathLike], start_date, rnday, elev2D=
 
                 xmin, xmax = np.min(blon), np.max(blon)
                 ymin, ymax = np.min(blat), np.max(blat)
+                bbox = Bbox.from_extents(xmin, ymin, xmax, ymax)
 
-                if date.strftime("%Y-%m-%d") >= datetime(2017, 2, 1).strftime("%Y-%m-%d") and  \
-                    date.strftime("%Y-%m-%d") < datetime(2017, 6, 1).strftime("%Y-%m-%d") or \
-                    date.strftime("%Y-%m-%d") >= datetime(2017, 10, 1).strftime("%Y-%m-%d"):
-                    xmin = xmin + 360. if xmin < 0 else xmin
-                    xmax = xmax + 360. if xmax < 0 else xmax
-                    bbox = Bbox.from_extents(xmin, ymin, xmax, ymax)
-                else:
-                    bbox = Bbox.from_extents(xmin, ymin, xmax, ymax)
-                #logger.info(f'xmin is {xmin}, xmax is {xmax}')
-
-                time_idx, lon_idx1, lon_idx2, lat_idx1, lat_idx2, x2, y2 = get_idxs(date, database, bbox)
+                time_idx, lon_idx1, lon_idx2, lat_idx1, lat_idx2, x2, y2, _ = get_idxs(date, database, bbox, lonc=blonc, latc=blatc)
 
-                if date.strftime("%Y-%m-%d") >= datetime.now().strftime("%Y-%m-%d"):
-                    date2 = datetime.now() - timedelta(days=1)
+                if date >= datetime.utcnow():
+                    date2 = datetime.utcnow() - timedelta(days=1)
                     url = f'https://tds.hycom.org/thredds/dodsC/{database}/FMRC/runs/GLBy0.08_930_FMRC_RUN_' + \
                         f'{date2.strftime("%Y-%m-%dT12:00:00Z")}?depth[0:1:-1],lat[{lat_idx1}:1:{lat_idx2}],' + \
                         f'lon[{lon_idx1}:1:{lon_idx2}],time[{time_idx}],' + \
@@ -374,7 +431,7 @@ def fetch_data(self, outdir: Union[str, os.PathLike], start_date, rnday, elev2D=
                 ds=Dataset(url)
                 dep=ds['depth'][:]
 
-                logger.info('****Interpolation starts****')
+                logger.info(f'****Interpolation starts for boundary {ibnd}****')
 
                 #ndt[it]=it*24*3600.
 
@@ -433,83 +490,79 @@ def fetch_data(self, outdir: Union[str, os.PathLike], start_date, rnday, elev2D=
 
 class Nudge:
 
-    def __init__(self):
-
-        self.include = None
+    def __init__(self, hgrid=None, ocean_bnd_ids=None):
 
+        if hgrid is None:
+            raise ValueError('No hgrid information!') 
+        else:
+            self.hgrid = hgrid
 
-    def gen_nudge(self, outdir: Union[str, os.PathLike], hgrid, rlmax = 1.5, rnu_day=0.25):
+        
+        if ocean_bnd_ids is None:
+            raise ValueError('Please specify indexes for ocean boundaries!') 
+        else:
+            self.ocean_bnd_ids = ocean_bnd_ids
 
-        @jit(nopython=True, parallel=True)
-        def compute_nudge(lon, lat, nnode, opbd, out):
-          
-            rnu_max = 1.0 / rnu_day / 86400.0
-            rnu = 0
 
-            for idn in prange(nnode):
-                if idn in opbd:
-                    rnu = rnu_max
-                    distmin = 0.
-                else:
-                    distmin = np.finfo(np.float64).max
-                    for j in opbd:
-                        tmp = np.square(lon[idn]-lon[j]) + np.square(lat[idn]-lat[j])
-                        rl2 = np.sqrt(tmp)
-                        if rl2 < distmin:
-                            distmin=rl2
-                rnu = 0.
-                if distmin <= rlmax:
-                    rnu = (1-distmin/rlmax)*rnu_max
-                    #idxs_nudge[idn]=1 #idn
-                out[idn] = rnu
+    def gen_nudge(self, outdir: Union[str, os.PathLike], rlmax = 1.5, rnu_day=0.25):
 
+          
         outdir = pathlib.Path(outdir)
 
-        #get nudge zone
-        lon=hgrid.coords[:,0]
-        lat=hgrid.coords[:,1]
+        rnu_max = 1.0 / rnu_day / 86400.0
 
-        #Get open boundary 
-        gdf=hgrid.boundaries.open.copy()
-        opbd=[]
-        for boundary in gdf.itertuples():
-            opbd.extend(list(boundary.indexes))
-        opbd = np.array(opbd)
+        #get nudge zone
+        lon = self.hgrid.coords[:,0]
+        lat = self.hgrid.coords[:,1]
+        gdf = self.hgrid.boundaries.open.copy()
+        elnode = self.hgrid.elements.array
+        NE, NP = elnode.shape[0],len(lon)
+        nudge_coeff = np.zeros(NP, dtype=float)
 
-        elnode=hgrid.elements.array
-        NE, NP = len(elnode), len(lon)
+        global_idxs = {}
 
-        out = np.zeros([NP])
-        idxs_nudge=np.zeros(NP, dtype=int)
         t0 = time()
-        #compute_nudge(lon, lat, NP, opbd2, idxs_nudge, out)
-        compute_nudge(lon, lat, NP, opbd, out)
-
-        idxs=np.where(out > 0)[0]
-        idxs_nudge[idxs]=1
-        #expand nudging marker to neighbor nodes
-        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_nudge=np.delete(idxs_nudge, np.where(idxs_nudge == -99))
-        idxs=np.where(idxs_nudge == 1)[0]
-        self.include=idxs
+        for i in self.ocean_bnd_ids:
+            print(f'boundary {i}')
+            bnd_idxs = gdf.iloc[i].indexes
+
+            dis = abs((lon + 1j*lat)[:, None] - (lon[bnd_idxs] + 1j*lat[bnd_idxs])[None, :]).min(axis=1)
+            out = (1-dis/rlmax)*rnu_max
+            out[out<0] = 0
+            out[out>rnu_max] = rnu_max
+            fp = out>0
+            nudge_coeff[fp] = out[fp]
+
+            idxs_nudge=np.zeros(NP, dtype=int)
+            idxs=np.where(out > 0)[0]
+            idxs_nudge[idxs]=1
+  
+            #expand nudging marker to neighbor nodes
+            i34 = self.hgrid.elements.i34
+            fp = i34==3
+            idxs=np.where(np.max(out[elnode[fp, 0:3]], axis=1) > 0)[0]
+            idxs_nudge[elnode[fp,0:3][idxs,:]]=1
+            idxs=np.where(np.max(out[elnode[~fp, :]], axis=1) > 0)[0]
+            idxs_nudge[elnode[~fp,:][idxs,:]]=1
+
+            idxs=np.where(idxs_nudge == 1)[0]
+            global_idxs[i] = idxs
+
+
         #logger.info(f'len of nudge idxs is {len(idxs)}')
         logger.info(f'It took {time() -t0} sencods to calcuate nudge coefficient')
 
-        nudge = [f"{rlmax}, {rnu_day}"]
+        nudge = [f"rlmax={rlmax}, rnu_day={rnu_day}"]
         nudge.extend("\n")
         nudge.append(f"{NE} {NP}")
         nudge.extend("\n")
-        hgrid = hgrid.to_dict()
+        hgrid = self.hgrid.to_dict()
         nodes = hgrid['nodes']
         elements = hgrid['elements']
         for idn, (coords, values) in nodes.items():
             line = [f"{idn}"]
             line.extend([f"{x:<.7e}" for x in coords])
-            line.extend([f"{out[int(idn)-1]:<.7e}"])
+            line.extend([f"{nudge_coeff[int(idn)-1]:<.7e}"])
             line.extend("\n")
             nudge.append(" ".join(line))
 
@@ -525,9 +578,9 @@ def compute_nudge(lon, lat, nnode, opbd, out):
 
         shutil.copy2(outdir / 'TEM_nudge.gr3', outdir / 'SAL_nudge.gr3')
 
-        return self.include
+        return global_idxs
 
-    def fetch_data(self, outdir: Union[str, os.PathLike], hgrid, vgrid, start_date, rnday):
+    def fetch_data(self, outdir: Union[str, os.PathLike], vgrid, start_date, rnday, restart = False, rlmax = None, rnu_day=None):
 
         outdir = pathlib.Path(outdir)
 
@@ -541,225 +594,281 @@ def fetch_data(self, outdir: Union[str, os.PathLike], hgrid, vgrid, start_date,
         vd=Vgrid.open(vgrid)
         sigma=vd.sigma
 
+        #define nudge zone and strength 
+        rlmax = 1.5 if rlmax is None else rlmax
+        rnu_day = 0.25 if rnu_day is None else rnu_day
+        logger.info(f'Max relax distance is {rlmax} degree, max relax strengh is {rnu_day} days.')
         #Get the index for nudge
-        include = self.gen_nudge(outdir,hgrid)
+        global_idxs = self.gen_nudge(outdir, rlmax = rlmax, rnu_day=rnu_day)
 
-        #get coords of SCHISM
-        loni=hgrid.nodes.coords[:,0]
-        lati=hgrid.nodes.coords[:,1]
 
         #get bathymetry
-        depth = hgrid.values
+        depth = self.hgrid.values
 
         #compute zcor
         zcor = depth[:,None]*sigma
         nvrt=zcor.shape[1]
-        #logger.info(f'zcor at node 1098677 is {zcor[1098676,:]}')
-
-        #Get open nudge array 
-        nlon = hgrid.coords[include, 0]
-        nlat = hgrid.coords[include, 1]
-        xi,yi = transform_ll_to_cpp(nlon, nlat)
-        bxy = np.c_[yi, xi]
-
-        zcor2=zcor[include,:]
-        idxs=np.where(zcor2 > 5000)
-        #logger.info(idxs)
-        zcor2[idxs]=5000.0-1.0e-6
-        #logger.info(f'zcor2 at node 200 is {zcor2[199,:]}')
-
-        #construct schism grid
-        x2i=np.tile(xi,[nvrt,1]).T
-        y2i=np.tile(yi,[nvrt,1]).T
-        bxyz=np.c_[zcor2.reshape(np.size(zcor2)),y2i.reshape(np.size(y2i)),x2i.reshape(np.size(x2i))]
-        logger.info('Computing SCHISM zcor is done!')
+
 
         #allocate output variables
-        nNode=len(include)
-        one=1
-        ntimes=self.rnday+1
+        include = global_idxs[0]
+        nbnd = len(self.ocean_bnd_ids)
+        if nbnd > 1:
+            include = np.concatenate((global_idxs[0], global_idxs[1],))
+            if nbnd > 2:
+                for i in self.ocean_bnd_ids[2:]:
+                    include = np.concatenate((include, global_idxs[i]))
+        else:
+            include = global_idxs[0]
+
+        nNode = include.shape[0]
+        one = 1
+        ntimes = self.rnday+1
+
+        #timeseries_s=np.zeros([ntimes,nNode,nvrt,one])
+        #timeseries_t=np.zeros([ntimes,nNode,nvrt,one])
+        #ndt=np.zeros([ntimes])
+        if restart:
+            dst_temp = Dataset(outdir / 'TEM_nu.nc', 'a', format='NETCDF4')
+            dst_salt = Dataset(outdir / 'SAL_nu.nc', 'a', format='NETCDF4')
+            time_idx_restart = dst_temp['time'][:].shape[0]
+        else:
+            dst_temp = Dataset(outdir / 'TEM_nu.nc', 'w', format='NETCDF4')
+            #dimensions
+            dst_temp.createDimension('node', nNode)
+            dst_temp.createDimension('nLevels', nvrt)
+            dst_temp.createDimension('one', one)
+            dst_temp.createDimension('time', None)
+            #variables
+            dst_temp.createVariable('time', 'f', ('time',))
+            #dst_temp['time'][:] = ndt
+
+            dst_temp.createVariable('map_to_global_node', 'i4', ('node',))
+            dst_temp['map_to_global_node'][:] = include+1
+
+            dst_temp.createVariable('tracer_concentration', 'f', ('time', 'node', 'nLevels', 'one'))
+            #dst_temp['tracer_concentration'][:,:,:,:] = timeseries_t
+
+            #salinity
+            dst_salt = Dataset(outdir / 'SAL_nu.nc', 'w', format='NETCDF4')
+            #dimensions
+            dst_salt.createDimension('node', nNode)
+            dst_salt.createDimension('nLevels', nvrt)
+            dst_salt.createDimension('one', one)
+            dst_salt.createDimension('time', None)
+            #variables
+            dst_salt.createVariable('time', 'f', ('time',))
+            #dst_salt['time'][:] = ndt
+
+            dst_salt.createVariable('map_to_global_node', 'i4', ('node',))
+            dst_salt['map_to_global_node'][:] = include+1
 
-        timeseries_s=np.zeros([ntimes,nNode,nvrt,one])
-        timeseries_t=np.zeros([ntimes,nNode,nvrt,one])
-        ndt=np.zeros([ntimes])
+            dst_salt.createVariable('tracer_concentration', 'f', ('time', 'node', 'nLevels', 'one'))
+            #dst_salt['tracer_concentration'][:,:,:,:] = timeseries_s
+        
 
         logger.info('**** Accessing GOFS data*****')
+        if restart:
+            #restart from one day earlier to make sure all files consistant
+            timevector = self.timevector[time_idx_restart-1:]
+            it0 = time_idx_restart-1
+        else:
+            timevector = self.timevector
+            it0 = 0
+
         t0=time()
-        for it, date in enumerate(self.timevector):
+        for it1, date in enumerate(timevector):
+
+            it = it0 + it1
 
             database=get_database(date)
             logger.info(f'Fetching data for {date} from database {database}')
 
-            xmin, xmax = np.min(nlon), np.max(nlon)
-            ymin, ymax = np.min(nlat), np.max(nlat)
 
-            if date.strftime("%Y-%m-%d") >= datetime(2017, 2, 1).strftime("%Y-%m-%d") and  \
-                date.strftime("%Y-%m-%d") < datetime(2017, 6, 1).strftime("%Y-%m-%d") or \
-                date.strftime("%Y-%m-%d") >= datetime(2017, 10, 1).strftime("%Y-%m-%d"):
-                logger.info('Convert xmin and xmax')
-                xmin = xmin + 360. if xmin < 0 else xmin
-                xmax = xmax + 360. if xmax < 0 else xmax
-                bbox = Bbox.from_extents(xmin, ymin, xmax, ymax)
-            else:
-                bbox = Bbox.from_extents(xmin, ymin, xmax, ymax)
-            #logger.info(f'xmin is {xmin}, xmax is {xmax}')
+            ind1 = 0
+            ind2 = 0
+            for ibnd in self.ocean_bnd_ids:
+                include = global_idxs[ibnd]
+
+                ind1 = ind2
+                ind2 = ind1 + include.shape[0]
+                #Get open nudge array 
+                nlon = self.hgrid.coords[include, 0]
+                nlat = self.hgrid.coords[include, 1]
+                nlonc = nlon.mean()
+                nlatc = nlat.mean()
+                xi,yi = transform_ll_to_cpp(nlon, nlat, nlonc, nlatc)
+                bxy = np.c_[yi, xi]
 
-            time_idx, lon_idx1, lon_idx2, lat_idx1, lat_idx2, x2, y2 = get_idxs(date, database, bbox)
+                zcor2=zcor[include,:]
+                idxs=np.where(zcor2 > 5000)
+                zcor2[idxs]=5000.0-1.0e-6
 
-            if date.strftime("%Y-%m-%d") >= datetime.now().strftime("%Y-%m-%d"):
-                date2 = datetime.now() - timedelta(days=1)
-                url = f'https://tds.hycom.org/thredds/dodsC/{database}/FMRC/runs/GLBy0.08_930_FMRC_RUN_' + \
-                    f'{date2.strftime("%Y-%m-%dT12:00:00Z")}?depth[0:1:-1],lat[{lat_idx1}:1:{lat_idx2}],' + \
-                    f'lon[{lon_idx1}:1:{lon_idx2}],time[{time_idx}],' + \
-                    f'water_temp[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \
-                    f'salinity[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}]'
+                #construct schism grid
+                x2i=np.tile(xi,[nvrt,1]).T
+                y2i=np.tile(yi,[nvrt,1]).T
+                bxyz=np.c_[zcor2.reshape(np.size(zcor2)),y2i.reshape(np.size(y2i)),x2i.reshape(np.size(x2i))]
+                logger.info('Computing SCHISM zcor is done!')
 
-            else:
-                url=f'https://tds.hycom.org/thredds/dodsC/{database}?lat[{lat_idx1}:1:{lat_idx2}],' + \
-                    f'lon[{lon_idx1}:1:{lon_idx2}],depth[0:1:-1],time[{time_idx}],' + \
-                    f'water_temp[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \
-                    f'salinity[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}]'
-            #logger.info(url)
+                xmin, xmax = np.min(nlon), np.max(nlon)
+                ymin, ymax = np.min(nlat), np.max(nlat)
+                bbox = Bbox.from_extents(xmin, ymin, xmax, ymax)
+                #logger.info(f'xmin is {xmin}, xmax is {xmax}')
 
-            ds=Dataset(url)
-            salt=np.squeeze(ds['salinity'][:,:,:])
-            temp=np.squeeze(ds['water_temp'][:,:,:])
-            #logger.info(f'The shape of temp is {temp.shape}')
+                time_idx, lon_idx1, lon_idx2, lat_idx1, lat_idx2, x2, y2, _ = get_idxs(date, database, bbox, lonc=nlonc, latc=nlatc)
 
-            #Convert temp to potential temp
-            dep=ds['depth'][:]
-            ptemp = ConvertTemp(salt, temp, dep)
+                if date >= datetime.utcnow():
+                    date2 = datetime.utcnow() - timedelta(days=1)
+                    url = f'https://tds.hycom.org/thredds/dodsC/{database}/FMRC/runs/GLBy0.08_930_FMRC_RUN_' + \
+                        f'{date2.strftime("%Y-%m-%dT12:00:00Z")}?depth[0:1:-1],lat[{lat_idx1}:1:{lat_idx2}],' + \
+                        f'lon[{lon_idx1}:1:{lon_idx2}],time[{time_idx}],' + \
+                        f'water_temp[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \
+                        f'salinity[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}]'
 
-            logger.info('****Interpolation starts****')
+                else:
+                    url=f'https://tds.hycom.org/thredds/dodsC/{database}?lat[{lat_idx1}:1:{lat_idx2}],' + \
+                        f'lon[{lon_idx1}:1:{lon_idx2}],depth[0:1:-1],time[{time_idx}],' + \
+                        f'water_temp[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \
+                        f'salinity[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}]'
+                #logger.info(url)
 
-            ndt[it]=it
-            #salt
-            salt_int = interp_to_points_3d(dep, y2, x2, bxyz, salt)
-            salt_int = salt_int.reshape(zcor2.shape)
-            timeseries_s[it,:,:,0]=salt_int
+                ds=Dataset(url)
+                salt=np.squeeze(ds['salinity'][:,:,:])
+                temp=np.squeeze(ds['water_temp'][:,:,:])
+                #logger.info(f'The shape of temp is {temp.shape}')
 
-            #temp
-            temp_int = interp_to_points_3d(dep, y2, x2, bxyz, ptemp)
-            temp_int = temp_int.reshape(zcor2.shape)
-            timeseries_t[it,:,:,0]=temp_int
+                #Convert temp to potential temp
+                dep=ds['depth'][:]
+                ptemp = ConvertTemp(salt, temp, dep)
 
-            ds.close()
- 
-        with Dataset(outdir / 'TEM_nu.nc', 'w', format='NETCDF4') as dst:
-        #dimensions
-            dst.createDimension('node', nNode)
-            dst.createDimension('nLevels', nvrt)
-            dst.createDimension('one', one)
-            dst.createDimension('time', None)
-        #variables
-            dst.createVariable('time', 'f', ('time',))
-            dst['time'][:] = ndt
-
-            dst.createVariable('map_to_global_node', 'i4', ('node',))
-            dst['map_to_global_node'][:] = include+1
-
-            dst.createVariable('tracer_concentration', 'f', ('time', 'node', 'nLevels', 'one'))
-            dst['tracer_concentration'][:,:,:,:] = timeseries_t
-
-        with Dataset(outdir / 'SAL_nu.nc', 'w', format='NETCDF4') as dst:
-        #dimensions
-            dst.createDimension('node', nNode)
-            dst.createDimension('nLevels', nvrt)
-            dst.createDimension('one', one)
-            dst.createDimension('time', None)
-        #variables
-            dst.createVariable('time', 'f', ('time',))
-            dst['time'][:] = ndt
-
-            dst.createVariable('map_to_global_node', 'i4', ('node',))
-            dst['map_to_global_node'][:] = include+1
-
-            dst.createVariable('tracer_concentration', 'f', ('time', 'node', 'nLevels', 'one'))
-            dst['tracer_concentration'][:,:,:,:] = timeseries_s
+                logger.info('****Interpolation starts****')
 
+                #ndt[it]=it
+                #salt
+                dst_salt['time'][it] = it
+                salt_int = interp_to_points_3d(dep, y2, x2, bxyz, salt)
+                salt_int = salt_int.reshape(zcor2.shape)
+                #timeseries_s[it,:,:,0]=salt_int
+                dst_salt['tracer_concentration'][it,ind1:ind2,:,0] = salt_int
+
+                #temp
+                dst_temp['time'][it] = it
+                temp_int = interp_to_points_3d(dep, y2, x2, bxyz, ptemp)
+                temp_int = temp_int.reshape(zcor2.shape)
+                #timeseries_t[it,:,:,0]=temp_int
+                dst_temp['tracer_concentration'][it,ind1:ind2,:,0] = temp_int
 
+                ds.close()
+ 
+        #dst_temp.close()
+        #dst_salt.close()
+        
         logger.info(f'Writing *_nu.nc takes {time()-t0} seconds')
 
 class DownloadHycom:
 
-    def __init__(self, hgrid):
-
-        xmin, xmax = hgrid.coords[:, 0].min(), hgrid.coords[:, 0].max()
-        ymin, ymax = hgrid.coords[:, 1].min(), hgrid.coords[:, 1].max()
-        xmin = xmin + 360. if xmin < 0 else xmin
-        xmax = xmax + 360. if xmax < 0 else xmax
-        self.bbox = Bbox.from_extents(xmin, ymin, xmax, ymax)
-
-    def fetch_data(self, date):
-
-        database=get_database(date)
-        logger.info(f'Fetching data for {date} from database {database}')
-
-        time_idx, lon_idx1, lon_idx2, lat_idx1, lat_idx2, x2, y2 = get_idxs(date, database, self.bbox)
-
-        url_ssh = f'https://tds.hycom.org/thredds/dodsC/{database}?lat[{lat_idx1}:1:{lat_idx2}],' + \
-            f'lon[{lon_idx1}:1:{lon_idx2}],depth[0:1:-1],time[{time_idx}],' + \
-            f'surf_el[{time_idx}][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}]'
-        foutname = f'SSH_{date.strftime("%Y%m%d")}.nc'
-        logger.info(f'filename is {foutname}')
-        ds = xr.open_dataset(url_ssh)
-        ds = convert_longitude(ds)
-        ds = ds.rename_dims({'lon':'xlon'})
-        ds = ds.rename_dims({'lat':'ylat'})
-        ds = ds.rename_vars({'lat':'ylat'})
-        ds = ds.rename_vars({'lon':'xlon'})
-        ds.to_netcdf(foutname, 'w', 'NETCDF3_CLASSIC', unlimited_dims='time')
-        ds.close()
-
-        url_uv = f'https://tds.hycom.org/thredds/dodsC/{database}?lat[{lat_idx1}:1:{lat_idx2}],' + \
-            f'lon[{lon_idx1}:1:{lon_idx2}],depth[0:1:-1],time[{time_idx}],' + \
-            f'water_u[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \
-            f'water_v[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}]'
-
-        foutname = f'UV_{date.strftime("%Y%m%d")}.nc'
-        logger.info(f'filename is {foutname}')
-        ds = xr.open_dataset(url_uv)
-        ds = convert_longitude(ds)
-        ds = ds.rename_dims({'lon':'xlon'})
-        ds = ds.rename_dims({'lat':'ylat'})
-        ds = ds.rename_vars({'lat':'ylat'})
-        ds = ds.rename_vars({'lon':'xlon'})
-        ds.to_netcdf(foutname, 'w', 'NETCDF3_CLASSIC', unlimited_dims='time')
-        ds.close()
-
-        url_ts = f'https://tds.hycom.org/thredds/dodsC/{database}?lat[{lat_idx1}:1:{lat_idx2}],' + \
-            f'lon[{lon_idx1}:1:{lon_idx2}],depth[0:1:-1],time[{time_idx}],' + \
-            f'water_temp[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \
-            f'salinity[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}]'
-
-        foutname = f'ST_{date.strftime("%Y%m%d")}.nc'
-        logger.info(f'filename is {foutname}')
-
-        ds = xr.open_dataset(url_ts)
-        
-        temp = ds.water_temp.values
-        salt = ds.salinity.values
-        dep = ds.depth.values
-
-        #convert in-situ temperature to potential temperature
-        ptemp = ConvertTemp(salt, temp, dep)
-
-        #drop water_temp variable and add new temperature variable
-        ds = ds.drop('water_temp')
-        ds['temperature']=(['time','depth','lat','lon'], ptemp)
-        ds.temperature.attrs = {
-            'long_name': 'Sea water potential temperature',
-            'standard_name': 'sea_water_potential_temperature',
-            'units': 'degC'
-        }
-
-        #ds.assign(water_temp2=ptemp)
-        #ds.assign.attrs = ds.water_temp.attrs
-
-        ds = convert_longitude(ds)
-        ds = ds.rename_dims({'lon':'xlon'})
-        ds = ds.rename_dims({'lat':'ylat'})
-        ds = ds.rename_vars({'lat':'ylat'})
-        ds = ds.rename_vars({'lon':'xlon'})
-        ds.to_netcdf(foutname, 'w', unlimited_dims='time', encoding={'temperature':{'dtype': 'h', '_FillValue': -30000.,'scale_factor': 0.001, 'add_offset': 20., 'missing_value': -30000.}})
-        ds.close()
+    def __init__(self, hgrid=None, bbox=None):
+
+        if hgrid is None and bbox is None:
+            raise ValueError('Either hgrid or bbox must be provided!') 
+
+        if hgrid is not None:        
+            xmin, xmax = hgrid.coords[:, 0].min(), hgrid.coords[:, 0].max()
+            ymin, ymax = hgrid.coords[:, 1].min(), hgrid.coords[:, 1].max()
+            self.bbox = Bbox.from_extents(xmin, ymin, xmax, ymax)
+        elif bbox is not None:
+            self.bbox = bbox
+
+    def fetch_data(self, start_date, rnday=1, fmt='schism', bnd=False, nudge=False, sub_sample=1, outdir=None):
+        '''
+        start_date: datetime.datetime
+        rnday: integer
+        fmt: 'schism' - for Fortran code; 'hycom' - raw netCDF from HYCOM
+        bnd: file names are SSH_*.nc, TS_*.nc, UV_*.nc used in gen_hot_3Dth_from_hycom.f90
+        nudge: file name is TS_*.nc used in gen_nudge_from_hycom.f90
+        outdir: directory for output files
+        '''
+        if rnday == 1:
+            timevector = [start_date]
+        else:
+            timevector = np.arange(
+                start_date, start_date + timedelta(days=rnday+1), timedelta(days=1)
+            ).astype(datetime)
+
+        for i, date in enumerate(timevector):
+            database=get_database(date)
+            logger.info(f'Fetching data for {date} from database {database}')
+
+            time_idx, lon_idx1, lon_idx2, lat_idx1, lat_idx2, x2, y2, isLonSame = get_idxs(date, database, self.bbox)
+
+            if fmt == 'schism':
+                url_ssh = f'https://tds.hycom.org/thredds/dodsC/{database}?lat[{lat_idx1}:{sub_sample}:{lat_idx2}],' + \
+                    f'lon[{lon_idx1}:{sub_sample}:{lon_idx2}],depth[0:1:-1],time[{time_idx}],' + \
+                    f'surf_el[{time_idx}][{lat_idx1}:{sub_sample}:{lat_idx2}][{lon_idx1}:{sub_sample}:{lon_idx2}],' + \
+                    f'water_u[{time_idx}][0:1:39][{lat_idx1}:{sub_sample}:{lat_idx2}][{lon_idx1}:{sub_sample}:{lon_idx2}],' + \
+                    f'water_v[{time_idx}][0:1:39][{lat_idx1}:{sub_sample}:{lat_idx2}][{lon_idx1}:{sub_sample}:{lon_idx2}],' + \
+                    f'water_temp[{time_idx}][0:1:39][{lat_idx1}:{sub_sample}:{lat_idx2}][{lon_idx1}:{sub_sample}:{lon_idx2}],' + \
+                    f'salinity[{time_idx}][0:1:39][{lat_idx1}:{sub_sample}:{lat_idx2}][{lon_idx1}:{sub_sample}:{lon_idx2}]'
+
+                foutname = f'hycom_{date.strftime("%Y%m%d")}.nc'
+                #foutname = f'TS_{i+1}.nc'
+                logger.info(f'filename is {foutname}')
+                ds = xr.open_dataset(url_ssh)
+
+                #convert in-situ temperature to potential temperature
+                temp = ds.water_temp.values
+                salt = ds.salinity.values
+                dep = ds.depth.values
+
+                ptemp = ConvertTemp(salt, temp, dep)
+                #drop water_temp variable and add new temperature variable
+                ds = ds.drop('water_temp')
+                ds['temperature']=(['time','depth','lat','lon'], ptemp)
+                ds.temperature.attrs = {
+                    'long_name': 'Sea water potential temperature',
+                    'standard_name': 'sea_water_potential_temperature',
+                    'units': 'degC'
+                }
+
+                if not isLonSame:
+                    logger.info('Lon is not the same!')
+                    ds = convert_longitude(ds, self.bbox)
+
+                ds = ds.rename_dims({'lon':'xlon'})
+                ds = ds.rename_dims({'lat':'ylat'})
+                ds = ds.rename_vars({'lat':'ylat'})
+                ds = ds.rename_vars({'lon':'xlon'})
+
+                t0 =  time()
+                logger.info(f'Start writing nc file!')
+                ds.to_netcdf(foutname, 'w', unlimited_dims='time', encoding={'temperature':{'dtype': 'h', '_FillValue': -30000.,'scale_factor': 0.001, 'add_offset': 20., 'missing_value': -30000.}})
+                ds.close()
+                logger.info(f'It took {time()-t0} seconds to write nc file!')
+
+                #link output file to Fortran required files
+                dir_path = os.path.abspath(outdir)
+                logger.info(f'current dir is {dir_path}')
+                src = f'{dir_path}/{foutname}'
+                if bnd:
+                    names = ['SSH', 'TS', 'UV']
+                    for name in names:
+                        dst = f'{dir_path}/{name}_{i+1}.nc'
+                        os.symlink(src, dst)
+                elif nudge:
+                        dst = f'{dir_path}/TS_{i+1}.nc'
+                        os.symlink(src, dst)
+
+            elif fmt == 'hycom':
+                url=f'https://tds.hycom.org/thredds/dodsC/{database}?lat[{lat_idx1}:1:{lat_idx2}],' + \
+                    f'lon[{lon_idx1}:1:{lon_idx2}],depth[0:1:-1],time[{time_idx}],' + \
+                    f'surf_el[{time_idx}][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \
+                    f'water_temp[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \
+                    f'salinity[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \
+                    f'water_u[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \
+                    f'water_v[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}]'
+
+                foutname = f'hycom_{date.strftime("%Y%m%d")}.nc' 
+                ds = xr.open_dataset(url)
+                ds.to_netcdf(foutname, 'w')
+
+                ds.close()
diff --git a/pyschism/forcing/nws/nws2/era5.py b/pyschism/forcing/nws/nws2/era5.py
index 7bd77217..930d085c 100644
--- a/pyschism/forcing/nws/nws2/era5.py
+++ b/pyschism/forcing/nws/nws2/era5.py
@@ -23,36 +23,42 @@
 
 class ERA5DataInventory:
 
-    def __init__(self, start_date=None, rnday: Union[float, timedelta] = 4, bbox=None):
+    def __init__(self, start_date=None, rnday: Union[float, timedelta] = 4, bbox=None, tmpdir=None):
 
         self.start_date = start_date
         self.rnday = rnday
         self.end_date = self.start_date+timedelta(self.rnday + 1)
         self.client=cdsapi.Client()
         self._bbox = bbox
+        if tmpdir is not None:
+            self.tmpdir = tmpdir
 
-        r = self.client.retrieve(
-            'reanalysis-era5-single-levels',
-            {
-            'variable':[
-                '10m_u_component_of_wind','10m_v_component_of_wind','mean_sea_level_pressure',
-                '2m_dewpoint_temperature','2m_temperature','mean_total_precipitation_rate',
-                'mean_surface_downward_long_wave_radiation_flux','mean_surface_downward_short_wave_radiation_flux'
-                ],
-            'product_type':'reanalysis',
-            'date':f"{self.start_date.strftime('%Y-%m-%d')}/{self.end_date.strftime('%Y-%m-%d')}",
-            'time':[
-                '00:00','01:00','02:00','03:00','04:00','05:00',
-                '06:00','07:00','08:00','09:00','10:00','11:00',
-                '12:00','13:00','14:00','15:00','16:00','17:00',
-                '18:00','19:00','20:00','21:00','22:00','23:00'
-                ],
-            'area': [self._bbox.ymax+0.5, self._bbox.xmin-0.5, self._bbox.ymin-0.5, self._bbox.xmax+0.5], # North, West, South, East. Default: global
-            'format':'netcdf'
-            })
- 
         filename = self.tmpdir / f"era5_{self.start_date.strftime('%Y%m%d')}.nc"
-        r.download(filename)
+
+        if filename.is_file() == False:
+
+            r = self.client.retrieve(
+                'reanalysis-era5-single-levels',
+                {
+                'variable':[
+                    '10m_u_component_of_wind','10m_v_component_of_wind','mean_sea_level_pressure',
+                    '2m_dewpoint_temperature','2m_temperature','mean_total_precipitation_rate',
+                    'mean_surface_downward_long_wave_radiation_flux','mean_surface_downward_short_wave_radiation_flux'
+                    ],
+                'product_type':'reanalysis',
+                'date':f"{self.start_date.strftime('%Y-%m-%d')}/{self.end_date.strftime('%Y-%m-%d')}",
+                'time':[
+                    '00:00','01:00','02:00','03:00','04:00','05:00',
+                    '06:00','07:00','08:00','09:00','10:00','11:00',
+                    '12:00','13:00','14:00','15:00','16:00','17:00',
+                    '18:00','19:00','20:00','21:00','22:00','23:00'
+                    ],
+                'area': [self._bbox.ymax+0.5, self._bbox.xmin-0.5, self._bbox.ymin-0.5, self._bbox.xmax+0.5], # North, West, South, East. Default: global
+                'format':'netcdf'
+                }
+            )
+ 
+            r.download(filename)
         
     @property
     def tmpdir(self):
@@ -60,6 +66,10 @@ def tmpdir(self):
             self._tmpdir = tempfile.TemporaryDirectory()
         return pathlib.Path(self._tmpdir.name)
 
+    @tmpdir.setter
+    def tmpdir(self, value):
+        self._tmpdir = value
+
     @property
     def files(self):
         return sorted(list(self.tmpdir.glob('**/era5_*.nc')))
@@ -81,13 +91,7 @@ def lat(self):
         return self._lat
 
     def xy_grid(self):
-        lon = []
-        for x in self.lon:
-            if x > 180:
-                lon.append(x-360)
-            else:
-                lon.append(x)
-        return np.meshgrid(np.array(lon), self.lat)
+        return np.meshgrid(self.lon, self.lat)
 
 def put_sflux_fields(iday, date, timevector, ds, nx_grid, ny_grid, air, rad, prc, output_interval, OUTDIR):
     rt=pd.to_datetime(str(date))
@@ -274,6 +278,7 @@ def write(
         bbox: Bbox = None,
         overwrite: bool=False,
         output_interval: int = 1,
+        tmpdir = None,
     ):
         self.start_date=start_date
         self.rnday=rnday
@@ -294,6 +299,7 @@ def write(
             self.start_date,
             self.rnday,
             bbox,
+            tmpdir = tmpdir
         )
 
         logger.info('Finished downloading ERA5')
diff --git a/pyschism/forcing/source_sink/nwm.py b/pyschism/forcing/source_sink/nwm.py
index 321ef261..4afcf2cb 100644
--- a/pyschism/forcing/source_sink/nwm.py
+++ b/pyschism/forcing/source_sink/nwm.py
@@ -31,6 +31,7 @@
 
 from pyschism import dates
 from pyschism.mesh.base import Gr3
+from pyschism.utils.jsonencoder import NpEncoder
 
 from pyschism.forcing.source_sink.base import SourceSink, Sources, Sinks
 
@@ -39,7 +40,6 @@
 
 logger = logging.getLogger(__name__)
 
-
 class NWMElementPairings:
     def __init__(self, hgrid: Gr3, nwm_file=None, workers=-1):
 
@@ -198,13 +198,13 @@ def save_json(self, sources=True, sinks=True):
         if sources:
             logger.info(f"Saving {sources}")
             with open(sources, "w") as fh:
-                json.dump(self.sources, fh)
+                json.dump(self.sources, fh, cls=NpEncoder)
 
         sinks = "sinks.json" if sinks is True else sinks
         if sinks:
             with open(sinks, "w") as fh:
                 logger.info(f"Saving {sinks}")
-                json.dump(self.sinks, fh)
+                json.dump(self.sinks, fh, cls=NpEncoder)
 
     @staticmethod
     def load_json(hgrid, sources=None, sinks=None):
@@ -261,11 +261,14 @@ def hgrid(self):
     def gdf(self):
         if not hasattr(self, "_gdf"):
             gdf_coll = []
-            for reach_layer in [
-                reach_layer
-                for reach_layer in fiona.listlayers(self.nwm_file)
-                if "reaches" in reach_layer
-            ]:
+            #for reach_layer in [
+            #    reach_layer
+            #    for reach_layer in fiona.listlayers(self.nwm_file)
+            #    if "reaches" in reach_layer
+            #]:
+            reach_layers = ['nwm_reaches_conus']
+            for reach_layer in reach_layers:
+                logger.info(f'layer is {reach_layer}')
                 layer_crs = gpd.read_file(
                     self.nwm_file, rows=1, layer=reach_layer).crs
                 bbox = self.hgrid.get_bbox(crs=layer_crs)
@@ -276,7 +279,8 @@ def gdf(self):
                         layer=reach_layer,
                     )
                 )
-            self._gdf = pd.concat(gdf_coll)
+            #self._gdf = pd.concat(gdf_coll)
+            self._gdf = gpd.GeoDataFrame(gdf_coll[0])
         return self._gdf
 
     @property
diff --git a/pyschism/forcing/source_sink/sflux2source/SourceSinkIn.py b/pyschism/forcing/source_sink/sflux2source/SourceSinkIn.py
new file mode 100644
index 00000000..cf11331a
--- /dev/null
+++ b/pyschism/forcing/source_sink/sflux2source/SourceSinkIn.py
@@ -0,0 +1,267 @@
+import os
+import copy
+
+import numpy as np
+from scipy import interpolate
+import pandas as pd
+from netCDF4 import Dataset
+
+from pyschism.forcing.source_sink.sflux2source.TimeHistory import TimeHistory
+
+def BinA(A, B):
+    x = A
+    y = B
+
+    index = np.argsort(x)
+    sorted_x = x[index]
+    sorted_index = np.searchsorted(sorted_x, y)
+
+    yindex = np.take(index, sorted_index, mode="clip")
+    mask = x[yindex] != y
+
+    result = np.ma.array(yindex, mask=mask)
+
+    return [np.logical_not(np.ma.getmask(result)), np.ma.getmask(result)]
+
+
+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
+        if filename is not None:
+            """Read a usgs data file and initialize from its content """
+            self.source_file = filename
+            self.np_group = []
+            self.ip_group = []
+
+            with open(self.source_file) as fin:
+                for k in range(0, self.n_group):
+                    self.np_group.append(int(fin.readline().split()[0]))
+                    print("Points in Group " + str(k+1) + ": " + str(self.np_group[k]))
+                    self.ip_group.append(np.empty((self.np_group[k]), dtype=int))
+                    for i in range(0, self.np_group[k]):
+                        self.ip_group[k][i] = int(fin.readline())
+                    fin.readline()  # empty line between groups
+                    if self.np_group[k] > 0:
+                        print("p first: " + str(self.ip_group[k][0]))
+                        print("p last: " + str(self.ip_group[k][-1]))
+        else:
+            self.np_group = [len(x) for x in ele_groups]
+            self.ip_group = [np.array(x) for x in ele_groups]
+
+    def writer(self, filename=None):
+        if filename is None:
+            filename = self.source_file
+
+        with open(filename, 'w') as fout:
+            for k in range(0, self.n_group):
+                print("Points in Group " + str(k+1) + ": " + str(self.np_group[k]))
+                fout.write(f"{self.np_group[k]}\n")
+                for i in range(0, self.np_group[k]):
+                    fout.write(f"{self.ip_group[k][i]}\n")
+                fout.write("\n")  # empty line
+
+    def add_points(self, added_points=np.array([]), group=None):
+        self.np_group[group] += len(added_points)
+        self.ip_group[group] = np.r_[self.ip_group[group], added_points]
+
+    def get_max(self):
+        """get max prop value"""
+
+    def get_ele(self, pval):
+        """return ele id with specified prop value"""
+
+
+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):
+
+        dummy_source_sink = TimeHistory(
+            file_name=None, start_time_str='2000-01-01 00:00:00',
+            data_array=np.c_[np.array(timedeltas), np.array(timedeltas)*0.0], columns=['datetime', '1']
+        )
+
+        if source_dir is None:
+            nt = len(timedeltas)
+            nsources = len(source_eles)
+            nsinks = len(sink_eles)
+
+            if vsource_data is None:
+                vsource_data = np.zeros([nt, nsources])
+            if vsink_data is None:
+                vsink_data = np.zeros([nt, nsinks])
+
+            if nsources > 0:
+                self.vsource = TimeHistory(file_name=None, start_time_str=start_time_str,
+                                           data_array=np.c_[np.array(timedeltas), vsource_data],
+                                           columns=['datetime']+source_eles)
+                self.msource = TimeHistory(file_name=None, start_time_str=start_time_str,
+                                           data_array=np.c_[np.array(timedeltas), -9999*np.ones([nt, nsources]), np.zeros([nt, nsources])],
+                                           columns=['datetime']+source_eles+source_eles)
+            else:
+                source_eles = [1]  # dummy with 0 vsource
+                self.vsource = copy.deepcopy(dummy_source_sink)
+                self.msource = TimeHistory(file_name=None, start_time_str=start_time_str,
+                                           data_array=np.c_[np.array(timedeltas), -9999*np.ones([nt, 1]), np.zeros([nt, 1])],
+                                           columns=['datetime', '1', '1'])
+
+            if nsinks > 0:
+                self.vsink = TimeHistory(file_name=None, start_time_str=start_time_str,
+                                         data_array=np.c_[np.array(timedeltas), vsink_data],
+                                         columns=['datetime']+sink_eles)
+            else:
+                sink_eles = [1]  # dummy with 0 vsource
+                self.vsink = copy.deepcopy(dummy_source_sink)
+
+            self.source_sink_in = SourceSinkIn(filename=None, number_of_groups=2, ele_groups=[source_eles, sink_eles])
+
+        else:  # read from existing files
+            self.source_sink_in = SourceSinkIn(f"{source_dir}/source_sink.in")
+            self.source_dir = source_dir
+
+            if self.source_sink_in.np_group[0] > 0:
+                print('reading vsource\n')
+                self.vsource = TimeHistory(f"{source_dir}/vsource.th", start_time_str=start_time_str, columns=['datetime']+self.source_sink_in.ip_group[0].tolist())
+                print('reading msource\n')
+                self.msource = TimeHistory(f"{source_dir}/msource.th", start_time_str=start_time_str, columns=['datetime']+self.source_sink_in.ip_group[0].tolist()+self.source_sink_in.ip_group[0].tolist())
+                source_eles = self.source_sink_in.ip_group[0]
+            else:
+                self.vsource = copy.deepcopy(dummy_source_sink)
+                self.msource = TimeHistory(file_name=None, start_time_str=start_time_str,
+                                           data_array=np.c_[np.array(timedeltas), -9999*np.ones([nt, 1]), np.zeros([nt, 1])],
+                                           columns=['datetime', '1', '1'])
+                source_eles = [1]
+
+
+            if self.source_sink_in.np_group[1] > 0:
+                print('reading vsink\n')
+                self.vsink = TimeHistory(f"{source_dir}/vsink.th", start_time_str=start_time_str, columns=['datetime']+self.source_sink_in.ip_group[1].tolist())
+                self.vsink.df.columns = self.vsink.df.columns.map(str)
+                sink_eles = self.source_sink_in.ip_group[1]
+            else:  # dummy vsink
+                self.vsink = copy.deepcopy(dummy_source_sink)
+                sink_eles = [1]
+
+            # accomodate dummy sources/sinks
+            self.source_sink_in = SourceSinkIn(filename=None, number_of_groups=2, ele_groups=[source_eles, sink_eles])
+
+        self.vsource.df.columns = self.vsource.df.columns.map(str)
+        self.msource.df.columns = self.msource.df.columns.map(str)
+        self.vsink.df.columns = self.vsink.df.columns.map(str)
+
+        self.nsource = self.source_sink_in.np_group[0]
+        self.nsink = self.source_sink_in.np_group[1]
+        self.source_eles = source_eles
+        self.sink_eles = sink_eles
+        self.ntracers = 2
+    
+    def update_vars(self):
+        self.vsource.df_propagate()
+        self.msource.df_propagate()
+        self.vsink.df_propagate()
+        self.nsource = self.source_sink_in.np_group[0]
+        self.nsink = self.source_sink_in.np_group[1]
+        self.source_eles = self.source_sink_in.ip_group[0]
+        self.sink_eles = self.source_sink_in.ip_group[1]
+        self.ntracers = 2
+        
+        self.vsource.data = self.vsource.df.values[:, 1:]
+
+    def __add__(self, other):
+        '''
+        Add source_sink B to source_sink A,
+        retaining A's time stamps
+        '''
+        A = copy.deepcopy(self)
+        B = other
+        for i, source_sink in enumerate(['vsource', 'vsink']):  # 0: source; 1: sink
+            print(f'Adding type {i}:')
+            #if vsink only contains dummy sink, then skip interpolation (todo: remove dummy sinks)
+            if i == 1 and len(A.source_sink_in.ip_group[i]) <= 1: continue
+            [_, new_eles_inds] = BinA(A.source_sink_in.ip_group[i],
+                                      B.source_sink_in.ip_group[i])
+            [existing_eles_inds, _] = BinA(B.source_sink_in.ip_group[i],
+                                           A.source_sink_in.ip_group[i])
+            new_eles = B.source_sink_in.ip_group[i][new_eles_inds]
+            existing_eles = A.source_sink_in.ip_group[i][existing_eles_inds]
+
+            A.source_sink_in.add_points(added_points=new_eles, group=i)
+
+            B_ss = eval(f"B.{source_sink}")
+            A_ss = eval(f"A.{source_sink}")
+
+            #use A's time orgin
+            time_diff = (B_ss.df.iloc[0, 0] - A_ss.df.iloc[0, 0]).total_seconds()
+
+            f_interp = interpolate.interp1d(B_ss.time + time_diff, B_ss.df.iloc[:, 1:].to_numpy().T)
+            B_df_interp = pd.DataFrame(data=f_interp(A_ss.time).T, columns=B_ss.df.columns[1:])
+
+            # clip the values in case of small negative/positive values due to trunction error during interpolation
+ 
+            if i == 0:  # source
+                B_df_interp = B_df_interp.clip(0.0, None)
+            else:
+                B_df_interp = B_df_interp.clip(None, 0.0)
+ 
+            A_ss.df = copy.deepcopy(A_ss.df.join(B_df_interp[new_eles.astype(str)]))
+            A_ss.df[existing_eles.astype('str')] += B_df_interp[existing_eles.astype('str')]
+
+            A_ss.n_time = A_ss.df.shape[0]
+            A_ss.n_station = A_ss.df.shape[1] - 1
+
+        A.msource = TimeHistory(
+            file_name=None,
+            data_array=np.c_[np.array([0.0, 86400*365*100]), -9999*np.ones([2, len(A.vsource.df.columns)-1]), np.zeros([2, len(A.vsource.df.columns)-1])],
+            columns=['datetime']+A.vsource.df.columns[1:].tolist()+A.vsource.df.columns[1:].tolist())
+
+        A.update_vars()
+
+        return A
+
+    def writer(self, dirname=None):
+        if dirname is None:
+            raise Exception("dirname is required.")
+        os.makedirs(dirname, exist_ok=True)
+
+    def nc_writer(self, dirname=None):
+        if dirname is None:
+            raise Exception("dirname is required.")
+        os.makedirs(dirname, exist_ok=True)
+
+        msource_data = self.msource.data.reshape([self.msource.n_time, self.ntracers, self.nsource])
+        fout = Dataset(f'{dirname}/source.nc', 'w', format='NETCDF4')
+        fout.createDimension('nsources', self.nsource)
+        fout.createDimension('nsinks', self.nsink)
+        fout.createDimension('ntracers', self.ntracers)
+        fout.createDimension('time_msource', self.msource.n_time)
+        fout.createDimension('time_vsource', self.vsource.n_time)
+        fout.createDimension('time_vsink', self.vsink.n_time)
+        fout.createDimension('one', np.atleast_1d(1))
+
+        fout.createVariable('source_elem', 'i', ('nsources'))
+        fout['source_elem'][:] = self.source_sink_in.ip_group[0]
+
+        fout.createVariable('vsource', 'f', ('time_vsource', 'nsources'))
+        fout['vsource'][:,:] = self.vsource.data
+
+        fout.createVariable('msource', 'f', ('time_msource', 'ntracers', 'nsources'))
+        fout['msource'][:,:,:] = msource_data
+
+        fout.createVariable('sink_elem', 'i', ('nsinks'))
+        fout['sink_elem'][:] = self.source_sink_in.ip_group[1]
+
+        fout.createVariable('vsink', 'f', ('time_vsink', 'nsinks'))
+        fout['vsink'][:,:] = self.vsink.data
+        #C.vars.extend(['time_step_vsource','time_step_msource','time_step_vsink'])
+   
+        fout.createVariable('time_step_vsource', 'f', ('one'))
+        fout['time_step_vsource'][:] = self.vsource.delta_t
+
+        fout.createVariable('time_step_msource', 'f', ('one'))
+        fout['time_step_msource'][:] = self.msource.delta_t
+
+        fout.createVariable('time_step_vsink', 'f', ('one'))
+        fout['time_step_vsink'][:] = self.vsink.delta_t
+
+        fout.close()
diff --git a/pyschism/forcing/source_sink/sflux2source/TimeHistory.py b/pyschism/forcing/source_sink/sflux2source/TimeHistory.py
new file mode 100644
index 00000000..aa9f4f93
--- /dev/null
+++ b/pyschism/forcing/source_sink/sflux2source/TimeHistory.py
@@ -0,0 +1,238 @@
+from datetime import datetime, timedelta
+import pickle
+import copy
+from numbers import Integral, Real
+
+from pandas import read_csv, to_datetime, to_numeric
+import pandas as pd
+from matplotlib import pyplot
+import numpy
+
+def running_mean1(X, N):
+    Y = X * numpy.nan
+    if (X.ndim > 1):  # multi-dimensional
+        for i, x in enumerate(X.transpose()):
+            Y[:, i] = numpy.convolve(x, numpy.ones((N,))/N)[(N-1):]
+    else:  # 1-d array
+        Y = numpy.convolve(X, numpy.ones((N,))/N)[(N-1):]
+
+    return Y
+
+def running_mean(X, N):
+    cumsum = numpy.cumsum(numpy.insert(X, 0, 0))
+    return (cumsum[N:] - cumsum[:-N]) / float(N)
+
+
+class TimeHistory():
+
+
+    """Class for manipulating *.th"""
+
+    def __init__(self, file_name=None, start_time_str="2000-01-01 00:00:00", mask_val=None, sec_per_time_unit=1.0, data_array=None, columns=None):
+
+        if mask_val is None:
+            mask_val = -9999999
+        self.mask_val = mask_val
+
+        """initialize """
+        self.source_file = file_name
+        self.start_time_str = start_time_str
+        self.n_time = -99999999
+        self.n_station = -99999999
+        self.time = None
+        self.data = None
+        self.delta_t = -99999999.9
+        self.datetime = None
+        self.sec_per_time_unit = sec_per_time_unit
+
+        if file_name is None:
+            self.df = pd.DataFrame(data_array)
+        else:
+            self.df = read_csv(file_name, delim_whitespace=True,  index_col=False, header=None)
+        if columns is not None:
+            self.df.columns = columns
+
+        self.df.rename(columns={0: 'datetime'}, inplace=True)
+        self.df_propagate()  # populate vars from dataframe
+
+    def df_propagate(self):
+        [self.n_time, self.n_station] = list(self.df.shape)
+        self.n_station -= 1  # first col is time
+        if isinstance(self.df['datetime'][0], (Real, Integral)):
+            self.time = self.df['datetime'].to_numpy()
+            self.datetime = [datetime.strptime(self.start_time_str, '%Y-%m-%d %H:%M:%S')
+                             + timedelta(seconds=x*self.sec_per_time_unit) for x in self.time]
+            self.df['datetime'] = self.datetime
+        elif isinstance(self.df['datetime'][0], datetime):
+            time_delta = self.df['datetime'] - self.df['datetime'][0]
+            self.time = numpy.array([x.total_seconds() for x in time_delta])
+            self.datetime = self.df['datetime']
+        else:
+            raise Exception('unknown datatype for the first column, neither datetime or float')
+
+
+        # mask invalid values
+        self.data = self.df.iloc[:, 1:].to_numpy(dtype=float)
+        self.data[abs(self.data-float(self.mask_val)) < 1e-5] = numpy.nan
+        print("Number of times: " + str(self.n_time))
+        print("Number of stations: " + str(self.n_station))
+
+        self.delta_t = self.time[1] - self.time[0]
+        if self.delta_t < 50:  # probably in days
+            self.time *= 86400.0
+            self.delta_t *= 86400.0
+            print("Time unit converted from days to seconds")
+        else:
+            print("Time unit: seconds")
+
+    def get_running_mean(self, station_idx, window):
+        """ sort time series by time average"""
+        if station_idx is None:
+            valid_idx = range(0, self.n_station)
+        else:
+            valid_idx = station_idx
+        print(valid_idx)
+
+        self.data = self.df.iloc[:, 1:].to_numpy(dtype=float)
+        self.data[abs(self.data-float(self.mask_val)) < 1e-5] = numpy.nan
+        data_rm = running_mean(self.data[:, valid_idx], window)
+
+        return [data_rm]
+
+    def get_time_average(self, station_idx, start_time_str=None, end_time_str=None):
+        if station_idx == []:
+            valid_idx = range(0, self.n_station)
+        else:
+            valid_idx = station_idx
+        print(valid_idx)
+
+        if start_time_str is None:
+            idx1 = 0
+        else:
+            idx1 = self.get_snapshot(start_time_str)
+
+        if end_time_str is None:
+            idx2 = self.n_time
+        else:
+            idx2 = self.get_snapshot(end_time_str)
+
+        self.data = self.df.iloc[:, 1:].to_numpy(dtype=float)
+        self.data[abs(self.data-float(self.mask_val)) < 1e-5] = numpy.nan
+        data_mean = numpy.mean(self.data[idx1:idx2, valid_idx], axis=0)
+
+        return data_mean
+
+    def sort_time_average(self, station_idx):
+        """ sort time series by time average"""
+        if station_idx == []:
+            valid_idx = range(0, self.n_station)
+        else:
+            valid_idx = station_idx
+        print(valid_idx)
+
+        self.data = self.df.iloc[:, 1:].to_numpy(dtype=float)
+        self.data[abs(self.data-float(self.mask_val)) < 1e-5] = numpy.nan
+        data_mean = numpy.mean(self.data[:, valid_idx], axis=0)
+        id_sort = numpy.argsort(data_mean)
+        print("Preview on sorted: ")
+        print(data_mean[id_sort[-1:-11:-1]])
+
+        return [id_sort, data_mean]
+
+    def get_snapshot(self, this_time_str=None, this_time_s_1970=None):
+        """ get a snapshot closest to the input time"""
+        if this_time_str is None:
+            raise Exception("needs this_time_str to get snapshot")
+        this_time_sec = timedelta(datetime.strftime(this_time_str, "%Y-%m-%d %H:%M:%S")
+                                  - self.datetime[0]).seconds
+
+        t_idx = numpy.argmin(numpy.abs(self.time - this_time_sec))
+        return(t_idx)
+
+    def get_max_station(self):
+        """ get max among all stations """
+        self.data = self.df.iloc[:, 1:].to_numpy(dtype=float)
+        self.data[abs(self.data-float(self.mask_val)) < 1e-5] = numpy.nan
+        data_sum = numpy.sum(self.data, axis=0)
+        max_idx = numpy.argmax(data_sum)
+        print("Station with the largest time-average value: " + str(max_idx))
+
+        return max_idx
+
+    def simple_plot(self, idx, **args):
+        self.data = self.df.iloc[:, 1:].to_numpy(dtype=float)
+        self.data[abs(self.data-float(self.mask_val)) < 1e-5] = numpy.nan
+        pyplot.plot(self.datetime, self.data[:, idx], **args)
+
+    def time_unit_conversion(self, unit_conv):
+        """convert time unit"""
+        """Col 0 is time: unit_conv from sec to day is 1/86400"""
+        self.time = self.time * unit_conv
+
+    def data_unit_conversion(self, unit_conv):
+        """convert data unit"""
+        self.data = self.df.iloc[:, 1:].to_numpy(dtype=float)
+        self.data[abs(self.data-float(self.mask_val)) < 1e-5] = numpy.nan
+        self.data = self.data * unit_conv
+
+    def export_subset(self, station_idx=None, time_idx=None, i_reset_time=False, subset_filename=None):
+        import os
+        """extract a subset from the original *.th"""
+        """by station_idx and time_idx"""
+        """[] idx means no operation"""
+        """if i_reset_time == True, then the subset *.th starts from 0 time"""
+
+        if subset_filename is None:
+            subset_filename = os.path.abspath(self.source_file) + '.subset'
+
+        if station_idx is None:
+            station_idx = numpy.array(range(self.n_station))
+        else:
+            station_idx = numpy.array(station_idx)
+
+        if time_idx is None:
+            time_idx = numpy.array(range(self.n_time))
+        else:
+            time_idx = numpy.array(time_idx)
+
+        self.data = self.df.iloc[:, 1:].to_numpy(dtype=float)
+        self.data[abs(self.data-float(self.mask_val)) < 1e-5] = numpy.nan
+
+        sub_data = copy.deepcopy(self.data)
+        sub_time = copy.deepcopy(self.time)
+        if (len(station_idx) > 0):
+            sub_data = sub_data[:, station_idx]
+        if (len(time_idx) > 0):
+            sub_data = sub_data[time_idx, :]
+            sub_time = self.time[time_idx]
+
+        if i_reset_time:
+            new_start_time_str = datetime.strftime(datetime.strptime(self.start_time_str, "%Y-%m-%d %H:%M:%S") +timedelta(seconds=sub_time[0]*self.sec_per_time_unit), "%Y-%m-%d %H:%M:%S")
+            sub_time = sub_time - sub_time[0]
+        else:
+            new_start_time_str = self.start_time_str
+
+        with open(subset_filename, 'w') as fout:
+            for i, _ in enumerate(sub_time):
+                fout.write(str(sub_time[i]) + " " +
+                           ' '.join(map(str, sub_data[i, :])) +
+                           "\n")
+        
+        return TimeHistory(file_name=subset_filename, start_time_str=new_start_time_str, mask_val=self.mask_val, sec_per_time_unit=self.sec_per_time_unit,
+                           data_array=None, columns=None)
+
+    def writer(self, out_file_name):
+        self.df_propagate()
+        """ write self """
+        with open(out_file_name, 'w') as fout:
+            for i, _ in enumerate(self.time):
+                fout.write(str(self.time[i]) + " " +
+                           ' '.join(map(str, self.data[i, :])) +
+                           "\n")
+
+    @classmethod
+    def timestamp_microsecond(cls, epoch, utc_time):
+        """ convert date to microseconds """
+        time_diff = utc_time - epoch
+        assert time_diff.resolution == timedelta(microseconds=1)
+        return (time_diff.days * 86400 + time_diff.seconds) * 10**6 + time_diff.microseconds
diff --git a/pyschism/mesh/base.py b/pyschism/mesh/base.py
index cf0dc0d7..ee1f9681 100644
--- a/pyschism/mesh/base.py
+++ b/pyschism/mesh/base.py
@@ -268,6 +268,66 @@ def get_node_ball(self):
         ine = np.array([np.array(ine[i]) for i in np.arange(NP)], dtype='O') 
         return nne, ine
 
+    def compute_centroid(self):
+        elnode = self.array
+        NE = self.elements.__len__()
+        depth = self.nodes.values
+
+        x_centr, y_centr, dp_centr = np.zeros([3, NE])
+
+        if np.any(elnode.mask):
+            mask = elnode.mask[:, -1]
+            #centroid for tris 
+            x_centr[mask] = self.nodes.coords[elnode[mask, :3], 0].mean(axis=1)
+            y_centr[mask] = self.nodes.coords[elnode[mask, :3], 1].mean(axis=1)
+            dp_centr[mask] = depth[elnode[mask, :3]].mean(axis=1)
+
+            #centroid for quads
+            x_centr[~mask] = self.nodes.coords[elnode[~mask, :], 0].mean(axis=1)
+            y_centr[~mask] = self.nodes.coords[elnode[~mask, :], 1].mean(axis=1)
+            dp_centr[~mask] = depth[elnode[~mask, :]].mean(axis=1)
+        else:
+            x_centr = self.nodes.coords[elnode, 0].mean(axis=1)
+            y_centr = self.nodes.coords[elnode, 1].mean(axis=1)
+            dp_centr = depth[elnode].mean(axis=1)
+
+        return x_centr, y_centr, dp_centr
+
+    def get_areas(self, crs: Union[CRS, str] = None):
+        xy = self.nodes.get_xy(crs)
+        x = xy[:, 0]
+        y = xy[:, 1]
+
+        elnode = self.array
+        x1 = x[elnode[:, 0]]; y1 = y[elnode[:, 0]]
+        x2 = x[elnode[:, 1]]; y2 = y[elnode[:, 1]]
+        x3 = x[elnode[:, 2]]; y3 = y[elnode[:, 2]]
+        if np.any(elnode.mask):
+            x4 = x[elnode[:, 3]]; y4 = y[elnode[:, 3]]
+            mask = elnode.mask[:, -1]
+            x4[mask] = x1[mask]; y4[mask] = y1[mask]
+        else:
+            x4 = x1; y4 = y1
+        area=((x2-x1)*(y3-y1)-(x3-x1)*(y2-y1)+(x3-x1)*(y4-y1)-(x4-x1)*(y3-y1))/2
+      
+        return area
+        
+        ##self.gdf method is very slow
+        #if self.nodes.crs.is_geographic:
+        #    elements = []
+        #    for row in self.gdf.itertuples():
+        #        aeqd = CRS.from_user_input(
+        #            f"+proj=aeqd +R=6371000 +units=m "
+        #            f"+lat_0={row.geometry.centroid.y} +lon_0={row.geometry.centroid.x}"
+        #        )
+        #        current_to_aeqd = Transformer.from_crs(
+        #            self.nodes.crs, aeqd, always_xy=True
+        #        ).transform
+        #        elements.append(ops.transform(current_to_aeqd, row.geometry))
+        #    return [element.area for element in elements]
+        #else:
+        #    return [row.geometry.area for row in self.gdf.itertuples()]
+
     def get_triangulation_mask(self, element_mask):
 
         triangulation_mask = []
@@ -288,21 +348,6 @@ def get_triangulation_mask(self, element_mask):
 
         return np.array(triangulation_mask)
 
-    def get_areas(self):
-        if self.nodes.crs.is_geographic:
-            elements = []
-            for row in self.gdf.itertuples():
-                aeqd = CRS.from_user_input(
-                    f"+proj=aeqd +R=6371000 +units=m "
-                    f"+lat_0={row.geometry.centroid.y} +lon_0={row.geometry.centroid.x}"
-                )
-                current_to_aeqd = Transformer.from_crs(
-                    self.nodes.crs, aeqd, always_xy=True
-                ).transform
-                elements.append(ops.transform(current_to_aeqd, row.geometry))
-            return [element.area for element in elements]
-        else:
-            return [row.geometry.area for row in self.gdf.itertuples()]
 
     @property
     def array(self):
@@ -316,6 +361,15 @@ def array(self):
             self._array = array
         return self._array
 
+    @property
+    def i34(self):
+        if not hasattr(self, "_i34"):
+            if np.any(self.array.mask):
+                self._i34 = np.sum(~self.array.mask, axis=1)
+            else:
+                self._i34 = np.full(self.array.shape[0], 3)
+        return self._i34
+
     @property
     def triangles(self):
         if not hasattr(self, "_triangles"):
diff --git a/pyschism/mesh/gridgr3.py b/pyschism/mesh/gridgr3.py
index e0bf28fb..a3606b21 100644
--- a/pyschism/mesh/gridgr3.py
+++ b/pyschism/mesh/gridgr3.py
@@ -112,7 +112,8 @@ def transform_to_cpp(coords, lonc, latc):
 class Shapiro(Gr3Field):
 
     @classmethod
-    def slope_filter(cls, hgrid, shapiro_vals1, depths, shapiro_max, threshold_slope, regions, shapiro_vals2, flags, lonc, latc):
+    def slope_filter(cls, hgrid, shapiro_vals1, depths, shapiro_max, threshold_slope, 
+        regions=None, shapiro_vals2=None, flags=None, lonc=None, latc=None):
         """
         https://github.com/wzhengui/pylibs/blob/1ee35efaa2d52fa682113126d84846ba33318f99/Utility/schism_file.py#L265
         """
@@ -208,57 +209,14 @@ class IcField(Gr3Field):
     pass
 
 
-class ElevIc(IcField):
+class ElevIc(Gr3):
     @classmethod
-    def default(cls, hgrid, offset=-0.1):
-        obj = cls.constant(hgrid, 0.0)
-        mask = np.logical_and(hgrid.values > 0.0, obj.values < hgrid.values)
-        idxs = np.where(mask)
-        obj.values[idxs] = hgrid.values[idxs] + offset
+    def default(cls, hgrid, offset=0.1):
+        #Use Gr3.open to read hgrid.gr3 thus keep original bathymetry 
+        obj = hgrid.copy()
+        obj.values[:] = np.maximum(0, -hgrid.values - offset)
         return obj
 
-    @classmethod
-    def modify_by_region(cls, hgrid, region=None, lonc = -77.07, latc = 24.0, offset=0.1):
-        #convert hgrid from lon/lat to cpp
-        hgrid = hgrid.copy()
-        hgrid.nodes.transform_to_cpp(lonc, latc)
-        #xy = hgrid.nodes.coords
-        #x = xy[:, 0]
-        #y = xy[:, 1]
-
-        obj = cls.constant(hgrid, 0)
-
-        if region is not None:
-            lines=[line.strip().split() for line in open(region, 'r').readlines()]
-            data=np.squeeze(np.array([lines[3:]])).astype('float')
-            x=data[:,0]
-            y=data[:,1]
-            coords = list( zip(x, y))
-            poly = Polygon(coords)
-
-            #region is in cpp projection 
-            gdf1 = gpd.GeoDataFrame(
-                    {'geometry': [poly]})
-
-            points = [Point(*coord) for coord in hgrid.nodes.coords]
-            gdf2 = gpd.GeoDataFrame(
-                     {'geometry': points, 'index': list(range(len(points)))})
-            gdf_in = gpd.sjoin(gdf2, gdf1, op="within")
-            picks = [i.index for i in gdf_in.itertuples()]
-
-            #generate include.gr3
-            obj = cls.constant(hgrid, 0)
-            obj.values[picks] = 1
-            obj.description = 'include.reg, in: 1, out: 0'
-            obj.write('include.gr3', overwrite=True)
-
-        elevic = hgrid.copy()
-        idxs = np.where(obj.values == 0)
-        elevic.nodes.values[idxs] = np.maximum(0, -elevic.nodes.values[idxs])
-        idxs = np.where(obj.values == 1)
-        elevic.nodes.values[idxs] = np.maximum(0, -elevic.nodes.values[idxs]-offset)
-        return elevic 
-
 class TempIc(IcField):
     @classmethod
     def from_hycom(cls, gr3: Gr3, hycom: Hycom, date):
diff --git a/pyschism/param/core.py b/pyschism/param/core.py
index 1d97276e..1c682d5e 100644
--- a/pyschism/param/core.py
+++ b/pyschism/param/core.py
@@ -7,7 +7,6 @@
 import f90nml
 
 from pyschism.enums import Stratification
-from pyschism.param.schism_init import GitParamTemplate
 
 
 # class IbcType(Enum):
@@ -202,7 +201,7 @@ def template(self):
     @template.setter
     def template(self, template: Union[str, os.PathLike, None]):
         if template is True:
-            template = GitParamTemplate().path
+            template = pathlib.Path(__file__) / 'param.nml'
         elif template is False or template is None:
             template = None
         else:
@@ -212,5 +211,5 @@ def template(self, template: Union[str, os.PathLike, None]):
     @property
     def defaults(self):
         if self.template is None:
-            return GitParamTemplate().core
+            return f90nml.read(pathlib.Path(__file__).parent / 'param.nml')["core"]
         return f90nml.read(self.template)['core']
diff --git a/pyschism/param/opt.py b/pyschism/param/opt.py
index cb2e4d6e..b6bcf365 100644
--- a/pyschism/param/opt.py
+++ b/pyschism/param/opt.py
@@ -1,7 +1,7 @@
 from datetime import datetime, timedelta
 import logging
 
-# import pathlib
+import pathlib
 from typing import Union
 
 import f90nml
@@ -10,9 +10,7 @@
 
 from pyschism import dates
 from pyschism.mesh.fgrid import NchiType
-from pyschism.param.schism_init import GitParamTemplate
 
-PARAM_DEFAULTS = f90nml.read(GitParamTemplate().path)["opt"]
 
 logger = logging.getLogger(__name__)
 
@@ -42,6 +40,8 @@
 
 class OptMeta(type):
     def __new__(meta, name, bases, attrs):
+
+        PARAM_DEFAULTS = f90nml.read(pathlib.Path(__file__).parent / 'param.nml')["opt"]
         for key, value in PARAM_DEFAULTS.items():
             if key not in attrs:
                 if isinstance(value, list):
diff --git a/pyschism/param/param.nml b/pyschism/param/param.nml
new file mode 100644
index 00000000..b9453459
--- /dev/null
+++ b/pyschism/param/param.nml
@@ -0,0 +1,1104 @@
+!parameter inputs via namelist convention.
+!(1) Use ' ' (single quotes) for chars;
+!(2) integer values are fine for real vars/arrays;
+!(3) if multiple entries for a parameter are found, the last one wins - please avoid this
+!(4) array inputs follow column major (like FORTRAN) and can spill to multiple lines. 
+!    Values can be separated by commas or spaces.
+!(5) space allowed before/after '='
+
+&CORE
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+! Core (mandatory) parameters; no defaults
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+! Pre-processing option. Useful for checking grid errors etc. Need to use 1 
+! core only for compute (plus necessary scribe cores). Under scribe I/O, the 
+! code (the scribe part) will hang but the outputs will be there. Just kill 
+! the job.
+  ipre = 0 !Pre-processor flag (1: on; 0: off)
+
+! Baroclinic/barotropic option. If ibc=0 (baroclinic model), ibtp is not used.
+  ibc = 0 !Baroclinic option
+  ibtp = 1 
+
+  rnday = 30 !total run time in days
+  dt = 100. !Time step in sec
+
+! Grid for WWM (USE_WWM)
+  msc2 = 24     !same as msc in .nml ... for consitency check between SCHISM and WWM
+  mdc2 = 30     !same as mdc in .nml
+
+! Define # of tracers in tracer modules (if enabled)
+  ntracer_gen = 2 !user defined module (USE_GEN)
+  ntracer_age = 4 !age calculation (USE_AGE). Must be =2*N where N is # of age tracers
+  sed_class = 5 !SED3D (USE_SED)
+  eco_class = 27 !EcoSim (USE_ECO): must be between [25,60]
+
+! Global output controls
+  nspool = 36 !output step spool
+  ihfskip = 864 !stack spool; every ihfskip steps will be put into 1_*, 2_*, etc...
+/
+
+&OPT
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+! Optional parameters. The values shown below are default unless otherwise noted
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+!-----------------------------------------------------------------------
+! Optional offline partitioning using NO_PARMETIS. If on, need partition.prop (like global_to_local.prop).
+!-----------------------------------------------------------------------
+
+
+!-----------------------------------------------------------------------
+! 2nd pre-proc flag. If ipre2/=0, code will output some diagnotic outputs and stop.
+! These outputs are: drag coefficients (Cdp)
+!-----------------------------------------------------------------------
+  ipre2 = 0
+
+!-----------------------------------------------------------------------
+! Option to only solve tracer transport (and bypass most b-tropic solver)
+! Usage: turn the flag on ('1' or '2') and turn on your tracer modules and make sure
+! (1) your inputs are consistent with the original hydro-only run (with additional tracers
+! of course); (2) hydro-only run results are in hydro_out/schout*.nc, which must
+! have 'hvel_side' (not normal hvel), 'elev', 'diffusivity', 'temp_elem', 'salt_elem' (for
+! new scribe outputs, use corresponding files/var names); (3) dt above is 
+! multiple of _output_ step used in the original hydro-only run 
+! (as found in in hydro_out/schout*.nc); e.g. dt = 1 hour. (4). When itransport_only=2,
+! additional variables ('totalSuspendedLoad','sedBedStress') are needed.
+! Hotstart should work also, but you'd probably not use an aggressively large dt especially
+! when air-sea exchange is involved.
+!-----------------------------------------------------------------------
+  itransport_only = 0
+
+!-----------------------------------------------------------------------
+! Option to add self-attracting and loading tide (SAL) into tidal potential 
+! (usually for basin-scale applications). 
+! If iloadtide=0, no SAL.
+! If iloadtide=1, needs inputs: loadtide_[FREQ].gr3, 
+! where [FREQ] are freq names (shared with tidal potential, in upper cases) 
+! and the _two_ 'depths' inside are amplitude (m) and phases (degrees behind GMT), 
+! interpolated from global tide model (e.g. FES2014). In this option, SAL is 
+! lumped into tidal potential so it shares some parameters with tidal potential 
+! in bctides.in (cut-off depth, frequencies).
+! If iloadtide=2 or 3, use a simple scaling for gravity approach (in this option,
+! SAL is applied everywhere and does not share parameters with tidal potential).
+! For iloadtide=2, a const scaling (1-loadtide_coef) is used; for iloadtide=3, the scaling is 
+! dependent on depth (Stepanov & Hughes 2004) with max of loadtide_coef.
+!-----------------------------------------------------------------------
+  iloadtide = 0
+  loadtide_coef = 0.1 !only used if iloadtide=2,3 (for '3', the default should be 0.12)
+
+! Starting time
+  start_year = 2000 !int
+  start_month = 1 !int
+  start_day = 1 !int
+  start_hour = 0 !double
+  utc_start = 8 !double
+
+!-----------------------------------------------------------------------
+! Coordinate option: 1: Cartesian; 2: lon/lat (hgrid.gr3=hgrid.ll in this case,
+! and orientation of element is outward of earth)
+! Notes for lon/lat: make sure hgrid.ll and grid in sflux are consistent in 
+! longitude range!
+!-----------------------------------------------------------------------
+  ics = 1 !Coordinate option
+
+!-----------------------------------------------------------------------
+! Hotstart option. 0: cold start; 1: hotstart with time reset to 0; 2: 
+! continue from the step in hotstart.nc
+!-----------------------------------------------------------------------
+  ihot = 0
+
+!-----------------------------------------------------------------------
+! Equation of State type used
+! ieos_type=0: UNESCO 1980 (nonlinear); =1: linear function of T ONLY, i.e. 
+! \rho=eos_b+eos_a*T, where eos_a<=0 in kg/m^3/C
+!-----------------------------------------------------------------------
+  ieos_type = 0
+  ieos_pres = 0 !used only if ieos_type=0. 0: without pressure effects 
+  eos_a = -0.1 !needed if ieos_type=1; should be <=0 
+  eos_b = 1001. !needed if ieos_type=1
+
+!-----------------------------------------------------------------------
+! Main ramp option
+!-----------------------------------------------------------------------
+!  nramp = 1 !ramp-up option (1: on; 0: off)
+  dramp = 1. !ramp-up period in days for b.c. etc (no ramp-up if <=0)
+
+!  nrampbc = 0 !ramp-up flag for baroclinic force
+  drampbc = 0. !ramp-up period in days for baroclinic force
+
+!-----------------------------------------------------------------------
+! Method for momentum advection. 0: ELM; 1: upwind (not quite working yet)
+!-----------------------------------------------------------------------
+  iupwind_mom = 0
+
+!-----------------------------------------------------------------------
+! Methods for computing velocity at nodes. 
+! If indvel=0, conformal linear shape function is used; if indvel=1, averaging method is used.
+! For indvel=0, a stabilization method is needed (see below). 
+!-----------------------------------------------------------------------
+  indvel = 0 
+ 
+!-----------------------------------------------------------------------
+! 2 stabilization methods, mostly for indvel=0.
+! (1) Horizontal viscosity option. ihorcon=0: no viscosity is used; =1: Lapacian;
+! =2: bi-harmonic. If ihorcon=1, horizontal viscosity _coefficient_ (<=1/8, related
+! to diffusion number) is given in hvis_coef0, and the diffusion # 
+! is problem dependent; [0.001-1/8] seems to work well.
+! If ihorcon=2, diffusion number is given by hvis_coef0 (<=0.025).
+! If indvel=1, no horizontal viscosity is needed. 
+! (2) Shapiro filter (see below)
+!
+! For non-eddying regime applications (nearshore, estuary, river), an easiest option is: 
+!  indvel=0, ishapiro=1 (shapiro0=0.5), ihorcon=inter_mom=0.
+! For applications that include eddying regime, refer to the manual.
+!-----------------------------------------------------------------------
+  ihorcon = 0
+  hvis_coef0 = 0.025 !const. diffusion # if ihorcon/=0; <=0.025 for ihorcon=2, <=0.125 for ihorcon=1
+!  cdh = 0.01 !needed only if ihorcon/=0; land friction coefficient - not active yet
+
+!-----------------------------------------------------------------------
+! 2nd stabilization method via Shapiro filter. This should normally be used 
+! if indvel=0. To transition between eddying/non-eddying regimes, use
+! indvel=0, ihorcon/=0, and ishapiro=-1 (requiring shapiro.gr3) or 2 (Smagorinsky-like filter).
+!-----------------------------------------------------------------------
+  ishapiro = 1 !options
+  niter_shap = 1 !needed if ishapiro/=0: # of iterations with Shapiro filter
+  !shapiro0: Shapiro filter strength, needed only if ishapiro=1 or 2 
+  !If ishapiro=1, shapiro0 is the filter strength (max is 0.5).
+  !If ishapiro=2, shapiro0 is the coefficient in tanh(). Experiences so far suggest ~1.e3
+  !If ishapiro=-1, this is not used as the filter strength is read in from shapiro.gr3
+  shapiro0 = 0.5 
+
+!-----------------------------------------------------------------------
+! Implicitness factor (0.5<thetai<=1).
+!-----------------------------------------------------------------------
+  thetai = 0.6 
+
+!-----------------------------------------------------------------------
+! If WWM is used, set coupling/decoupling flag 'icou_elfe_wwm'. 
+! No effects if USE_WWM is distabled in Makefile.
+! Note that all these parameters must be present in this file (even when not used).
+!       0: no elevation and no currents in wwm, no wave force in SCHISM (but wave turbulecne, WBL etc are still in SCHISM);
+!       1: full coupled (elevation, vel, and wind are all passed to WWM); 
+!       2: elevation and currents in wwm, no wave force in SCHISM;  
+!       3: no elevation and no currents in wwm, wave force in SCHISM;
+!       4: elevation but no currents in wwm, wave force in SCHISM;
+!       5: elevation but no currents in wwm, no wave force in SCHISM;
+!       6: no elevation but currents in wwm, wave force in SCHISM;
+!       7: no elevation but currents in wwm, no wave force in SCHISM;
+! If WWM is used and RADFLAG='VOR' in wwminput.nml, some parameters (fwvor_*)
+! are added to take into account (1) or not (0) the different terms in vortex
+! force expression.
+! To get SCHISM-only result withno feedback from WWM, compile without WWM.
+!-----------------------------------------------------------------------
+  icou_elfe_wwm = 0 
+  nstep_wwm = 1  !call WWM every this many time steps 
+  iwbl = 0 !wave boundary layer formulation (used only if USE_WMM and 
+           !icou_elfe_wwm/=0 and nchi=1. If icou_elfe_wwm=0, set iwbl=0): 
+           !1-modified Grant-Madsen formulation; 2-Soulsby (1997)
+  hmin_radstress = 1. !min. total water depth used only in radiation stress calculation [m]
+!  nrampwafo = 0       !ramp-up option for the wave forces (1: on; 0: off)
+  drampwafo = 0.      !ramp-up period in days for the wave forces (no ramp-up if <=0)
+  turbinj = 0.15      !% of depth-induced wave breaking energy injected in turbulence 
+                      !(default: 0.15 (15%), as proposed by Feddersen, 2012)
+  turbinjds = 1.0     !% of wave energy dissipated through whitecapping injected in turbulence 
+                      !(default: 1 (100%), as proposed by Paskyabi et al. 2012)
+  alphaw = 0.5        !for itur=4 : scaling parameter for the surface roughness z0s = alphaw*Hm0. 
+                      !If negative z0s = abs(alphaw) e.g. z0s=0.2 m (Feddersen and Trowbridge, 2005)
+                         ! Vortex Force terms (off/on:0/1)
+  fwvor_advxy_stokes = 1 !  --> Stokes drift advection (xy), Coriolis
+  fwvor_advz_stokes = 1  !  --> Stokes drift advection (z) , Coriolis
+  fwvor_gradpress = 1    !  --> Pressure term
+  fwvor_breaking = 1     !  --> Wave breaking
+  fwvor_streaming = 1    !  --> Wave streaming (works with iwbl /= 0)
+  fwvor_wveg = 0         !  --> Wave dissipation by vegetation acceleration term
+  fwvor_wveg_NL = 0      !  --> Non linear intrawave vegetation force (see Dean and Bender, 2006 or van Rooijen et al., 2016 for details) 
+  cur_wwm = 0            ! Coupling current in WWM
+                         ! 0: surface layer current
+                         ! 1: depth-averaged current
+                         ! 2: computed according to Kirby and Chen (1989)
+  wafo_obcramp = 0       ! Ramp on wave forces at open boundary (1: on / 0: off)
+                         !  --> this option requires the input file wafo_ramp.gr3
+                         !      which defines the ramp value (between 0 and 1)
+                         !      at nodes over the whole domain
+
+!-----------------------------------------------------------------------
+! Bed deformation option (0: off; 1: vertical deformation only; 2: 3D bed deformation). 
+! If imm=1, bdef.gr3 is needed; if imm=2, user needs to update depth info etc
+! in the code (not working for ics=2 yet).
+!-----------------------------------------------------------------------
+  imm = 0
+  ibdef = 10 !needed if imm=1; # of steps used in deformation
+
+!-----------------------------------------------------------------------
+! Reference latitude for beta-plane approximation when ncor=1 (not used if ics=2)
+!-----------------------------------------------------------------------
+  slam0 = -124  !lon - not really used
+  sfea0 = 45 !lat
+
+!-----------------------------------------------------------------------
+! Option to deal with under resolution near steep slopes in deeper depths
+! 0: use h[12,_bcc below; /=0: use hw_* below
+!-----------------------------------------------------------------------
+  iunder_deep = 0 
+
+!-----------------------------------------------------------------------
+! Baroclinicity calculation in off/nearshore with iunder_deep=ibc=0. 
+! The 'below-bottom' gradient
+! is zeroed out if h>=h2_bcc (i.e. like Z) or uses const extrap
+! (i.e. like terrain-following) if h<=h1_bcc(<h2_bcc) (and linear
+! transition in between based on local depth)
+!-----------------------------------------------------------------------
+  h1_bcc = 50. ![m]
+  h2_bcc = 100. ![m]; >h1_bcc
+
+!-----------------------------------------------------------------------
+! Hannah-Wright-like ratio & depth used to account for under-resolution
+! in a b-clinic model. Used only if ibc=0 and iunder_deep/=0.
+! The b-clinic force at prism centers is calculated with a reconstruction
+! method in horizontal that has a stencil of an element and its adjacent elements.
+! If the depths change is too much between the elem and its adjacent elem
+! the under-resolution occurs (with steep bottom slope) and b-clinic force
+! needs to be zeroed out below the (higher) bottom, specifically, if
+! max(2 depths)>=hw_depth and abs(diff(2 depths))>=hw_ratio*max(2 depths).
+!-----------------------------------------------------------------------
+  hw_depth = 1.e6 !threshold depth in [m]
+  hw_ratio = 0.5 !ratio
+
+!-----------------------------------------------------------------------
+! Hydraulic model option. If ihydraulics/=0, hydraulics.in 
+! is required. This option cannot be used with non-hydrostatic model.
+!-----------------------------------------------------------------------
+  ihydraulics = 0
+
+!-----------------------------------------------------------------------
+! Point sources/sinks option (0: no; 1: ASCII inputs; -1: netcdf). 
+! If =1, needs source_sink.in (list of elements),
+! vsource,th, vsink.th, and msource.th (the source/sink values must be single precision). 
+! If =-1, all info is in source.nc
+! and each type of volume/mass source/sink can have its own time step and
+! # of records.
+! Source/sinks can be specified at an elem more
+! than once, and the code will accumulate the volumes; for mass conc, 
+! values are applied at _net_ source elem (no summation for conc).
+
+! meth_sink: options to treat sinks @ dry elem: no speacial treatment if meth_sink=0;
+! zero out sink if the elem is dry if meth_sink=1
+!-----------------------------------------------------------------------
+  if_source = 0
+!  nramp_ss = 1 !needed if if_source=1; ramp-up flag for source/sinks
+  dramp_ss = 2 !needed if if_source/=0; ramp-up period in days for source/sinks (no ramp-up if <=0)
+  meth_sink = 1 !options to treat sinks @ dry elem
+
+! Specify vertical level to inject source concentration for each tracer _module_.
+! Code will extrapolate below bottom/above surface unless level 0 is specified. In
+! the latter case, the incoming tracer mass will be distributed across all vertical layers - 
+! this approach works best in shallow waters. 
+  lev_tr_source(1) = -9 !T
+  lev_tr_source(2) = -9 !S
+  lev_tr_source(3) = -9 !GEN
+  lev_tr_source(4) = -9 !AGE
+  lev_tr_source(5) = -9 !SED3D
+  lev_tr_source(6) = -9 !EcoSim
+  lev_tr_source(7) = -9 !ICM
+  lev_tr_source(8) = -9 !CoSINE
+  lev_tr_source(9) = -9 !Feco
+  lev_tr_source(10) = -9 !TIMOR
+  lev_tr_source(11) = -9 !FABM
+  lev_tr_source(12) = -9 !DVD
+
+!-----------------------------------------------------------------------
+! Horizontal diffusivity option. if ihdif=1, horizontal diffusivity is given in hdif.gr3
+!-----------------------------------------------------------------------
+  ihdif = 0 
+
+!-----------------------------------------------------------------------
+! Bottom friction. 
+!           nchi=0: drag coefficients specified in drag.gr3; nchi=-1: Manning's 
+!           formulation (even for 3D prisms) with n specified in manning.gr3. 
+!           nchi=1: bottom roughness (in meters) specified in rough.gr3 (and in this case, negative
+!           or 0 depths in rough.gr3 indicate time-independent Cd, not roughness!).
+!           Cd is calculated using the log law, when dzb>=dzb_min; when dzb<dzb_min,
+!           Cd=Cdmax, where Cdmax=Cd(dzb=dzb_min).
+!           If iwbl/=0, nchi must =1.
+!-----------------------------------------------------------------------
+  nchi = 0 
+  dzb_min = 0.5 !needed if nchi=1; min. bottom boundary layer thickness [m].
+!  dzb_decay = 0. !needed if nchi=1; a decay const. [-]. should =0
+  hmin_man = 1. !needed if nchi=-1: min. depth in Manning's formulation [m]
+
+!-----------------------------------------------------------------------
+! Coriolis. If ncor=-1, specify "rlatitude" (in degrees); if ncor=0,
+! specify Coriolis parameter in "coricoef"; if ncor=1, model uses
+! lat/lon in hgrid.ll for beta-plane approximation if ics=1, and in this case,
+! the latitude specified in CPP projection ('sfea0') is used. If ncor=1 and ics=2,
+! Coriolis is calculated from local latitude, and 'sfea0' is not used.
+!-----------------------------------------------------------------------
+  ncor = 0 !should usually be 1 if ics=2
+  rlatitude = 46 !if ncor=-1
+  coricoef = 0 !if ncor=0
+
+!-----------------------------------------------------------------------
+! Elevation initial condition flag for cold start only. If ic_elev=1, elev.ic (in *.gr3 format) is needed
+! to specify the initial elevations; otherwise elevation is initialized to 0 everywhere 
+!-----------------------------------------------------------------------
+  ic_elev = 0
+
+!-----------------------------------------------------------------------
+! Elevation boundary condition ramp-up flag. =0: ramp up from 0; =1: ramp up from
+! elev. values read in from elev.ic or hotstart.nc - if neither is present, from 0.
+! This flag is mainly used to start the simulation from non-zero elev.
+! The ramp-up period is 'dramp' (so if dramp=0, full strength is applied).
+!-----------------------------------------------------------------------
+  nramp_elev = 0
+
+!-----------------------------------------------------------------------
+! Optional inverse barometric effects on the elev. b.c.
+! If inv_atm_bnd=1, the elev.'s at boundary are corrected by the difference
+! between the actual atmos. pressure and a reference pressure (prmsl_ref below)
+!-----------------------------------------------------------------------
+  inv_atm_bnd = 0 !0: off; 1: on
+  prmsl_ref = 101325. !reference atmos. pressure on bnd [Pa]
+
+!-----------------------------------------------------------------------
+! Initial condition for T,S. This value only matters for ihot=0 (cold start).
+! If flag_ic(1:2)=1, the initial T,S field is read in from temp.ic and salt.ic (horizontally varying).
+! If flag_ic(1:2)=2, the initial T,S field is read in from ts.ic (vertical varying).
+! If ihot=0 && flag_ic(1:2)=2 || ibcc_mean=1, ts.ic is used for removing mean density profile.
+! flag_ic(1) must =flag_ic(2)
+!-----------------------------------------------------------------------
+  flag_ic(1) = 1 !T
+  flag_ic(2) = 1 !S
+
+! initial conditions for other tracers.
+! 1: needs inputs [MOD]_hvar_[1,2,...].ic ('1...' is tracer id); format of each file is similar to salt.ic;
+!    i.e. horizontally varying i.c. is used for each tracer.
+! 2: needs [MOD]_vvar_[1,2,...].ic. Format of each file (for each tracer in tis MOD) is similar to ts.ic
+!    (i.e. level #, z-coord., tracer value). Verically varying i.c. is used for each tracer.
+! 0: model sets own i.c. (EcoSim; TIMOR)
+  flag_ic(3) = 1 !GEN (user defined module)
+!  flag_ic(4) = 1 !Age i.c. flag set inside code
+  flag_ic(5) = 1 !SED3D
+  flag_ic(6) = 1 !EcoSim
+  flag_ic(7) = 1 !ICM
+  flag_ic(8) = 1 !CoSINE
+  flag_ic(9) = 1 !FIB
+  flag_ic(10) = 1 !TIMOR
+  flag_ic(11) = 1 !FABM
+  flag_ic(12) = 0 !DVD (must=0)
+
+!-----------------------------------------------------------------------
+! Settling vel [m/s] for GEN module (positive downward)
+!-----------------------------------------------------------------------
+  gen_wsett = 0 !1.e-4
+
+!-----------------------------------------------------------------------
+! Mean T,S profile option. If ibcc_mean=1 (or ihot=0 and flag_ic(1)=2), mean profile
+! is read in from ts.ic, and will be removed when calculating baroclinic force.
+! No ts.ic is needed if ibcc_mean=0.
+!-----------------------------------------------------------------------
+  ibcc_mean = 0 
+
+!-----------------------------------------------------------------------
+! Max. horizontal velocity magnitude, used mainly to prevent problem in 
+! bulk aerodynamic module
+!-----------------------------------------------------------------------
+  rmaxvel = 5.
+
+!-----------------------------------------------------------------------
+!  Following parameters control backtracking
+!-----------------------------------------------------------------------
+!-----------------------------------------------------------------------
+!  min. vel for invoking btrack and for abnormal exit in quicksearch
+!-----------------------------------------------------------------------
+  velmin_btrack = 1.e-4
+!-----------------------------------------------------------------------
+! Nudging factors for starting side/node - add noise to avoid underflow
+! The starting location is nudged to: old*(1-btrack_nudge)+btrack_nudge*centroid
+! Suggested value: btrack_nudge=9.013e-3
+!-----------------------------------------------------------------------
+  btrack_nudge= 9.013e-3 
+
+!-----------------------------------------------------------------------
+! Behavior when trajectory hits open bnd. If ibtrack_openbnd=0, slide with
+! tangential vel; otherwise, stop and exit btrack 
+!-----------------------------------------------------------------------
+!  ibtrack_openbnd = 1 !hardwired
+
+!-----------------------------------------------------------------------
+! Wetting and drying. 
+! - if ihhat=1, \hat{H} is made non-negative to enhance robustness near 
+! wetting and drying; if ihhat=0, no retriction is imposed for this quantity. 
+! - inunfl=0 is used for normal cases and inunfl=1 is used for more accurate wetting
+! and drying if grid resolution is sufficiently fine.
+! - if shorewafo=1, we impose radiation stress R_s = g*grad(eta) (balance between radiation stress
+! gradients and the barotropic gradients) at the numerical shoreline (boundary between
+! dry and wet elements). This option ensures that the shallow depth in dry elements does not
+! create unphysical and very high wave forces at the shoreline (advised for morphodynamics runs).
+!-----------------------------------------------------------------------
+  ihhat = 1 
+  inunfl = 0
+  h0 = 0.01     !min. water depth for wetting/drying [m]
+  shorewafo = 0 !Matters only if USE_WWM
+
+!-----------------------------------------------------------------------
+! Solver options
+! USE_PETSC controls the solver type. If it's diabled, the default JCG 
+! solver is used. If it's enabled, use PetSc lib. Some of the parameters
+! have different meanings under these 2 options. Also with PetSc one can
+! use cmd line options to choose solver etc.
+!-----------------------------------------------------------------------
+  moitn0 = 50 !output spool for solver info; used only with JCG
+  mxitn0 = 1500 !max. iteration allowed
+  rtol0 = 1.e-12 !error tolerance
+
+!-----------------------------------------------------------------------
+! Advection (ELM) option. If nadv=1, backtracking is done using Euler method; 
+! nadv=2, using 2nd order Runge-Kutta; if nadv=0, advection in momentum 
+! is turned off/on in adv.gr3 (the depths=0,1, or 2 also control methods 
+! in backtracking as above). dtb_max/min are the max/min steps allowed -
+! actual step is calculated adaptively based on local gradient.
+!-----------------------------------------------------------------------
+  nadv = 1
+  dtb_max = 30. !in sec
+  dtb_min = 10.
+
+!-----------------------------------------------------------------------
+! If inter_mom=0, linear interpolation is used for velocity at foot of char. line.
+! If inter_mom=1 or -1, Kriging is used, and the choice of covariance function is
+! specified in 'kr_co'. If inter_mom=1, Kriging is applied to whole domain;
+! if inter_mom=-1, the regions where Kriging is used is specified in krvel.gr3 
+! (depth=0: no kriging; depth=1: with kriging). 
+!-----------------------------------------------------------------------
+  inter_mom = 0 
+  kr_co = 1 !not used if inter_mom=0
+
+!-----------------------------------------------------------------------
+! Tracer transport method. TVD or WENO method requires tvd.prop.
+! TVD or WENO method is used on an element/prism if the total depth (at all nodes of the elem.)>=h_tvd and the flag in
+! tvd.prop = 1 for the elem; otherwise upwind is used for efficiency. 
+! itr_met=3 (horizontal TVD) or 4 (horizontal WENO): implicit TVD in the vertical dimension. 
+! Also if itr_met==3 and h_tvd>=1.e5, some parts of the code are bypassed for efficiency.
+! Controls for WENO are not yet in place.
+!-----------------------------------------------------------------------
+  itr_met = 3 
+  h_tvd = 5. !cut-off depth (m) 
+  !If itr_met=3 or 4, need the following 2 tolerances of convergence. The convergence
+  !is achieved when sqrt[\sum_i(T_i^s+1-T_i^s)^2]<=eps1_tvd_imp*sqrt[\sum_i(T_i^s)^2]+eps2_tvd_imp
+  eps1_tvd_imp = 1.e-4 !suggested value is 1.e-4, but for large suspended load, need to use a smaller value (e.g. 1.e-9)
+  eps2_tvd_imp = 1.e-14  
+
+  !Optional hybridized ELM transport for efficiency
+  ielm_transport = 0 !1: turn on 
+  max_subcyc = 10 !used only if ielm_transport/=0. Max # of subcycling per time step in transport allowed
+
+  !if itr_met = 4, the following parameters are needed
+  !if itr_met=4 and ipre=1, diagnostic outputs are generated for weno accuracy and stencil quality, 
+  !  see subroutine weno_diag in src/Hydro/misc_subs.F90 for details
+  ip_weno = 2   !order of accuracy: 0- upwind; 1- linear polynomial, 2nd order; 2- quadratic polynomial, 3rd order
+  courant_weno=0.5 !Courant number for weno transport
+  nquad = 2  !number of quad points on each side, nquad= 1 or 2
+  ntd_weno = 1 !order of temporal discretization: (1) Euler (default); (3): 3rd-order Runge-Kutta (only for benchmarking)
+  epsilon1 = 1.e-15   !coefficient for 2nd order weno smoother (larger values are more prone to numerical dispersion)
+  epsilon2 = 1.e-10  !1st coefficient for 3rd order weno smoother (larger values are more prone to numerical dispersion
+                     !, 1.e-10 should be fairly safe, recommended values: 1.e-8 ~ 1.e-6 (from the real applications so far)
+  i_prtnftl_weno = 0 !option for writing nonfatal errors on invalid temp. or salinity for density: (0) off; (1) on.
+
+  !Inactive at the moment:
+  epsilon3 = 1.e-25  !2nd coefficient for 3rd order weno smoother (inactive at the moment)
+  !Elad filter has not been implemented yet; preliminary tests showed it might not be necessary
+  ielad_weno = 0      !ielad, if ielad=1, use ELAD method to suppress dispersion (inactive at the moment)
+  small_elad = 1.e-4  !small (inactive at the moment)
+
+!-----------------------------------------------------------------------
+! Atmos. option. nws=3 is reserved for coupling with atmospheric model.
+! If nws=0, no atmos. forcing is applied. If nws=1, atmos.
+! variables are read in from wind.th. If nws=2, atmos. variables are
+! read in from sflux_ files.
+! If nws=4, ascii format is used for wind and atmos. pressure at each node (see source code).
+! If nws=-1 (requires USE_PAHM), use Holland parametric wind model (barotropic only with wind and atmos. pressure).
+!  In this case, the Holland model is called every step so wtiminc is not used. An extra 
+!  input is needed: hurricane-track.dat
+!
+! Stress calculation:
+! If nws=2, ihconsv=1 and iwind_form=0, the stress is calculated from heat exchange
+! routine; 
+! Otherwise if iwind_form=-1, the stress is calculated from Pond & Pichard formulation;
+! if iwind_form=1, Hwang (2018) formulation (Cd tapers off at high wind).
+! If WWM is enabled and icou_elfe_wwm>0 and iwind_form=-2 or -3, stress is overwritten by WWM:
+! If iwind_form=-2, stress=rho_air*ufric^2; scaled by rho_water
+! If iwind_form=-3, the stress is calculated according to the param. of Donelan et al. (1993) based on the wave age.
+! In all cases, if USE_ICE the stress in ice-covered portion is calculated by ICE routine.
+!-----------------------------------------------------------------------
+  nws = 0 
+  wtiminc = 150. !time step for atmos. forcing. Default: same as dt
+!  nrampwind = 1 !ramp-up option for atmos. forcing
+  drampwind = 1. !ramp-up period in days for wind (no ramp-up if <=0)
+  iwindoff = 0 !needed only if nws/=0; '1': needs windfactor.gr3
+  iwind_form = 1 !needed if nws/=0
+  !If IMPOSE_NET_FLUX is on and nws=2, read in net _surface_ heat flux as var 'dlwrf' 
+  !(Downward Long Wave) in sflux_rad (solar radiation is still used separately),
+  !and if PREC_EVAP is on, also read in net P-E as 'prate' (Surface Precipitation Rate) in sflux_prc.
+
+!-----------------------------------------------------------------------
+! Heat and salt exchange using Zeng's (default) or Fairall (USE_BULK_FAIRALL) scheme.
+! isconsv=1 needs ihconsv=1; ihconsv=1 needs nws=2.
+! If isconsv=1, need to compile with precip/evap module turned on.
+! If at least one of ihconsv and isconsv is set to 1,
+! locally turning off air-sea exchange in shallow waters (<0.25 m) is recommended.
+! Options for locally turning off heat/salt exchange are specified by
+! i_hmin_airsea_ex, i_hmin_salt_ex respectively:
+!   0: not turned off locally;
+!   1: locally turned off, based on grid depth, i.e.,
+!      off when the local grid depth (z in hgrid) < hmin_airsea_ex or hmin_salt_ex;
+!   2 (recommended): locally turned off, based on local total water depth, i.e.,
+!      off when the local total water depth < hmin_airsea_ex or hmin_salt_ex,
+! hmin_airsea_ex, hmin_salt_ex:
+!   0.2 m is recommended for both "1" and "2"
+!-----------------------------------------------------------------------
+  ihconsv = 0 !heat exchange option
+  isconsv = 0 !evaporation/precipitation model
+  i_hmin_airsea_ex = 2 ! no effect if ihconsv=0 
+  hmin_airsea_ex = 0.2 ![m], no effect if ihconsv=0 
+  i_hmin_salt_ex = 2 ! no effect if isconsv=0 
+  hmin_salt_ex = 0.2 ![m], no effect if isconsv=0 
+  iprecip_off_bnd = 0 !if /=0, precip will be turned off near land bnd
+
+!-----------------------------------------------------------------------
+! Turbulence closure.
+!-----------------------------------------------------------------------
+  itur = 3 !Default: 0
+  dfv0 = 1.e-2 !needed if itur=0
+  dfh0 = 1.e-4 !needed if itur=0
+  mid = 'KL' !needed if itur=3,5. Use KE if itur=5 
+  stab = 'KC' !needed if itur=3 or 5. Use 'GA' if turb_met='MY'; otherwise use 'KC'. 
+  xlsc0 = 0.1 !needed if itur=3 or 5. Scale for surface & bottom mixing length (>0)
+
+!-----------------------------------------------------------------------
+! Sponge layer for elevation and vel.
+! If inu_elev=0, no relaxation is applied to elev.
+! If inu_elev=1, relax. constants (in 1/sec, i.e. relax*dt<=1) are specified in elev_nudge.gr3
+! and applied to eta=0 (thus a depth=0 means no relaxation).
+! Similarly for inu_uv (with input uv_nudge.gr3)
+!-----------------------------------------------------------------------
+  inu_elev = 0
+  inu_uv = 0
+
+!-----------------------------------------------------------------------
+! Nudging options for tracers. If inu_[MOD]=0, no nudging is used. If inu_tr=1,
+! nudge to initial condition according to relaxation constants specified.
+! If inu_tr=2, nudge to values in [MOD]_nu.nc (with step 'step_nu_tr').
+! The relaxation constants = [horizontal relax (specified in [MOD]_nudge.gr3) + or x
+! vertical relax] times dt, where vertical relax is a linear function of 
+! vnh[1,2] and vnf[1,2], and [MOD] are tracer model names. 'nu_sum_mult' decides
+! '+' or 'x' in the calculation of final relax.
+! Code will ignore junk values (<=-99) inside [MOD]_nu.nc, so 1 way to avoid 
+! nudging for a tracer is to set its nudged values to -9999
+!-----------------------------------------------------------------------
+  inu_tr(1) = 0 !T
+  inu_tr(2) = 0 !S
+  inu_tr(3) = 0 !GEN
+  inu_tr(4) = 0 !Age
+  inu_tr(5) = 0 !SED3D
+  inu_tr(6) = 0 !EcoSim 
+  inu_tr(7) = 0 !ICM 
+  inu_tr(8) = 0 !CoSINE 
+  inu_tr(9) = 0 !FIB
+  inu_tr(10) = 0 !TIMOR 
+  inu_tr(11) = 0 !FABM 
+  inu_tr(12) = 0 !DVD (must=0)
+
+  nu_sum_mult=1 !1: final relax is sum of horizontal&vertical; 2: product
+  vnh1 = 400 !vertical nudging depth 1
+  vnf1 = 0. !vertical relax \in [0,1]
+  vnh2 = 500 !vertical nudging depth 2 (must >vnh1)
+  vnf2 = 0. !vertical relax
+
+  step_nu_tr = 86400. !time step [sec] in all [MOD]_nu.nc (for inu_[MOD]=2)
+
+!-----------------------------------------------------------------------
+! Cut-off depth for cubic spline interpolation near bottom when computing horizontal gradients
+! e.g. using hgrad_nodes() (radiation stress, and gradients of qnon and qhat in non-hydro model). 
+! If depth > h_bcc1 ('deep'),
+! a min. (e.g. max bottom z-cor for the element) is imposed in the spline and so a more
+! conservative method is used without extrapolation beyond bottom; 
+! otherwise constant extrapolation below bottom is used.
+!-----------------------------------------------------------------------
+  h_bcc1 = 100. !h_bcc1
+
+!-----------------------------------------------------------------------
+! Dimensioning parameters for inter-subdomain btrack. 
+! If error occurs like 'bktrk_subs: overflow' or 'MAIN: nbtrk > mxnbt'
+! gradually increasing these will solve the problem
+!-----------------------------------------------------------------------
+  s1_mxnbt = 0.5
+  s2_mxnbt = 3.5
+
+!-----------------------------------------------------------------------
+! Flag for harmonic analysis for elevation. If used , need to turn on USE_HA
+! in Makefile, and input harm.in. Otherwise set it to 0. Hotstart ihot=2 is not working with HA
+! Outputs are harme_* and use combine_outHA to combine.
+!-----------------------------------------------------------------------
+  iharind = 0
+
+!-----------------------------------------------------------------------
+! Conservation check option. If iflux=1, some fluxes are computed
+! in regions specified in fluxflag.prop (regional number from -1 to an arbitrary integer).
+! in output flux.out, positive means flux from region n to region n-1 (n>=1)
+! output file flux.dat need to be saved before continuing hotstart
+!-----------------------------------------------------------------------
+  iflux = 0
+  iflux_out_format = 0 !0: basic flux output; 1: more detailed outputs
+
+!-----------------------------------------------------------------------
+! Test flags for debugging. These flags should be turned off normally.
+!-----------------------------------------------------------------------
+! Williamson test #5 (zonal flow over an isolated mount); if
+! on, ics must =2
+!-----------------------------------------------------------------------
+  izonal5 = 0 !"0" - no test; otherwise on
+
+!-----------------------------------------------------------------------
+! Rotating Gausshill test with stratified T,S (1: on; 0: off)
+! Surface T,S read in from *.ic; code generates stratification
+!-----------------------------------------------------------------------
+  ibtrack_test = 0
+
+!-----------------------------------------------------------------------
+! Rouse profile test (1: on; 0: off)
+! If on, must turn on USE_TIMOR
+!-----------------------------------------------------------------------
+  irouse_test = 0
+
+!-----------------------------------------------------------------------
+! Flag to choose FIB model for bacteria decay (used with USE_FIB)
+! flag_fib = 1 - Constant decay rate (/day) in .gr3 format
+!                (kkfib_1.gr3 and kkfib_2.gr3)
+! flag_fib = 2 - Decay rate computed from Canteras et al., 1995
+! flag_fib = 3 - Decay rate computed from Servais et al., 2007
+!----------------------------------------------------------------------
+  flag_fib = 1
+
+!----------------------------------------------------------------------
+! Marsh model parameters (only if USE_MARSH is on)
+!----------------------------------------------------------------------
+  slr_rate = 120. !sea-level rise rate in mm/yr, times morphological acceleration if used
+
+!----------------------------------------------------------------------
+! Vegetation model
+! If isav=1, need 4 extra inputs: (1) sav_D.gr3 (depth is stem diameter in meters);
+! (2) sav_N.gr3 (depth is # of stems per m^2);
+! (3) sav_h.gr3 (height of canopy in meters);
+! (4) sav_cd.gr3 (drag coefficient).									
+! If one of these depths=0 at a node, the code will set all to 0. 
+! If USE_MARSH is on and isav=1, all .gr3 must have constant depths!
+!----------------------------------------------------------------------
+  isav = 0 !on/off flag
+
+!----------------------------------------------------------------------
+! Coupling step with ICE module.
+!----------------------------------------------------------------------
+  nstep_ice = 1 !call ice module every nstep_ice steps of SCHISM
+
+
+!----------------------------------------------------------------------
+! Specify level #'s if age module is invoked (USE_AGE), for 1st half of tracers only
+!----------------------------------------------------------------------
+  level_age = 9, -999 !default: -999 (all levels)
+
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+! Physical constants
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+!-----------------------------------------------------------------------
+! Earth's radii at pole and equator (to define an ellipsoid)
+!-----------------------------------------------------------------------
+  rearth_pole = 6378206.4
+  rearth_eq = 6378206.4
+
+!-----------------------------------------------------------------------
+! Specific heat of water (C_p) in J/kg/K
+!-----------------------------------------------------------------------
+  shw = 4184.d0
+
+!-----------------------------------------------------------------------
+! Reference water density for Boussinesq approximation
+!-----------------------------------------------------------------------
+  rho0 = 1000.d0 !kg/m^3
+
+!-----------------------------------------------------------------------
+! Fraction of vertical flux closure adjustment applied at surface, then subtracted
+! from all vertical fluxes. This is currently done for T,S only
+! 0.0 <= vclose_surf_frac < 1.0
+! 1: fully from surface (i.e. no correction as before); 0: fully from bottom 
+!-----------------------------------------------------------------------
+  vclose_surf_frac=1.0
+
+!-----------------------------------------------------------------------
+! Option to enforce strict mass conservation for each tracer model (only works with itr_met=3,4)
+! At moment the scheme has not accounted for bottom 'leaking' (e.g. in SED), 
+! so iadjust_mass_consv0(5) must =0
+!-----------------------------------------------------------------------
+  iadjust_mass_consv0(1)=0 !T
+  iadjust_mass_consv0(2)=0 !S
+  iadjust_mass_consv0(3)=0 !GEN
+  iadjust_mass_consv0(4)=0 !AGE
+  iadjust_mass_consv0(5)=0 !SED3D (code won't allow non-0 for this module)
+  iadjust_mass_consv0(6)=0 !EcoSim
+  iadjust_mass_consv0(7)=0 !ICM
+  iadjust_mass_consv0(8)=0 !CoSiNE
+  iadjust_mass_consv0(9)=0 !Feco
+  iadjust_mass_consv0(10)=0 !TIMOR
+  iadjust_mass_consv0(11)=0 !FABM
+  iadjust_mass_consv0(12)=0 !DVD (must=0)
+
+! For ICM, impose mass conservation for depths larger than a threshold by considering prism 
+! volume change from step n to n+1. rinflation_icm is the max ratio btw H^{n+1} and H^n allowed.
+  h_massconsv = 2. ![m]
+  rinflation_icm = 1.e-3
+
+/
+
+&SCHOUT
+!-----------------------------------------------------------------------
+! Output section - all optional. Values shown are default unless otherwise stated,
+! and default for most global outputs is off
+!-----------------------------------------------------------------------
+
+!-----------------------------------------------------------------------
+! Main switch to control netcdf. If =0, SCHISM won't output nc files 
+! at all (useful for other programs like ESMF to output)
+!-----------------------------------------------------------------------
+  nc_out = 1
+
+!-----------------------------------------------------------------------
+! UGRID option for _3D_ outputs under scribed IO (out2d*.nc always has meta
+! data info). If iof_ugrid/=0, 3D outputs will also have UGRID metadata (at
+! the expense of file size). 
+!-----------------------------------------------------------------------
+  iof_ugrid = 0
+
+!-----------------------------------------------------------------------
+! Option for hotstart outputs
+!-----------------------------------------------------------------------
+  nhot = 0 !1: output *_hotstart every 'hotout_write' steps
+  nhot_write = 8640 !must be a multiple of ihfskip if nhot=1
+
+!-----------------------------------------------------------------------
+! Station output option. If iout_sta/=0, need output skip (nspool_sta) and
+! a station.in. If ics=2, the cordinates in station.in must be in lon., lat,
+! and z (positive upward; not used for 2D variables). 
+!-----------------------------------------------------------------------
+  iout_sta = 0
+  nspool_sta = 10 !needed if iout_sta/=0; mod(nhot_write,nspool_sta) must=0
+
+!-----------------------------------------------------------------------
+! Global output options
+! The variable names that appear in nc output are shown in {}
+!-----------------------------------------------------------------------
+  iof_hydro(1) = 1 !0: off; 1: on - elev. [m]  {elevation}  2D
+  iof_hydro(2) = 0 !air pressure [Pa]  {airPressure}  2D
+  iof_hydro(3) = 0 !air temperature [C] {airTemperature}  2D
+  iof_hydro(4) = 0 !Specific humidity [-] {specificHumidity}  2D
+  iof_hydro(5) = 0 !Net downward solar (shortwave) radiation after albedo [W/m/m] {solarRadiation}  2D
+  iof_hydro(6) = 0 !sensible flux (positive upward) [W/m/m]  {sensibleHeat}  2D
+  iof_hydro(7) = 0 !latent heat flux (positive upward) [W/m/m] {latentHeat}  2D
+  iof_hydro(8) = 0 !upward longwave radiation (positive upward) [W/m/m] {upwardLongwave}  2D
+  iof_hydro(9) = 0 !downward longwave radiation (positive downward) [W/m/m] {downwardLongwave}  2D
+  iof_hydro(10) = 0 !total flux=-flsu-fllu-(radu-radd) [W/m/m] {totalHeat}  2D
+  iof_hydro(11) = 0 !evaporation rate [kg/m/m/s] {evaporationRate}  2D
+  iof_hydro(12) = 0 !precipitation rate [kg/m/m/s] {precipitationRate}  2D
+  iof_hydro(13) = 0 !Bottom stress vector [kg/m/s^2(Pa)] {bottomStressX,Y}  2D vector
+  iof_hydro(14) = 0 !wind velocity vector [m/s] {windSpeedX,Y}  2D vector
+  iof_hydro(15) = 0 !wind stress vector [m^2/s/s] {windStressX,Y}  2D vector
+  iof_hydro(16) = 0 !depth-averaged vel vector [m/s] {depthAverageVelX,Y}  2D vector
+  iof_hydro(17) = 0 !vertical velocity [m/s] {verticalVelocity}  3D
+  iof_hydro(18) = 0 !water temperature [C] {temperature}  3D
+  iof_hydro(19) = 0 !water salinity [PSU] {salinity}  3D
+  iof_hydro(20) = 0 !water density [kg/m^3] {waterDensity}  3D
+  iof_hydro(21) = 0 !vertical eddy diffusivity [m^2/s] {diffusivity}  3D
+  iof_hydro(22) = 0 !vertical eddy viscosity [m^2/s] {viscosity}  3D
+  iof_hydro(23) = 0 !turbulent kinetic energy {turbulentKineticEner}   3D
+  iof_hydro(24) = 0 !turbulent mixing length [m] {mixingLength}  3D
+!  iof_hydro(25) = 1 !z-coord {zCoordinates} 3D - this flag should be on for visIT etc
+  iof_hydro(26) = 1 !horizontal vel vector [m/s] {horizontalVelX,Y} 3D vector
+  iof_hydro(27) = 0 !horizontal vel vector defined @side [m/s] {horizontalSideVelX,Y} 3D vector
+  iof_hydro(28) = 0 !vertical vel. @elem [m/s] {verticalVelAtElement} 3D
+  iof_hydro(29) = 0 !T @prism centers [C] {temperatureAtElement} 3D
+  iof_hydro(30) = 0 !S @prism centers [PSU] {salinityAtElement} 3D
+  iof_hydro(31) = 0 !Barotropic pressure gradient force vector (m.s-2) @side centers  {pressure_gradient} 2D vector
+
+!-----------------------------------------------------------------------
+! Outputs from optional modules. Only uncomment these if the USE_* is on
+!-----------------------------------------------------------------------
+!-----------------------------------------------------------------------
+! Outputs from DVD (USE_DVD must be on in Makefile)
+!-----------------------------------------------------------------------
+!  iof_dvd(1) = 1 !num mixing for S (PSU^2/s) {DVD_1}   3D elem
+
+!-----------------------------------------------------------------------
+! Outputs from WWM (USE_WWM must be on in Makefile)
+!-----------------------------------------------------------------------
+!  iof_wwm(1)  = 0 !sig. height (m) {sigWaveHeight}   2D
+!  iof_wwm(2)  = 0 !Mean average period (sec) - TM01 {meanWavePeriod}  2D
+!  iof_wwm(3)  = 0 !Zero down crossing period for comparison with buoy (s) - TM02 {zeroDowncrossPeriod}  2D
+!  iof_wwm(4)  = 0 !Average period of wave runup/overtopping - TM10 {TM10}  2D
+!  iof_wwm(5)  = 0 !Mean wave number (1/m) {meanWaveNumber}  2D
+!  iof_wwm(6)  = 0 !Mean wave length (m) {meanWaveLength}  2D
+!  iof_wwm(7)  = 0 !Mean average energy transport direction (degr) - MWD in NDBC? {meanWaveDirection}  2D
+!  iof_wwm(8) = 0 !Mean directional spreading (degr) {meanDirSpreading}  2D
+!  iof_wwm(9) = 0 !Discrete peak period (sec) - Tp {peakPeriod}  2D
+!  iof_wwm(10) = 0 !Continuous peak period based on higher order moments (sec) {continuousPeakPeriod}  2D
+!  iof_wwm(11) = 0 !Peak phase vel. (m/s) {peakPhaseVel}  2D
+!  iof_wwm(12) = 0 !Peak n-factor {peakNFactor}   2D
+!  iof_wwm(13) = 0 !Peak group vel. (m/s) {peakGroupVel}   2D
+!  iof_wwm(14) = 0 !Peak wave number {peakWaveNumber}  2D
+!  iof_wwm(15) = 0 !Peak wave length {peakWaveLength}  2D
+!  iof_wwm(16) = 0 !Peak (dominant) direction (degr) {dominantDirection}  2D
+!  iof_wwm(17) = 0 !Peak directional spreading {peakSpreading}  2D
+!  iof_wwm(18) = 0 !Discrete peak direction (radian?) {discretePeakDirectio}  2D
+!  iof_wwm(19) = 0 !Orbital vel. (m/s) {orbitalVelocity}  2D
+!  iof_wwm(20) = 0 !RMS Orbital vel. (m/s) {rmsOrbitalVelocity}  2D
+!  iof_wwm(21) = 0 !Bottom excursion period (sec?) {bottomExcursionPerio}  2D
+!  iof_wwm(22) = 0 !Bottom wave period (sec) {bottomWavePeriod}  2D
+!  iof_wwm(23) = 0 !Uresell number based on peak period {UresellNumber}  2D
+!  iof_wwm(24) = 0 !Friction velocity (m/s?) {frictionalVelocity}  2D
+!  iof_wwm(25) = 0 !Charnock coefficient {CharnockCoeff}  2D
+!  iof_wwm(26) = 0 !Rougness length {rougnessLength}  2D
+
+!  iof_wwm(27) = 0 !Roller energy dissipation rate (W/m²) @nodes {Drol} 2D
+!  iof_wwm(28) = 0 !Total wave energy dissipation rate by depth-induced breaking (W/m²) @nodes {wave_sbrtot}  2D
+!  iof_wwm(29) = 0 !Total wave energy dissipation rate by bottom friction (W/m²) @nodes {wave_sbftot} 2D
+!  iof_wwm(30) = 0 !Total wave energy dissipation rate by whitecapping (W/m²) @nodes {wave_sdstot} 2D
+!  iof_wwm(31) = 0 !Total wave energy dissipation rate by vegetation (W/m²) @nodes {wave_svegtot} 2D
+!  iof_wwm(32) = 0 !Total wave energy input rate from atmospheric forcing (W/m²) @nodes {wave_sintot} 2D
+!  iof_wwm(33) = 0 !WWM_energy vector {waveEnergyDirX,Y}  2D vector
+
+!  iof_wwm(34) = 0 !Vertical Stokes velocity (m.s-1) @sides and whole levels {stokes_wvel}  3D
+!  iof_wwm(35) = 0 !Wave force vector (m.s-2) computed by wwm @side centers and whole levels {waveForceX,Y}   3D vector
+
+!  iof_wwm(36) = 0 !Horizontal Stokes velocity (m.s-1) @nodes and whole levels {stokes_hvel} 3D vector
+!  iof_wwm(37) = 0 !Roller contribution to horizontal Stokes velocity (m.s-1) @nodes and whole levels {roller_stokes_hvel} 3D vector 
+
+!-----------------------------------------------------------------------
+! Tracer module outputs. In most cases, actual # of outputs depends on # of tracers used
+!-----------------------------------------------------------------------
+! Outputs for user-defined tracer module (USE_GEN)
+!-----------------------------------------------------------------------
+!  iof_gen(1) = 0 !1st tracer {GEN_1}  3D
+!  iof_gen(2) = 0 !2nd tracer {GEN_2}  3D
+
+!-----------------------------------------------------------------------
+! Outputs for (age). Indices from "1" to "ntracer_age/2"; [days]
+!-----------------------------------------------------------------------
+!  iof_age(1) = 0 !{AGE_1}  3D
+!  iof_age(2) = 0 !{AGE_2}  3D
+
+!-----------------------------------------------------------------------
+! Specific outputs in SED3D (USE_SED must be on in Makefile;
+! otherwise these are not needed)
+!-----------------------------------------------------------------------
+!  iof_sed(1) = 0 ! total bed thickness @elem (m) {sedBedThickness}  2D
+!  iof_sed(2) = 0 ! total bed age over all layers @elem (sec) {sedBedAge}  2D
+!  iof_sed(3) = 0 ! Sediment transport roughness length @elem (m) (sedTransportRough) {z0st}  2D
+!  iof_sed(4) = 0 !current-ripples roughness length @elem (m) (sedRoughCurrentRippl) {z0cr}  2D
+!  iof_sed(5) = 0 !sand-waves roughness length (m) @elem (z0sw_elem) {sedRoughSandWave}  2D
+!  iof_sed(6) = 0 !wave-ripples roughness length @elem (m) (z0wr_elem) {sedRoughWaveRipple}  2D
+
+!  iof_sed(7) = 0 !bottom depth _change_ from init. condition (m) {sedDepthChange}  2D
+!  iof_sed(8) = 0 ! Bed median grain size in the active layer (mm) {sedD50}  2D
+!  iof_sed(9) = 0 ! Bottom shear stress (Pa) {sedBedStress}  2D
+!  iof_sed(10) = 0 ! Bottom roughness lenghth (mm) {sedBedRoughness}  2D
+!  iof_sed(11) = 0 ! sediment porosity in the top layer (-) {sedPorocity}  2D
+!  iof_sed(12) = 0 ! erosion flux for suspended transport (kg/m/m/s) {sedErosionalFlux}  2D
+!  iof_sed(13) = 0 ! deposition flux for suspended transport (kg/m/m/s) {sedDepositionalFlux}  2D
+!  iof_sed(14) = 0 ! Bedload transport rate vector due to wave acceleration (kg/m/s) {sedBedloadTransportX,Y}  2D vector
+
+!  Example of using 2 classes
+!  iof_sed(15) = 0 !Bedload transport rate _vector_ (kg.m-1.s-1) for 1st tracer (one output per class) {sedBedload[X,Y]_1}  2D vector
+!  iof_sed(16) = 0 !Bedload transport of 2nd class {sedBedFraction_[X,Y]_2}  2D vector
+!  iof_sed(17) = 0 !Bed fraction 1st tracer (one output per class) [-] {sedBedFraction_1}   2D
+!  iof_sed(18) = 0 !Bed fraction of 2nd class {sedBedFraction_2}   2D
+
+!  iof_sed(19) = 0 !conc. of 1st class (one output need by each class) [g/L] {sedConcentration_1}   3D
+!  iof_sed(20) = 0 !conc. of 2nd class {sedConcentration_2}   3D
+
+!  iof_sed(21) = 0 !total suspended concentration (g/L) {totalSuspendedLoad}  3D
+
+!-----------------------------------------------------------------------
+! EcoSim outputs 
+!-----------------------------------------------------------------------
+!  iof_eco(1) = 0 {ECO_1}  3D
+
+!-----------------------------------------------------------------------
+! ICM outputs 
+!-----------------------------------------------------------------------
+!  !core Module
+!  iof_icm_core(1)  = 1 !PB1
+!  iof_icm_core(2)  = 1 !PB2
+!  iof_icm_core(3)  = 1 !PB3
+!  iof_icm_core(4)  = 1 !RPOC
+!  iof_icm_core(5)  = 1 !LPOC
+!  iof_icm_core(6)  = 1 !DOC
+!  iof_icm_core(7)  = 1 !RPON
+!  iof_icm_core(8)  = 1 !LPON
+!  iof_icm_core(9)  = 1 !DON
+!  iof_icm_core(10) = 1 !NH4
+!  iof_icm_core(11) = 1 !NO3
+!  iof_icm_core(12) = 1 !RPOP
+!  iof_icm_core(13) = 1 !LPOP
+!  iof_icm_core(14) = 1 !DOP
+!  iof_icm_core(15) = 1 !PO4
+!  iof_icm_core(16) = 1 !COD
+!  iof_icm_core(17) = 1 !DOX
+!
+!  !silica module
+!  iof_icm_silica(1) = 1 !SU
+!  iof_icm_silica(2) = 1 !SA
+!
+!  !zooplankton module
+!  iof_icm_zb(1)  = 1 !ZB1
+!  iof_icm_zb(2)  = 1 !ZB2
+!
+!  !pH model
+!  iof_icm_ph(1) = 1 !TIC
+!  iof_icm_ph(2) = 1 !ALK
+!  iof_icm_ph(3) = 1 !CA
+!  iof_icm_ph(4) = 1 !CACO3
+!
+!  !CBP model
+!  iof_icm_cbp(1) = 1 !SRPOC
+!  iof_icm_cbp(2) = 1 !SRPON
+!  iof_icm_cbp(3) = 1 !SRPOP
+!  iof_icm_cbp(4) = 1 !PIP
+!
+!  !SAV model
+!  iof_icm_sav(1) = 1 !stleaf: total leaf biomass @elem [gC/m^2]
+!  iof_icm_sav(2) = 1 !ststem: total stem biomass @elem [gC/m^2]
+!  iof_icm_sav(3) = 1 !stroot: total root biomass @elem [gC/m^2]
+!  iof_icm_sav(4) = 1 !sht:    canopy height @elem [m]
+!
+!  !VEG model
+!  iof_icm_veg(1)  = 1 !vtleaf1: leaf biomass group 1 [gC/m^2]
+!  iof_icm_veg(2)  = 1 !vtleaf2: leaf biomass group 2 [gC/m^2]
+!  iof_icm_veg(3)  = 1 !vtleaf3: leaf biomass group 3 [gC/m^2]
+!  iof_icm_veg(4)  = 1 !vtstem1: stem biomass group 1 [gC/m^2]
+!  iof_icm_veg(5)  = 1 !vtstem2: stem biomass group 2 [gC/m^2]
+!  iof_icm_veg(6)  = 1 !vtstem3: stem biomass group 3 [gC/m^2]
+!  iof_icm_veg(7)  = 1 !vtroot1: root biomass group 1 [gC/m^2]
+!  iof_icm_veg(8)  = 1 !vtroot2: root biomass group 2 [gC/m^2]
+!  iof_icm_veg(9)  = 1 !vtroot3: root biomass group 3 [gC/m^2]
+!  iof_icm_veg(10) = 1 !vht1:    canopy height group 1 [m]
+!  iof_icm_veg(11) = 1 !vht2:    canopy height group 2 [m]
+!  iof_icm_veg(12) = 1 !vht3:    canopy height group 3 [m]
+!
+!  !SFM model
+!  iof_icm_sed(1)   = 1 !bPOC1 (g.m-3)
+!  iof_icm_sed(2)   = 1 !bPOC2 (g.m-3)
+!  iof_icm_sed(3)   = 1 !bPOC3 (g.m-3)
+!  iof_icm_sed(4)   = 1 !bPON1 (g.m-3)
+!  iof_icm_sed(5)   = 1 !bPON2 (g.m-3)
+!  iof_icm_sed(6)   = 1 !bPON3 (g.m-3)
+!  iof_icm_sed(7)   = 1 !bPOP1 (g.m-3)
+!  iof_icm_sed(8)   = 1 !bPOP2 (g.m-3)
+!  iof_icm_sed(9)   = 1 !bPOP3 (g.m-3)
+!  iof_icm_sed(10)  = 1 !bNH4  (g.m-3)
+!  iof_icm_sed(11)  = 1 !bNO3  (g.m-3)
+!  iof_icm_sed(12)  = 1 !bPO4  (g.m-3)
+!  iof_icm_sed(13)  = 1 !bH2S  (g.m-3)
+!  iof_icm_sed(14)  = 1 !bCH4  (g.m-3)
+!  iof_icm_sed(15)  = 1 !bPOS  (g.m-3)
+!  iof_icm_sed(16)  = 1 !bSA   (g.m-3)
+!  iof_icm_sed(17)  = 1 !bstc: surface transfer coeff. (m/day)
+!  iof_icm_sed(18)  = 1 !bSTR: benthic stress      (day)
+!  iof_icm_sed(19)  = 1 !bThp: consective days of hypoxia (day)
+!  iof_icm_sed(20)  = 1 !bTox: consective days of oxic condition after hypoxia (day)
+!  iof_icm_sed(21)  = 1 !SOD   (g.m-2.day-1)
+!  iof_icm_sed(22)  = 1 !JNH4  (g.m-2.day-1)
+!  iof_icm_sed(23)  = 1 !JNO3  (g.m-2.day-1)
+!  iof_icm_sed(24)  = 1 !JPO4  (g.m-2.day-1)
+!  iof_icm_sed(25)  = 1 !JSA   (g.m-2.day-1)
+!  iof_icm_sed(26)  = 1 !JCOD  (g.m-2.day-1)
+!
+!  !Benthic Algae model
+!  iof_icm_ba(1)    = 1 !BA    (g[C].m-2)
+!
+!  !ICM Debug Outputs (need coding, for developers)
+!  iof_icm_dbg(1)   = 1 !2D ICM debug variables
+!  iof_icm_dbg(2)   = 1 !3D ICM debug variables
+           
+!-----------------------------------------------------------------------
+! CoSINE outputs: all 3D
+!-----------------------------------------------------------------------
+!  iof_cos(1)  = 0 !NO3 (uM)
+!  iof_cos(2)  = 0 !SiO4(uM)
+!  iof_cos(3)  = 0 !NH4 (uM)
+!  iof_cos(4)  = 0 !S1  (uM)
+!  iof_cos(5)  = 0 !S2  (uM)
+!  iof_cos(6)  = 0 !Z1  (uM)
+!  iof_cos(7)  = 0 !Z2  (uM)
+!  iof_cos(8)  = 0 !DN  (uM) 
+!  iof_cos(9)  = 0 !DSi (uM) 
+!  iof_cos(10) = 0 !PO4 (uM) 
+!  iof_cos(11) = 0 !DOX (uM) 
+!  iof_cos(12) = 0 !CO2 (uM) 
+!  iof_cos(13) = 0 !ALK (uM) 
+
+!-----------------------------------------------------------------------
+! Fecal indicating bacteria module
+!-----------------------------------------------------------------------
+!  iof_fib(1) = 0 ! {FIB_1}  3D
+
+!-----------------------------------------------------------------------
+! Specific outputs in SED2D (USE_SED2D must be on in Makefile;
+! otherwise these are not needed) - not implemented yet
+!-----------------------------------------------------------------------
+!  iof_sed2d(1)  = 0 !bottom depth _change_ from init. condition (m) {SED2D_depth_change}
+!  iof_sed2d(2)  = 0 !drag coefficient used in transport formulae SED2D_Cd{}
+!  iof_sed2d(3) = 0 !Courant number (b.qtot.dt / h.dx) {SED2D_cflsed}
+!  iof_sed2d(4)    = 0 !Top layer d50 (m) {SED2D_d50}
+!  iof_sed2d(5)   = 0 !total transport rate vector (kg/m/s) {SED2D_total_transport}
+!  iof_sed2d(6)   = 0 !suspended tranport rate vector (kg/m/s) {SED2D_susp_load}
+!  iof_sed2d(7)   = 0 !bedload transport rate vector (kg/m/s) {SED2D_bed_load}
+!  iof_sed2d(8)    = 0 !time averaged total transport rate vector (kg/m/s) {SED2D_average_transport}
+!  iof_sed2d(9)  = 0 !bottom slope vector (m/m); negative uphill {SED2D_bottom_slope}
+!  iof_sed2d(10) = 0 !Total roughness length @elem (m) (z0eq) {z0eq}
+!  iof_sed2d(11) = 0 !current-ripples roughness length @elem (m) (z0cr) {z0cr}
+!  iof_sed2d(12) = 0 !sand-waves roughness length @elem (m) (z0sw) {z0sw}
+!  iof_sed2d(13) = 0 !wave-ripples roughness length @elem (m) (z0wr) {z0wr}
+ 
+!-----------------------------------------------------------------------
+!  marsh flags (USE_MARSH on)
+!-----------------------------------------------------------------------
+!  iof_marsh(1) = 0 ! {marshFlag}  2D elem
+
+!-----------------------------------------------------------------------
+! Ice module outputs (if USE_ICE is on)
+!-----------------------------------------------------------------------
+!  iof_ice(1)= 0 !strain rate @ elem [1/sec] {iceStrainRate}  2D
+!  iof_ice(2)= 0 !ice advective velcoity vector [m/s] {iceVelocityX,Y}  2D vector
+!  iof_ice(3)= 0 !net heat flux to ocean (>0 warm up SST) [W/m/m] {iceNetHeatFlux}  2D
+!  iof_ice(4)= 0 !net fresh water flux to ocean (>0 freshens up SSS) [kg/s/m/m] {iceFreshwaterFlux}  2D
+!  iof_ice(5)= 0 !ice temperature [C] at air-ice interface {iceTopTemperature}  2D
+
+!  iof_ice(6)= 0 !ice volume [m] {iceTracer_1}   2D
+!  iof_ice(7)= 0 !ice concentration [-] {iceTracer_2}  2D
+!  iof_ice(8)= 0 !snow volume [m] {iceTracer_3}  2D
+
+!-----------------------------------------------------------------------
+! Analysis module outputs (USE_ANALYSIS)
+!-----------------------------------------------------------------------
+!  iof_ana(1) = 0 !min time step at each element over all subcycles in horizontal transport solver [s]   {minTransportTimeStep}  2D
+!  iof_ana(2) = 0 !x-component of \nabla air_pres /\rho_0 [m/s/s] {airPressureGradientX}  2D
+!  iof_ana(3) = 0 !y-component of \nabla air_pres /\rho_0 [m/s/s] {airPressureGradientY}  2D
+!  iof_ana(4) = 0  !\alpha*g*\nabla \Psi [m/s/s] (gradient of tidal potential) {tidePotentialGradX}  2D
+!  iof_ana(5) = 0  !\alpha*g*\nabla \Psi [m/s/s] {tidePotentialGradY}  2D
+!  iof_ana(6) = 0 !\nabla \cdot (\mu \nabla u) [m/s/s] (horizontal momentum mixing) {horzontalViscosityX}  3D side
+!  iof_ana(7) = 0 !\nabla \cdot (\mu \nabla v) [m/s/s] {horzontalViscosityY}   3D side
+!  iof_ana(8) = 0 !-g/rho0* \int_z^\eta dr_dx dz  [m/s/s] (b-clinic gradient) {baroclinicForceX}  3D side
+!  iof_ana(9) = 0 !-g/rho0* \int_z^\eta dr_dy dz  [m/s/s] {baroclinicForceY}  3D side
+!  iof_ana(10) = 0 !d (\nu du/dz)/dz [m/s/s] - no vegetation effects (vertical momentum mixing) {verticalViscosityX}  3D side
+!  iof_ana(11) = 0 !d (\nu dv/dz)/dz [m/s/s] - no vegetation effects {verticalViscosityY}  3D side
+!  iof_ana(12) = 0 !(u \cdot \nabla) u [m/s/s] (momentum advection) {mommentumAdvectionX}  3D side
+!  iof_ana(13) = 0 !(u \cdot \nabla) u [m/s/s] {mommentumAdvectionY}  3D side
+!  iof_ana(14) = 0 !gradient Richardson number [-] {gradientRichardson}   3D
+/
diff --git a/pyschism/param/param.py b/pyschism/param/param.py
index 87873893..e56b3176 100644
--- a/pyschism/param/param.py
+++ b/pyschism/param/param.py
@@ -1,7 +1,8 @@
 from datetime import timedelta
+from typing import Union
 import logging
 import pathlib
-from typing import Union
+import tempfile
 
 import f90nml
 
@@ -10,11 +11,9 @@
 from pyschism.param.core import CORE
 from pyschism.param.opt import OPT
 from pyschism.param.schout import SCHOUT
-
+from pyschism.param import schism_init
 # from pyschism.stations import Stations
 
-# PARAM_TEMPLATE = pathlib.Path(__file__).parent / 'param.nml.template'
-from pyschism.param.schism_init import GitParamTemplate
 
 logger = logging.getLogger(__name__)
 
@@ -48,8 +47,12 @@ def write(self, path, overwrite=False, use_template=False):
         if path.is_file() and not overwrite:
             raise IOError(f"File {path} exists and overwrite=False")
         if use_template:
-            PARAM_TEMPLATE = GitParamTemplate().path if use_template is True else use_template
-            f90nml.patch(PARAM_TEMPLATE, self.to_dict(), path)
+            PARAM_TEMPLATE = pathlib.Path(__file__).parent / 'param.nml' if use_template is True else use_template
+            schism_param_sample = schism_init.read_schism_param_sample_patched(PARAM_TEMPLATE)
+            tmpfile = tempfile.NamedTemporaryFile()
+            with open(tmpfile.name, 'w') as fh:
+                fh.write(schism_param_sample)
+            f90nml.patch(tmpfile.name, self.to_dict(), path)
         else:
             with open(path, "w") as f:
                 f.write(str(self))
diff --git a/pyschism/param/schism_init.py b/pyschism/param/schism_init.py
index fee5ea09..d8da0abb 100644
--- a/pyschism/param/schism_init.py
+++ b/pyschism/param/schism_init.py
@@ -111,7 +111,8 @@ def schout(self):
         return SchoutParser(self.schism_init_f90)
 
 
-def patch_sample_param_to_avoid_f90nml_bug(schism_param_sample):
+def read_schism_param_sample_patched(schism_param_file):
+    schism_param_sample = pathlib.Path(schism_param_file).read_text().decode('utf-8').split('\n')
     pattern = r"^[! _A-Za-z0-9]+\([*\d]*(\d+)[^\d]*\)[ ]*="
     var_collection = []
     for i, line in enumerate(schism_param_sample):
@@ -132,7 +133,7 @@ def patch_sample_param_to_avoid_f90nml_bug(schism_param_sample):
                 else:
                     schism_param_sample.insert(i + 1, f"{var}({varcount}) = 0")
                 break
-
+    return '\n'.join(schism_param_sample)
 
 class Singleton(type):
     def __init__(cls, name, bases, dict):
@@ -149,8 +150,11 @@ class GitParamTemplate(metaclass=Singleton):
     def __init__(self, branch="master"):
         url = f"https://raw.githubusercontent.com/schism-dev/schism/{branch}/sample_inputs/param.nml"
         response = urllib.request.urlopen(url)
-        schism_param_sample = response.read().decode("utf-8").split("\n")
-        patch_sample_param_to_avoid_f90nml_bug(schism_param_sample)
+        tmpfile = tempfile.NamedTemporaryFile()
+        with open(tmpfile.name) as fh:
+            fh.write(response.read().decode("utf-8"))
+        # schism_param_sample = response.read().decode("utf-8").split("\n")
+        schism_param_sample = read_schism_param_sample_patched(tmpfile.name)
         tmpfile1 = tempfile.NamedTemporaryFile()
         with open(tmpfile1.name, "w") as f:
             f.write("\n".join(schism_param_sample))
diff --git a/pyschism/utils/jsonencoder.py b/pyschism/utils/jsonencoder.py
new file mode 100644
index 00000000..b2822b49
--- /dev/null
+++ b/pyschism/utils/jsonencoder.py
@@ -0,0 +1,13 @@
+import json
+import numpy as np
+ 
+#https://java2blog.com/object-of-type-int64-is-not-json-serializable/
+class NpEncoder(json.JSONEncoder):
+    def default(self, obj):
+        if isinstance(obj, np.integer):
+            return int(obj)
+        if isinstance(obj, np.floating):
+            return float(obj)
+        if isinstance(obj, np.ndarray):
+            return obj.tolist()
+        return super(NpEncoder, self).default(obj)


Test patch:
diff --git a/examples/example_4/test_hycom.py b/examples/example_4/test_hycom.py
deleted file mode 100644
index b9aa1949..00000000
--- a/examples/example_4/test_hycom.py
+++ /dev/null
@@ -1,21 +0,0 @@
-import numpy as np
-from datetime import datetime, timedelta
-import logging
-from time import time
-
-from netCDF4 import Dataset
-from matplotlib.transforms import Bbox
-
-from pyschism.mesh import Hgrid
-from pyschism.forcing.hycom.hycom import HotStartInventory
-
-logger = logging.getLogger(__name__)
-
-if __name__ ==  '__main__':
-    logger.info('Starting!')
-    start_date = datetime.now()
-    hgrid = Hgrid.open('../ICOGS3D_dep/hgrid.gr3',crs='EPSG:4326')
-    hycom = HotStartInventory()
-    t0 = time()
-    hycom.fetch_data(hgrid, start_date)
-    print(f'Took {time()-t0} seconds')