diff --git a/src/grid/zgr.py b/src/grid/zgr.py index 128fdb86..3f1a51ed 100644 --- a/src/grid/zgr.py +++ b/src/grid/zgr.py @@ -212,7 +212,6 @@ def fill_zgrid_vars(zgr_type, grid, hgr_type, e_dict, missing): grid (dict) : vertical grid data dictionary """ # gdep - t_done = "gdept" not in missing if t_done is False: # Fill in the 3D gdept data from 1D gdept_0 @@ -248,10 +247,25 @@ def fill_zgrid_vars(zgr_type, grid, hgr_type, e_dict, missing): if w_done is False: if "e3w" not in missing: grid["gdepw"] = np.cumsum(grid["e3w"], axis=1) + elif zgr_type == "zco": + # Don't waste time calculating gdepw everywhere if it all the same. + gdepw_0 = calc_gdepw( + grid["gdept_0"][np.newaxis, :, np.newaxis, np.newaxis] + ).squeeze() + grid["gdepw"] = np.tile( + gdepw_0, + ( + grid["mbathy"].shape[0], + grid["mbathy"].shape[1], + grid["mbathy"].shape[2], + 1, + ), + ) + grid["gdepw"] = np.transpose(grid["gdepw"], axes=[0, 3, 1, 2]) else: # Fill in the 3D gdepw data from gdept grid["gdepw"] = calc_gdepw(grid["gdept"]) - missing = sorted(list(set(missing) - set(["gdepw"]))) + missing = sorted(list(set(missing) - set(["gdepw"]))) # Calculate other gdep values gdep = horiz_interp_lev(grid["gdept"], grid["gdepw"], zgr_type, hgr_type) @@ -264,7 +278,6 @@ def fill_zgrid_vars(zgr_type, grid, hgr_type, e_dict, missing): grid[vi] = gdep[vi[4:]] # e3 - e3t_done = "e3t" not in missing if e3t_done is False: grid["e3t"] = vert_calc_e3(grid["gdept"], grid["gdepw"], "e3t") @@ -301,8 +314,10 @@ def calc_gdepw(gdept): dep_out = np.zeros((gdept.shape)) indxt = np.arange(gdept.shape[1]) indxw = indxt[:-1] + 0.5 + for j in range(gdept.shape[2]): for i in range(gdept.shape[3]): + # dep_out[0, 1:, j, i] = np.interp(indxw, indxt, gdept[0, :, j, i]) func = interp.interp1d(indxt, gdept[0, :, j, i], kind="cubic") # top gdepw is zero dep_out[0, 1:, j, i] = func(indxw) diff --git a/src/pybdy/nemo_bdy_chunk.py b/src/pybdy/nemo_bdy_chunk.py index 1821b6dc..d47882bf 100644 --- a/src/pybdy/nemo_bdy_chunk.py +++ b/src/pybdy/nemo_bdy_chunk.py @@ -361,8 +361,8 @@ def chunk_large(ibdy, jbdy, chunk_number): for i in range(len(all_chunk)): i_chk = ibdy[chunk_number == all_chunk[i]] j_chk = jbdy[chunk_number == all_chunk[i]] - i_range = np.max(i_chk) - np.min(i_chk) - j_range = np.max(j_chk) - np.min(j_chk) + i_range = (np.max(i_chk) - np.min(i_chk)) + 1 + j_range = (np.max(j_chk) - np.min(j_chk)) + 1 ratio_box[i] = np.sum(chunk_number == all_chunk[i]) / (i_range * j_range) # Calculate which boundaries are too big diff --git a/src/pybdy/nemo_bdy_setup.py b/src/pybdy/nemo_bdy_setup.py index 1f084858..18776800 100644 --- a/src/pybdy/nemo_bdy_setup.py +++ b/src/pybdy/nemo_bdy_setup.py @@ -86,9 +86,9 @@ def _get_var_name_value(self, line): index = "-1" value = "" if name_prefix == "ln": - if value_str.find("true") is not -1: + if value_str.find("true") != -1: value = True - elif value_str.find("false") is not -1: + elif value_str.find("false") != -1: value = False else: raise ValueError( @@ -242,11 +242,11 @@ def _get_val(vars_dictionary, bool_vars_dictionary, line): value = line[1].strip() if name_prefix == "ln": - if value.find("true") is not -1: + if value.find("true") != -1: if (name in vars_dictionary) is not True: vars_dictionary[name] = True bool_vars_dictionary[name] = True - elif value.find("false") is not -1: + elif value.find("false") != -1: if (name in vars_dictionary) is not True: vars_dictionary[name] = False bool_vars_dictionary[name] = False diff --git a/src/pybdy/profiler.py b/src/pybdy/profiler.py index 8225bd85..debd4572 100644 --- a/src/pybdy/profiler.py +++ b/src/pybdy/profiler.py @@ -90,7 +90,6 @@ def process_bdy(setup_filepath=0, mask_gui=False): """ # Start Logger - logger.info("Start NRCT Logging: " + time.asctime()) logger.info("============================================") @@ -100,10 +99,11 @@ def process_bdy(setup_filepath=0, mask_gui=False): Setup = setup.Setup(setup_filepath) # default settings file settings = Setup.settings - logger.info("Reading grid completed") + logger.info("Reading setup completed") bdy_msk = _get_mask(Setup, mask_gui) DstCoord.bdy_msk = bdy_msk == 1 + logger.info("Reading mask completed") DstCoord.hgr = hgr.H_Grid(settings["dst_hgr"], settings["nme_map"], logger, dst=1) DstCoord.zgr = zgr.Z_Grid( @@ -114,8 +114,7 @@ def process_bdy(setup_filepath=0, mask_gui=False): logger, dst=1, ) - - logger.info("Reading mask completed") + logger.info("Reading dst grid completed") bdy_ind = {} # define a dictionary to hold the grid information @@ -136,38 +135,11 @@ def process_bdy(setup_filepath=0, mask_gui=False): logger.info("File: coordinates.bdy.nc generated and populated") # Idenitify number of boundary points - nbdy = {} for grd in ["t", "u", "v"]: nbdy[grd] = len(bdy_ind[grd].bdy_i[:, 0]) - # Gather grid information - - logger.info("Gathering grid information") - SourceCoord.hgr = hgr.H_Grid( - settings["src_hgr"], settings["nme_map"], logger, dst=0 - ) - SourceCoord.zgr = zgr.Z_Grid( - settings["src_zgr"], - settings["nme_map"], - SourceCoord.hgr.grid_type, - SourceCoord.hgr.grid, - logger, - dst=0, - ) - - # Fill horizontal grid information - - try: # if they are masked array convert them to normal arrays - SourceCoord.hgr.grid["glamt"] = SourceCoord.hgr.grid["glamt"].filled() # lon - except Exception: - logger.debug("Not a masked array.") - try: - SourceCoord.hgr.grid["gphit"] = SourceCoord.hgr.grid["gphit"].filled() # lat - except Exception: - logger.debug("Not a masked array.") - # Assign horizontal grid data DstCoord.lonlat = {"t": {}, "u": {}, "v": {}} @@ -196,6 +168,32 @@ def process_bdy(setup_filepath=0, mask_gui=False): logger.info("BDY lons/lats identified from %s", settings["dst_hgr"]) + # Gather grid information + + logger.info("Gathering src grid information") + SourceCoord.hgr = hgr.H_Grid( + settings["src_hgr"], settings["nme_map"], logger, dst=0 + ) + SourceCoord.zgr = zgr.Z_Grid( + settings["src_zgr"], + settings["nme_map"], + SourceCoord.hgr.grid_type, + SourceCoord.hgr.grid, + logger, + dst=0, + ) + + # Fill horizontal grid information + + try: # if they are masked array convert them to normal arrays + SourceCoord.hgr.grid["glamt"] = SourceCoord.hgr.grid["glamt"].filled() # lon + except Exception: + logger.debug("Not a masked array.") + try: + SourceCoord.hgr.grid["gphit"] = SourceCoord.hgr.grid["gphit"].filled() # lat + except Exception: + logger.debug("Not a masked array.") + # Define z at t/u/v points DstCoord.depths = {"t": {}, "u": {}, "v": {}}