diff --git a/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/AVISO/README b/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/AVISO/README index 45f855ba6..7565b1c8e 100644 --- a/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/AVISO/README +++ b/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/AVISO/README @@ -1,22 +1,32 @@ -1. Download adt data - Install the tool "copernicusmarine" following https://help.marine.copernicus.eu/en/articles/7970514-copernicus-marine-toolbox-installation. +Install the tool "copernicusmarine" following https://help.marine.copernicus.eu/en/articles/7970514-copernicus-marine-toolbox-installation. - Sample download command: - copernicusmarine subset --dataset-id cmems_obs-sl_glo_phy-ssh_nrt_allsat-l4-duacs-0.25deg_P1D -v adt -v sla -v err_sla -t 2024-04-01 -T 2024-04-10 -x -66.5 -X -52. -y 6.75 -Y 53.25 -f adt_test.nc --force-download --username xxx --password XXXXXX - "cmems_obs-sl_glo_phy-ssh_nrt_allsat-l4-duacs-0.25deg_P1D" covers period from 31 Dec 2021 to now. - Change dataset-id to "c3s_obs-sl_glo_phy-ssh_my_twosat-l4-duacs-0.25deg_P1D" to download period from "31 Dec 1992 to 6 Jun 2023" +Optional: + - Save login info using "copernicusmarine login"; otherwise you need to provide this info for every download. + - Sample script using python api: - see download_adt.py - - Sample checking the availability of a variable: + - Sample checking the availability of a variable: copernicusmarine describe --include-datasets -c adt | less +Download: + + - CLI: copernicusmarine subset --dataset-id cmems_obs-sl_glo_phy-ssh_nrt_allsat-l4-duacs-0.25deg_P1D -v adt -v sla -v err_sla -t 2024-04-01 -T 2024-04-10 -x -66.5 -X -52. -y 6.75 -Y 53.25 -f adt_test.nc --force-download --username xxx --password XXXXXX + + + - Sample script using python api: download_adt.py, see more instructions at the beginning of the script. + The usage is straightforward. + + For the STOFS-3D-Atlantic setup: + Use "minimum_longitude = -60, maximum_longitude = -53, minimum_latitude = 1, maximum_latitude = 55," + This is a slice only covering the ocean boundary, mainly for elev2D.th.nc + + Use "minimum_longitude = -105, maximum_longitude = -50, minimum_latitude = 1, maximum_latitude = 55," + This covers the entire domain. This is mainly used for setting initial elevation, but you don't need to worry about this for now. + + + +Temporal coverage, change dataset-id: + 31 Dec 2021 to present: "cmems_obs-sl_glo_phy-ssh_nrt_allsat-l4-duacs-0.25deg_P1D" + 31 Dec 1992 to 6 Jun 2023: "c3s_obs-sl_glo_phy-ssh_my_twosat-l4-duacs-0.25deg_P1D" (discontinued on Nov 26, 2024; they probably will have a new id for it) -Replace eta2 in hotstart.nc with ADT: - replace_eta2_aviso.py -Generate elev2D.th.nc with ADT: - gen_bnd_aviso.py (check file name in aviso2schism.py) - modify_elev.py (to apply adjustment in elev) diff --git a/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/Bctides/bctides/bctides.py b/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/Bctides/bctides/bctides.py index 9d8e6075f..44e682e6e 100644 --- a/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/Bctides/bctides/bctides.py +++ b/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/Bctides/bctides/bctides.py @@ -17,49 +17,36 @@ class Bctides: def __init__( self, hgrid, - flags: list = None, + bc_flags: list = None, constituents: Union[str, list] = 'major', database: str = 'tpxo', add_earth_tidal: bool = True, cutoff_depth: float = 50.0, - ethconst: list = None, - vthconst: list = None, - tthconst: list = None, - sthconst: list = None, - tobc: list = None, - sobc: list = None, - relax: list = None, + bc_const: list = None, + bc_relax: list = None, ): """Initialize Bctides ojbect Parameters -------- hgrid: Hgrid object - flags: nested list of bctypes + bc_flags: nested list of bctypes constituents: str or list database: str ('tpxo' or 'fes2014') add_earth_tidal: bool cutoff_depth: float - ethconst: list (constant elevation value for each open boundary) - vthconst: list (constant discharge value for each open boundary) - tthconst: list (constant temperature value for each open boundary) - sthconst: list (constant salinity value for each open boundary) - tobc: list (nuding factor of temperature for each open boundary) - sobc: list (nuding factor of salinity for each open boundary) - realx: list (relaxation constants for inflow and outflow) + bc_const: constant values for each variable at each open boundary + (None if not applicable) + bc_relax: relaxation constants for each variable at each open boundary + (None if not applicable) """ self.hgrid = hgrid self.add_earth_tidal = add_earth_tidal self.cutoff_depth = cutoff_depth self.tides = Tides(constituents=constituents, tidal_database=database) - self.flags = flags - self.ethconst = ethconst - self.vthconst = vthconst - self.tthconst = tthconst - self.sthconst = sthconst - self.tobc = tobc - self.sobc = sobc - self.relax = relax + self.bc_flags = bc_flags + self.bc_const = bc_const + self.bc_relax = bc_relax def __str__(self): @@ -103,9 +90,9 @@ def __str__(self): #get amplitude and phase for each open boundary f.append(f"{len(self.gdf)} !nope") - if len(self.gdf) != len(self.flags): - raise ValueError(f'Number of open boundary {len(self.gdf)} is not consistent with number of given bctypes {len(self.flags)}!') - for ibnd, (boundary, flag) in enumerate(zip(self.gdf.itertuples(), self.flags)): + if len(self.gdf) != len(self.bc_flags): + raise ValueError(f'Number of open boundary {len(self.gdf)} is not consistent with number of given bctypes {len(self.bc_flags)}!') + for ibnd, (boundary, flag) in enumerate(zip(self.gdf.itertuples(), self.bc_flags)): logger.info(f"Processing boundary {ibnd+1}:") #number of nodes and flags line = [ @@ -127,7 +114,7 @@ def __str__(self): logger.warning(f'time history of elevation is read in from elev.th (ASCII)!') elif iettype == 2: logger.info("You are choosing type 2 for elevation, value is {selfethconst[ibnd]} ") - f.append(f"{self.ethconst[ibnd]}") + f.append(f"{self.bc_const[ibnd][0]}") elif iettype == 4: logger.warning('time history of elevation is read in from elev2D.th.nc (netcdf)') elif iettype == 3 or iettype == 5: @@ -152,7 +139,7 @@ def __str__(self): logger.warning("time history of discharge is read in from flux.th (ASCII)!") elif ifltype == 2: logger.info("You are choosing type 2 for velocity, value is {self.vthconst[ibnd]} ") - f.append(f"{self.vthconst[ibnd]}") + f.append(f"{self.bc_const[ibnd][1]}") elif ifltype == 3 or ifltype == 5: if ifltype == 5: logger.warning(f'Combination of 3 and 4, time history of velocity is read in from uv.3D.th.nc!') @@ -165,8 +152,8 @@ def __str__(self): elif ifltype == 4 or -4: logger.warning("time history of velocity (not discharge!) is read in from uv3D.th.nc (netcdf)") if ifltype == -4: - logger.info(f"You are using type -4, relaxation constants for inflow is {self.relax[0]}, for outflow is {self.relax[1]}") - f.append(f"{self.relax[0]} {self.relax[1]} !relaxation constant") + logger.info(f"You are using type -4, relaxation constants for inflow is {self.bc_relax[ibnd][0]}, for outflow is {self.bc_relax[ibnd][1]}") + f.append(f"{self.bc_relax[ibnd][0]} {self.bc_relax[ibnd][1]} !relaxation constant") elif ifltype == -1: raise NotImplementedError(f"Velocity type {ifltype} not implemented yet!") #logger.info(f"Flather type radiation b.c. (iettype must be 0 in this case)!") @@ -179,25 +166,25 @@ def __str__(self): #temperature logger.info(f"Temperature type: {itetype}") if itetype == 0: - logger.warning("Temperature is not sepcified, not input needed!") + logger.warning("Temperature is not sepcified, no input needed!") elif itetype == 1: logger.warning("time history of temperature will be read in from TEM_1.th!") - logger.info(f"Nudging factor for T at boundary {ibnd+1} is {self.tobc[ibnd]}") - f.append(f"{self.tobc[ibnd]} !nudging factor for T") + logger.info(f"Nudging factor for T at boundary {ibnd+1} is {self.bc_relax[ibnd][2]}") + f.append(f"{self.bc_relax[ibnd][2]} !nudging factor for T") elif itetype == 2: - logger.info("You are choosing type 2 for temperature, value is {self.tthconst[ibnd]} ") - f.append(f"{self.tthconst[ibnd]} !T") + logger.info(f"You are choosing type 2 for temperature, value is {self.bc_const[ibnd]}[2] ") + f.append(f"{self.bc_const[ibnd][2]} !T") - logger.info(f"Nudging factor for T at boundary {ibnd+1} is {self.tobc[ibnd]}") - f.append(f"{self.tobc[ibnd]} !nudging factor for T") + logger.info(f"Nudging factor for T at boundary {ibnd+1} is {self.bc_relax[ibnd][2]}") + f.append(f"{self.bc_relax[ibnd][2]} !nudging factor for T") elif itetype == 3: logger.info("Using initial temperature profile for inflow") - logger.info(f"Nudging factor for T at boundary {ibnd+1} is{self.tobc[ibnd]}") - f.append(f"{self.tobc[ibnd]} !nudging factor for T") + logger.info(f"Nudging factor for T at boundary {ibnd+1} is{self.bc_relax[ibnd][2]}") + f.append(f"{self.bc_relax[ibnd][2]} !nudging factor for T") elif itetype == 4: logger.warning("time history of temperature is read in from TEM_3D.th.nc (netcdf)!") - logger.info(f"Nudging factor for T at boundary {ibnd+1} is{self.tobc[ibnd]}") - f.append(f"{self.tobc[ibnd]} !nudging factor for T") + logger.info(f"Nudging factor for T at boundary {ibnd+1} is{self.bc_relax[ibnd][2]}") + f.append(f"{self.bc_relax[ibnd][2]} !nudging factor for T") else: raise IOError(f"Invalid type {itetype} for salinity!") @@ -207,22 +194,22 @@ def __str__(self): logger.info("Salinity is not sepcified, not input needed!") elif isatype == 1: logger.warning("time history of salinity will be read in from SAL_1.th!") - logger.info(f"Nudging factor for salt at boundary {ibnd+1} is {self.sobc[ibnd]}") - f.append(f"{self.sobc[ibnd]} !nudging factor for S") + logger.info(f"Nudging factor for salt at boundary {ibnd+1} is {self.bc_relax[ibnd][3]}") + f.append(f"{self.bc_relax[ibnd][3]} !nudging factor for S") elif isatype == 2: - logger.info("Yor are choosing type 2 for salinity, value is {self.sthconst[ibnd]} ") - f.append(f"{self.sthconst[ibnd]} !S") + logger.info(f"Yor are choosing type 2 for salinity, value is {self.bc_relax[ibnd][3]} ") + f.append(f"{self.bc_relax[ibnd][3]} !S") - logger.info(f"Nudging factor for salt at boundary {ibnd+1} is {self.sobc[ibnd]}") - f.append(f"{self.sobc[ibnd]} !nudging factor for S") + logger.info(f"Nudging factor for salt at boundary {ibnd+1} is {self.bc_relax[ibnd][3]}") + f.append(f"{self.bc_relax[ibnd][3]} !nudging factor for S") elif isatype == 3: logger.info("Using initial salinity profile for inflow") - logger.info(f"Nudging factor for salt at boundary {ibnd+1} is {self.sobc[ibnd]}") - f.append(f"{self.sobc[ibnd]} !nudging factor for S") + logger.info(f"Nudging factor for salt at boundary {ibnd+1} is {self.bc_relax[ibnd][3]}") + f.append(f"{self.bc_relax[ibnd][3]} !nudging factor for S") elif isatype == 4: logger.warning("time history of salinity is read in from SAL_3D.th.nc (netcdf)!") - logger.info(f"Nudging factor for salt at boundary {ibnd+1} is {self.sobc[ibnd]}") - f.append(f"{self.sobc[ibnd]} !nudging factor for S") + logger.info(f"Nudging factor for salt at boundary {ibnd+1} is {self.bc_relax[ibnd][3]}") + f.append(f"{self.bc_relax[ibnd][3]} !nudging factor for S") else: raise IOError(f"Invalid type {isatype} for salinity!") diff --git a/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/Source_sink/Relocate/relocate_source_feeder.py b/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/Source_sink/Relocate/relocate_source_feeder.py index 50e42478c..5f1e102ad 100755 --- a/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/Source_sink/Relocate/relocate_source_feeder.py +++ b/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/Source_sink/Relocate/relocate_source_feeder.py @@ -182,6 +182,7 @@ def cpp2lonlat(x, y, lon0=0, lat0=0): def relocate_sources2( old_ss_dir=None, outdir=None, + no_feeder=False, feeder_info_file=None, hgrid_fname=None, allow_neglection=True, max_search_radius=1000.0, mandatory_sources_coor=np.empty((0, 4)), region_list=None, @@ -191,6 +192,8 @@ def relocate_sources2( :old_ss_dir: the directory of the original source_sink files, which should contain the sources.json file generated by pyschism + :no_feeder: if true, feeder_bases (which should still be present in a mesh without feeders) + will be used; otherwise, feeder_heads will be used. :outdir: the output directory :feeder_info_file: the feeder channel information generated by make_feeder_channel.py :hgrid_fname: the hgrid file name of the new grid with feeders @@ -276,7 +279,13 @@ def relocate_sources2( # -------------------------------------relocation------------------------------------- # Find matching source point at mandatory_sources_coor and feeders # These are the desired new source locations - new_sources_coor = np.r_[mandatory_sources_coor[:, :2], feeder_heads[:, :2]] + if no_feeder: + # If feeders are not implemented in the mesh, use "feeder_bases" instead of "feeder_heads", + # because "feeder_bases" are the locations where the feeders would be placed. + # The relocation still works in the sense that only resolved rivers receive new sources. + new_sources_coor = np.r_[mandatory_sources_coor[:, :2], feeder_bases[:, :2]] + else: + new_sources_coor = np.r_[mandatory_sources_coor[:, :2], feeder_heads[:, :2]] # These are the locations used to search for the closest old source # ; in other words, to link the new source to the old source. # For the mandatory sources, these may be different from the new source locations @@ -426,7 +435,7 @@ def relocate_sources( which is useful for the operational forecast to save time, in this case no need to provide the arguments below; otherwise, the relocation map will be generated. - : no_feeder: if true, feeder_bases (which should still be present in a mesh without feeders) + :no_feeder: if true, feeder_bases (which should still be present in a mesh without feeders) will be used; otherwise, feeder_heads will be used. :allow_neglection: Allow neglecting channels that cannot be matched to a new source location. This is useful when the relocation is restricted to a specific region. diff --git a/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/stofs3d_atl_config.py b/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/stofs3d_atl_config.py index efddfdf92..cbeb4e4f2 100644 --- a/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/stofs3d_atl_config.py +++ b/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/stofs3d_atl_config.py @@ -21,11 +21,14 @@ def __init__( self, startdate=datetime(2017, 12, 1), # start date of the model rnday=60, # number of days to run the model + elev2d_uniform_shift=0.0, # add a uniform shift to elev2D nudging_zone_width=1.5, # in degrees nudging_day=1.0, # in days shapiro_zone_width=2.5, # in degrees shapiro_tilt=2.0, # more abrupt transition in the shapiro zone - bctides_flags=None, # a list of lists, each sublist is a set of flags for an open boundary + bc_flags=None, # a list of lists, each sublist is a set of flags for an open boundary + bc_const=None, # a list of lists, each sublist is a set of constants for Eta, Vel, T, S at teach open boundary + bc_relax=None, # a list of lists, each sublist is a set of relaxations for Eta, Vel, T, S at each open boundary nwm_cache_folder=None, relocate_source=True, feeder_info_file=None, # file containing feeder info, @@ -37,6 +40,7 @@ def __init__( self.startdate = startdate self.rnday = rnday + self.elev2d_uniform_shift = elev2d_uniform_shift self.nudging_zone_width = nudging_zone_width self.nudging_day = nudging_day self.shapiro_zone_width = shapiro_zone_width @@ -45,8 +49,9 @@ def __init__( self.nwm_cache_folder = nwm_cache_folder self.feeder_info_file = feeder_info_file self.mandatory_sources_coor = mandatory_sources_coor - if bctides_flags is None: - self.bctides_flags = [ + + if bc_flags is None: + self.bc_flags = [ [3, 3, 0, 0], # tides for elev and vel [3, 3, 0, 0], # tides for elev and vel [3, 3, 0, 0], # tides for elev and vel @@ -54,7 +59,29 @@ def __init__( [3, 3, 0, 0], # tides for elev and vel ] else: - self.bctides_flags = bctides_flags + self.bc_flags = bc_flags + + if bc_const is None: + self.bc_const = [ + [None, None, None, None], + [None, None, None, None], + [None, None, None, None], + [None, None, None, None], + [None, None, None, None], + ] + else: + self.bc_const = bc_const + + if bc_relax is None: + self.bc_relax = [ + [None, None, None, None], + [None, None, None, None], + [None, None, None, None], + [None, None, None, None], + [None, None, None, None], + ] + else: + self.bc_relax = bc_relax if gr3_values is None: self.gr3_values = { # uniform gr3 values @@ -83,7 +110,7 @@ def v6(cls): '/sciclone/schism10/feiye/STOFS3D-v5/Inputs/v14/Parallel/' 'SMS_proj/feeder/feeder.pkl'), mandatory_sources_coor=rsf.v19p2_mandatory_sources_coor, - bctides_flags=[[5, 5, 4, 4]], + bc_flags=[[5, 5, 4, 4]], tvd_regions=[ 'tvd0_1.reg', 'tvd0_2.reg', 'tvd0_3.reg', 'tvd0_4.reg', 'tvd0_5.reg', 'tvd0_6.reg', 'tvd0_7.reg', @@ -95,6 +122,7 @@ def v6(cls): def v7(cls): '''Factory method to create a configuration for STOFS3D-v7''' return cls( + elev2d_uniform_shift=-0.42, # add a uniform shift to elev2D nudging_zone_width=7.3, # default nudging zone shapiro_zone_width=11.5, # default shapiro zone shapiro_tilt=3.5, # default abrupt transition in the shapiro zone @@ -103,11 +131,21 @@ def v7(cls): 'feeder_heads_bases.xy'), mandatory_sources_coor=rsf.v19p2_mandatory_sources_coor, nwm_cache_folder=Path('/sciclone/schism10/whuang07/schism20/NWM_v2.1/'), - bctides_flags=[ + bc_flags=[ [5, 5, 4, 4], # Atlantic Ocean [5, 5, 4, 4], # Gulf of St. Lawrence [0, 1, 2, 2], # St. Lawrence River ], + bc_const=[ + [None, None, None, None], # Atlantic Ocean + [None, None, None, None], # Gulf of St. Lawrence + [None, None, 10.0, 0.0], # St. Lawrence River + ], + bc_relax=[ # relaxation timescale for each boundary variable + [None, None, 0.5, 0.5], # Atlantic Ocean + [None, None, 0.5, 0.5], # Gulf of St. Lawrence + [None, None, 0.01, 1.0], # St. Lawrence River + ], tvd_regions=[ 'tvd0_1.reg', 'tvd0_2.reg', 'tvd0_3.reg', 'tvd0_4.reg', 'tvd0_5.reg', 'tvd0_6.reg', 'tvd0_7.reg', @@ -117,8 +155,9 @@ def v7(cls): @classmethod def v7_hercules_test(cls): - '''Factory method to create a configuration for STOFS3D-v7''' + '''Factory method to create a configuration for STOFS3D-v7 2D setup''' return cls( + elev2d_uniform_shift=-0.42, # add a uniform shift to elev2D nudging_zone_width=7.3, # default nudging zone shapiro_zone_width=11.5, # default shapiro zone shapiro_tilt=3.5, # default abrupt transition in the shapiro zone @@ -128,7 +167,7 @@ def v7_hercules_test(cls): 'Feeder_channels/feeder_heads_bases.xy'), mandatory_sources_coor=rsf.v19p2_mandatory_sources_coor, nwm_cache_folder=None, - bctides_flags=[ + bc_flags=[ [3, 3, 0, 0], # Atlantic Ocean [3, 3, 0, 0], # Gulf of St. Lawrence [0, 1, 0, 0], # St. Lawrence River @@ -150,13 +189,14 @@ def v8_louisianna(cls): feeder_info_file='', relocate_source=False, nwm_cache_folder=Path('/sciclone/schism10/whuang07/schism20/NWM_v2.1/'), - bctides_flags=[[5, 5, 4, 4]] + bc_flags=[[5, 5, 4, 4]] ) @classmethod def v8(cls): - '''Factory method to create a configuration for STOFS3D-v8''' + '''Factory method to create a configuration for STOFS3D-v8 3D setup''' return cls( + elev2d_uniform_shift=-0.42, # add a uniform shift to elev2D nudging_zone_width=7.3, # default nudging zone shapiro_zone_width=11.5, # default shapiro zone shapiro_tilt=3.5, # default abrupt transition in the shapiro zone @@ -165,10 +205,55 @@ def v8(cls): 'feeder_heads_bases.xy'), mandatory_sources_coor=rsf.v23p3_mandatory_sources_coor, nwm_cache_folder=None, - bctides_flags=[ - [3, 3, 0, 0], # Atlantic Ocean - [3, 3, 0, 0], # Gulf of St. Lawrence - [0, 1, 0, 0], # St. Lawrence River + bc_flags=[ + [5, 5, 4, 4], # Atlantic Ocean + [5, 5, 4, 4], # Gulf of St. Lawrence + [0, 1, 2, 2], # St. Lawrence River + ], + bc_const=[ + [None, None, None, None], # Atlantic Ocean + [None, None, None, None], # Gulf of St. Lawrence + [None, None, 10.0, 0.0], # St. Lawrence River + ], + bc_relax=[ # relaxation timescale for each boundary variable + [None, None, 0.5, 0.5], # Atlantic Ocean + [None, None, 0.5, 0.5], # Gulf of St. Lawrence + [None, None, 0.01, 1.0], # St. Lawrence River + ], + tvd_regions=[ + 'tvd0_1.reg', 'tvd0_2.reg', 'tvd0_3.reg', 'tvd0_4.reg', + 'tvd0_5.reg', 'tvd0_6.reg', 'tvd0_7.reg', + 'upwind_east_Carribbean.rgn', 'upwind_west_Carribbean.rgn' + ] + ) + + @classmethod + def v8_Hercules(cls): + '''Factory method to create a configuration for STOFS3D-v8 3D setup''' + return cls( + elev2d_uniform_shift=-0.42, # add a uniform shift to elev2D + nudging_zone_width=7.3, # default nudging zone + shapiro_zone_width=11.5, # default shapiro zone + shapiro_tilt=3.5, # default abrupt transition in the shapiro zone + feeder_info_file=( + '/work/noaa/nosofs/feiye/STOFS-3D-Atl-Example-Setup/DATA/' + 'Feeder_channels/feeder_heads_bases.xy'), + mandatory_sources_coor=rsf.v19p2_mandatory_sources_coor, + nwm_cache_folder=None, + bc_flags=[ + [5, 5, 4, 4], # Atlantic Ocean + [5, 5, 4, 4], # Gulf of St. Lawrence + [0, 1, 2, 2], # St. Lawrence River + ], + bc_const=[ + [None, None, None, None], # Atlantic Ocean + [None, None, None, None], # Gulf of St. Lawrence + [None, None, 10.0, 0.0], # St. Lawrence River + ], + bc_relax=[ # relaxation timescale for each boundary variable + [None, None, 0.5, 0.5], # Atlantic Ocean + [None, None, 0.5, 0.5], # Gulf of St. Lawrence + [None, None, 0.01, 1.0], # St. Lawrence River ], tvd_regions=[ 'tvd0_1.reg', 'tvd0_2.reg', 'tvd0_3.reg', 'tvd0_4.reg', diff --git a/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/stofs3d_atl_driver.py b/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/stofs3d_atl_driver.py index 2b024e3c3..e8298a869 100755 --- a/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/stofs3d_atl_driver.py +++ b/src/Utility/Pre-Processing/STOFS-3D-Atl-shadow-VIMS/Pre_processing/stofs3d_atl_driver.py @@ -4,6 +4,12 @@ Simpler tasks such as generating uniform *.gr3 are included in this script. More complex tasks such as generating source_sink are imported as modules from subfolders. +Some scripts are written in Fortran and C++ for speed. +You may need to compile them before running this driver script. +See compile instructions at the beginning of each script. + +- Vgrid/: gen_vqs.f90, change_vgrid.f90 + Usage: Under the current schism git directory, @@ -16,7 +22,9 @@ Alternatively, you can copy the current directory (Pre_processing/) to your preferred location and run the script there. -Default settings for the STOFS-3D-ATL model are set in the ConfigStofs3dAtlantic class. +Default settings for the STOFS-3D-ATL model are set in the ConfigStofs3dAtlantic class +in the stofs3d_atl_config.py file. + Read the docstring of the main function for how to use a configuration preset. Also, set paths at the beginning of the main function for each new setup. @@ -33,14 +41,16 @@ import shutil from datetime import datetime import copy +import logging import numpy as np from scipy.spatial import cKDTree import shapefile # from pyshp +import netCDF4 # self-defined modules import pylib -from pylib import inside_polygon, read_schism_reg +from pylib import inside_polygon, read_schism_reg, read_schism_vgrid from pylib_experimental.schism_file import source_sink if 'sciclone' in socket.gethostname(): from pylib_experimental.schism_file import cread_schism_hgrid as schism_read @@ -49,6 +59,8 @@ from pylib import schism_grid as schism_read print('Using python function to read hgrid') from pyschism.mesh import Hgrid, Vgrid +from pyschism.forcing.hycom.hycom2schism import OpenBoundaryInventory as OpenBoundaryInventoryHYCOM +from AVISO.aviso2schism import OpenBoundaryInventory as OpenBoundaryInventoryAVISO from pyschism.forcing.hycom.hycom2schism import Nudge # Import from the sub folders. These are not from installed packages. @@ -76,10 +88,11 @@ # --------------------------------------------------------------------- # Utility functions # --------------------------------------------------------------------- -def mkcd_new_dir(path): +def mkcd_new_dir(path, remove=True): '''Make a new directory and change to it''' - remove_folder(path) - os.makedirs(path) + if remove: + remove_folder(path) + os.makedirs(path, exist_ok=True) os.chdir(path) @@ -202,6 +215,105 @@ def gen_nudge_coef(hgrid: pylib.schism_grid, rlmax=1.5, rnu_day=0.25, open_bnd_l return nudge_coeff +def gen_3dbc(hgrid_fname, vgrid_fname, outdir, start_date, rnday): + """ + Generate 3D boundary conditions based on the HYCOM data + + Sometimes the server can be busy, and the connection may be lost. + """ + # check input files + if not os.path.exists(hgrid_fname): + raise FileNotFoundError(f'{hgrid_fname} not found') + if not os.path.exists(vgrid_fname): + raise FileNotFoundError(f'{vgrid_fname} not found') + # enable pyschism's logging + logging.basicConfig( + format="[%(asctime)s] %(name)s %(levelname)s: %(message)s", + force=True, + ) + logger = logging.getLogger('pyschism') + logger.setLevel(logging.INFO) + + max_attempts = 10 + restart = False # first attempt doesn't need "restart" + + hgrid = Hgrid.open(hgrid_fname, crs='epsg:4326') + bnd = OpenBoundaryInventoryHYCOM(hgrid, vgrid_fname) + + nattempts = max_attempts + while nattempts > 0: + try: + if nattempts < max_attempts: + restart = True # continue from the last attempt + bnd.fetch_data( + outdir, start_date, rnday, + elev2D=True, TS=True, UV=True, + ocean_bnd_ids=[0, 1], restart=restart + ) + break # Exit the loop if fetch_data is successful + except Exception as e: + nattempts -= 1 + print(f"Attempt failed: {e}. Retrying... ({nattempts} attempts left)") + + if nattempts == 0: + raise Exception(f"Failed to fetch 3D boundary conditions after {max_attempts} attempts") + + +def gen_elev2d(hgrid_fname, outdir, start_date, rnday, uniform_shift=0.0): + ''' + Generate 2D elevation boundary conditions based on AVISO (CMEMS) data + + Needs to pre-download data from AVISO website, because + some clusters like Hercules and SciClone do not allow direct downloading + via copernicusmarine. + + See instructions in AVISO/README and sample download script in AVISO/download_aviso*.py + + Save the downloaded data in the current diretory as aviso.nc + + A uniform shift can be applied to the elevation data, e.g., -0.42 m in STOFS-3D v7 + ''' + if not os.path.exists(hgrid_fname): + raise FileNotFoundError(f'{hgrid_fname} not found') + # enable pyschism's logging + logging.basicConfig( + format="[%(asctime)s] %(name)s %(levelname)s: %(message)s", + force=True, + ) + logger = logging.getLogger('pyschism') + logger.setLevel(logging.INFO) + + while True: # search for aviso.nc until it is provided + if not os.path.exists('aviso.nc'): + print('\n' + '-'*50) + print(f'aviso.nc not found, please download the data to {outdir}') + print(f'See instructions in {script_path}/AVISO/README and ' + f'{script_path}/AVISO/download_aviso*.py') + input('Press Enter to continue after downloading aviso.nc') + print('\n' + '-'*50) + else: + break + + hgrid = Hgrid.open(hgrid_fname, crs='epsg:4326') + bnd = OpenBoundaryInventoryAVISO(hgrid) + bnd.fetch_data(outdir, start_date, rnday, ocean_bnd_ids=[0, 1], elev2D=True) + + # modify the elevation data based on the uniform shift + if np.abs(uniform_shift) > 1e-6: + print(f'shifting {uniform_shift}') + + with netCDF4.Dataset(f'{outdir}/elev2D.th.nc', mode="r+") as nc_file: + time_series = nc_file.variables["time_series"] + print("Original values:", time_series[:10]) + time_series[:] = time_series[:] + uniform_shift + print("Updated values:", time_series[:10]) + + # os.system(f'mv {outdir}/elev2D.th.nc {outdir}/elev2D.th.nc.aviso') + # ds = xr.open_dataset(f'{outdir}/elev2D.th.nc.aviso', engine="netcdf4") + # ds['time_series'] += uniform_shift + # ds.to_netcdf(f'{outdir}/elev2D.th.nc', engine="netcdf4") + + def gen_nudge_stofs(hgrid_fname, vgrid_fname, outdir, start_date, rnday): ''' Generate nudge coefficient, @@ -372,26 +484,28 @@ def main(): # -----------------input--------------------- # hgrid generated by SMS, pre-processed, and converted to *.gr3 - hgrid_path = ('/sciclone/schism10/Hgrid_projects/STOFS3D-v8/v23.1/hgrid_with_bnd.gr3') + # hgrid_path = ('/sciclone/schism10/feiye/STOFS3D-v8/I11/hgrid.gr3') + hgrid_path = ('/work/noaa/nosofs/feiye/Runs/R16x/hgrid.gr3') # get a configuration preset and make changes if necessary # alternatively, you can set the parameters directly on an # new instance of ConfigStofs3dAtlantic config = ConfigStofs3dAtlantic.v8() - config.rnday = 36 + config.rnday = 3 config.startdate = datetime(2024, 3, 5) config.nwm_cache_folder = Path( '/sciclone/schism10/feiye/STOFS3D-v8/I09/Source_sink/original_source_sink/20240305/') # define the project dir, where the run folders are located - project_dir = '/sciclone/schism10/feiye/STOFS3D-v8/' + # project_dir = '/sciclone/schism10/feiye/STOFS3D-v8/' + project_dir = '/work/noaa/nosofs/feiye/Runs/' # run ID. e.g, 00, 01a, 10b, etc.; the run folders are named using this ID as follows: # I{runid}: input directory for holding model input files and scripts; # R{runid}: run directory, where the run will be submitted to queue; # O{runid}: output directory for holding raw outputs and post-processing. # under project_dir - runid = '11a' + runid = '16y' # swithes to generate different input files input_files = { @@ -404,8 +518,9 @@ def main(): 'elev_ic': False, 'source_sink': False, 'hotstart.nc': False, - '*D.th.nc': False, - '*nu.nc': True, + '3D.th.nc': False, + 'elev2D.th.nc': True, + '*nu.nc': False, '*.prop': False, } # -----------------end input--------------------- @@ -431,8 +546,10 @@ def main(): mkcd_new_dir(f'{model_input_path}/{sub_dir}') Bctides( hgrid=hgrid_pyschism, # needs hgrid in pyschism's Hgrid class - flags=config.bctides_flags, - database='fes2014' + bc_flags=config.bc_flags, + bc_const=config.bc_const, + bc_relax=config.bc_relax, + database='fes2014', ).write( output_directory=f'{model_input_path}/{sub_dir}', start_date=config.startdate, @@ -452,28 +569,17 @@ def main(): # call a fortran program to generate vgrid.in print(f'compile the fortran program {script_path}/Vgrid/gen_vqs if necessary') - # fortran_script_path = '/path/to/your/fortran_script.f90' - # fortran_command = ['gfortran', fortran_script_path, '-o', 'fortran_script'] fortran_process = subprocess.Popen( f'{script_path}/Vgrid/gen_vqs', stdin=subprocess.PIPE - ) # , stdout=subprocess.PIPE, stderr=subprocess.PIPE) + ) # the command line argument 0 means no outputs of a sample vgrid along a given transect fortran_process.communicate(input=str(0).encode()) print(f'{DRIVER_PRINT_PREFIX}converting the format of vgrid.in ...') - # rename vgrid.in to vgrid.in.old - try_remove(f'{model_input_path}/{sub_dir}/vgrid.in.old') + vg = read_schism_vgrid(f'{model_input_path}/{sub_dir}/vgrid.in') os.rename('vgrid.in', 'vgrid.in.old') - - # convert the format of the vgrid.in file using a fortran script - subprocess.run( - [f'{script_path}/Vgrid/change_vgrid'], check=True - ) # , stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL) - - # rename vgrid.in.new to vgrid.in - try_remove(f'{model_input_path}/{sub_dir}/vgrid.in') - os.rename('vgrid.in.new', 'vgrid.in') + vg.save(f'{model_input_path}/{sub_dir}/vgrid.in') os.chdir(model_input_path) os.system(f'ln -sf {sub_dir}/vgrid.in .') @@ -673,14 +779,35 @@ def main(): # generate hotstart.nc raise NotImplementedError('hotstart.nc is not implemented yet') - # -----------------*D.th.nc--------------------- - if input_files['*D.th.nc']: - sub_dir = 'D.th' - print(f'{DRIVER_PRINT_PREFIX}Generating *D.th.nc files ...') + # -----------------3D.th.nc--------------------- + if input_files['3D.th.nc']: + sub_dir = '3Dth' + print(f'{DRIVER_PRINT_PREFIX}Generating *3D.th.nc files ...') mkcd_new_dir(f'{model_input_path}/{sub_dir}') # generate *D.th.nc - raise NotImplementedError('*D.th.nc is not implemented yet') + gen_3dbc( + hgrid_fname=f'{model_input_path}/hgrid.gr3', + vgrid_fname=f'{model_input_path}/vgrid.in', + outdir=f'{model_input_path}/{sub_dir}', + start_date=config.startdate, rnday=config.rnday + ) + + os.chdir(run_dir) + os.system(f'ln -sf ../I{runid}/{sub_dir}/*3D.th.nc .') + + # -----------------elev2D.th.nc--------------------- + if input_files['elev2D.th.nc']: + sub_dir = 'Elev2Dth' + print(f'{DRIVER_PRINT_PREFIX}Generating elev2D.th.nc ...') + mkcd_new_dir(f'{model_input_path}/{sub_dir}', remove=False) + + gen_elev2d( + hgrid_fname=f'{model_input_path}/hgrid.gr3', + outdir=f'{model_input_path}/{sub_dir}', + start_date=config.startdate, rnday=config.rnday, + uniform_shift=config.elev2d_uniform_shift + ) # -----------------*nu.nc--------------------- if input_files['*nu.nc']: