wgh1 dev

dot_in/wgh1_t0_xn0/forcing_list.csv

@@ -0,0 +1,4 @@
+atm,atm00
+ocn,ocnN
+riv,riv00
+open,NSW

dot_in/wgh1_t0_xn0/make_dot_in.py

@@ -0,0 +1,249 @@
+"""
+This creates and populates directories for ROMS runs on mox or similar.
+Its main product is the .in file used for a ROMS run.
+
+It is designed to work with the BLANK.in file, replacing things like
+$whatever$ with meaningful values.
+
+To test from ipython on mac:
+run make_dot_in -g wgh0 -t t0 -x xn0 -r backfill -s continuation -d 2021.07.04 -bu 0 -np 64
+
+If you call with -short_roms True it will create a .in that runs for a shorter time or
+writes history files more frequently (exact behavior is in the code below).  This can
+be really useful for debugging.
+
+"""
+
+# NOTE: we limit the imports to modules that exist in python3 on mox
+from pathlib import Path
+import sys
+from datetime import datetime, timedelta
+
+pth = Path(__file__).absolute().parent.parent.parent / 'lo_tools' / 'lo_tools'
+if str(pth) not in sys.path:
+    sys.path.append(str(pth))
+import dot_in_argfun as dfun
+import Lfun
+
+dot_in_dir = Path(__file__).absolute().parent
+
+Ldir = dfun.intro() # this handles all the argument passing
+
+# check that the folder name matches the arguments (cludgey)
+ppth = Path(__file__).absolute().parent.name
+gridname_check, tag_check, ex_name_check = (str(ppth)).split('_')
+if Ldir['gridname'] != gridname_check:
+    print('WARNING: gridname mismatch')
+    sys.exit()
+if Ldir['tag'] != tag_check:
+    print('WARNING: tag mismatch')
+    sys.exit()
+if Ldir['ex_name'] != ex_name_check:
+    print('WARNING: ex_name mismatch')
+    sys.exit()
+
+fdt = datetime.strptime(Ldir['date_string'], Lfun.ds_fmt)
+fdt_yesterday = fdt - timedelta(days=1)
+
+print(' --- making dot_in for ' + Ldir['date_string'])
+# initialize dict to hold values that we will substitute into the dot_in file.
+D = dict()
+
+D['EX_NAME'] = Ldir['ex_name'].upper()
+
+#### USER DEFINED VALUES ####
+
+# which ROMS code to use
+D['roms_name'] = 'LO_roms_source'
+
+multi_core = True # use more than one core
+
+if Ldir['run_type'] == 'backfill':
+    days_to_run = 1.0
+elif Ldir['run_type'] == 'forecast':
+    days_to_run = float(Ldir['forecast_days'])
+
+# time step in seconds (should fit evenly into 3600 sec)
+if Ldir['blow_ups'] == 0:
+    dtsec = 20
+elif Ldir['blow_ups'] == 1:
+    dtsec = 15
+elif Ldir['blow_ups'] == 2:
+    dtsec = 10
+elif Ldir['blow_ups'] == 3:
+    dtsec = 6
+elif Ldir['blow_ups'] == 4:
+    dtsec = 4
+elif Ldir['blow_ups'] == 5:
+    dtsec = 2
+else:
+    print('Unsupported number of blow ups: %d' % (Ldir['blow_ups']))
+
+D['ndtfast'] = 20
+
+his_interval = 3600 # seconds to define and write to history files
+rst_interval = 1 # days between writing to the restart file (e.g. 5)
+
+# Find which forcings to look for (search the csv file in this directory).
+# We use the csv file because driver_roms_mox.py also uses it to copy forcing
+# without extra stuff.
+this_dir = Path(__file__).absolute().parent
+with open(this_dir / 'forcing_list.csv', 'r') as f:
+    for line in f:
+        which_force, force_choice = line.strip().split(',')
+        D[which_force] = force_choice
+        # Note that the last of these will not be a force, but will
+        # instead encode the open boundaries, e.g.: D['open'] = 'NSW'
+
+# populate open boundary choices
+for O in list('NSEW'):
+    if O in list(D['open']):
+        D[O+'fs'] = 'Cha'
+        D[O+'2']  = 'Fla'
+        D[O+'3']  = 'RadNud'
+    else:
+        D[O+'fs'] = 'Clo'
+        D[O+'2']  = 'Clo'
+        D[O+'3']  = 'Clo'
+
+#### END USER DEFINED VALUES ####
+
+# DERIVED VALUES
+if multi_core:
+    if Ldir['np_num'] == 64: # for new mox nodes 2*32=64 2019_02
+        ntilei = '8' # number of tiles in I-direction
+        ntilej = '8' # number of tiles in J-direction
+    elif Ldir['np_num'] == 72:
+        ntilei = '6' # number of tiles in I-direction
+        ntilej = '12' # number of tiles in J-direction
+    elif Ldir['np_num'] == 112:
+        ntilei = '8' # number of tiles in I-direction
+        ntilej = '14' # number of tiles in J-direction
+    elif Ldir['np_num'] == 144:
+        ntilei = '8' # number of tiles in I-direction
+        ntilej = '18' # number of tiles in J-direction
+    elif Ldir['np_num'] == 196:
+        ntilei = '14' # number of tiles in I-direction
+        ntilej = '14' # number of tiles in J-direction
+    elif Ldir['np_num'] == 252:
+        ntilei = '14' # number of tiles in I-direction
+        ntilej = '18' # number of tiles in J-direction
+    elif Ldir['np_num'] == 392:
+        ntilei = '14' # number of tiles in I-direction
+        ntilej = '28' # number of tiles in J-direction
+    elif Ldir['np_num'] == 588:
+        ntilei = '21' # number of tiles in I-direction
+        ntilej = '28' # number of tiles in J-direction
+    elif Ldir['np_num'] == 400: # klone
+        ntilei = '20' # number of tiles in I-direction
+        ntilej = '20' # number of tiles in J-direction
+    elif Ldir['np_num'] == 200: # klone
+        ntilei = '10' # number of tiles in I-direction
+        ntilej = '20' # number of tiles in J-direction
+    elif Ldir['np_num'] == 40: # klone
+        ntilei = '5' # number of tiles in I-direction
+        ntilej = '8' # number of tiles in J-direction
+    elif Ldir['np_num'] == 4: # klone
+        ntilei = '2' # number of tiles in I-direction
+        ntilej = '2' # number of tiles in J-direction
+    else:
+        print('Unsupported number of processors: %d' % (Ldir['np_num']))
+else:
+    ntilei = '1'
+    ntilej = '1'
+D['ntilei'] = ntilei
+D['ntilej'] = ntilej
+
+# a string version of dtsec, for the .in file
+if dtsec == int(dtsec):
+    dt = str(dtsec) + '.0d0'
+else:
+    dt = str(dtsec) + 'd0'
+D['dt'] = dt
+
+if Ldir['short_roms']:
+    print(' --- running short roms')
+    his_interval = 10 * dtsec
+    D['ntimes'] = int(his_interval/dtsec) # run for his_interval
+else:
+    D['ntimes'] = int(days_to_run*86400/dtsec)
+
+D['ninfo'] = int(his_interval/dtsec) # how often to write info to the log file (# of time steps)
+D['nhis'] = int(his_interval/dtsec) # how often to write to the history files
+D['ndefhis'] = D['nhis'] # how often to create new history files
+D['nrst'] = int(rst_interval*86400/dtsec)
+
+# file location stuff
+date_string_yesterday = fdt_yesterday.strftime(Lfun.ds_fmt)
+D['dstart'] = int(Lfun.datetime_to_modtime(fdt) / 86400.)
+
+# Paths to forcing various file locations
+D['grid_dir'] = Ldir['grid']
+force_dir = Ldir['LOo'] / 'forcing' / Ldir['gridname'] / ('f' + Ldir['date_string'])
+D['force_dir'] = force_dir
+D['roms_varinfo_dir'] = Ldir['parent'] / 'LO_roms_source_alt' / 'varinfo'
+
+# get horizontal coordinate info
+with open(Ldir['grid'] / 'XY_COORDINATE_INFO.csv','r') as xyf:
+    for line in xyf:
+        ltup = line.split(',')
+        D[ltup[0]] = int(ltup[1]) - 2
+
+# get vertical coordinate info
+with open(Ldir['grid'] / 'S_COORDINATE_INFO.csv','r') as sf:
+    for line in sf:
+        ltup = line.split(',')
+        if ltup[0] != 'ITEMS':
+            D[ltup[0]] = int(ltup[1])
+
+# the output directory and the one from the day before
+out_dir = Ldir['roms_out'] / Ldir['gtagex'] / ('f' + Ldir['date_string'])
+D['out_dir'] = out_dir
+out_dir_yesterday = Ldir['roms_out'] / Ldir['gtagex'] / ('f' + date_string_yesterday)
+Lfun.make_dir(out_dir, clean=True) # make sure it exists and is empty
+
+if Ldir['start_type'] == 'perfect':
+    nrrec = '-1' # '-1' for a perfect restart
+    ininame = 'ocean_rst.nc' # start from restart file
+    ini_fullname = out_dir_yesterday / ininame
+elif Ldir['start_type'] == 'continuation':
+    nrrec = '0' # '0' for a history or ini file
+    ininame = 'ocean_his_0025.nc' # restart from a history file
+    ini_fullname = out_dir_yesterday / ininame
+elif Ldir['start_type'] == 'new':
+    nrrec = '0' # '0' for a history or ini file
+    ininame = 'ocean_ini.nc' # start from ini file
+    ini_fullname = force_dir / D['ocn'] / ininame
+D['nrrec'] = nrrec
+D['ini_fullname'] = ini_fullname
+
+# END DERIVED VALUES
+
+## create liveocean.in ##########################
+f = open(dot_in_dir / 'BLANK.in','r')
+f2 = open(out_dir / 'liveocean.in','w')
+for line in f:
+    for var in D.keys():
+        if '$'+var+'$' in line:
+            line2 = line.replace('$'+var+'$', str(D[var]))
+            line = line2 # needed because we loop over all "var" for line
+        else:
+            line2 = line
+    f2.write(line2)
+f.close()
+f2.close()
+
+# # create bio_Banas.in #######################
+# f = open(dot_in_dir / 'bio_Fennel_BLANK.in','r')
+# bio_dot_in_name = 'bio_Banas.in'
+# f3 = open(out_dir / bio_dot_in_name,'w')
+# for line in f:
+#     for var in D.keys():
+#         if '$'+var+'$' in line:
+#             line2 = line.replace('$'+var+'$', str(D[var]))
+#             line = line2
+#         else:
+#             line2 = line
+#     f3.write(line2)
+# f.close()
+# f3.close()

forcing/ocnN/get_one_time.py

@@ -55,7 +55,7 @@ def get_bounds(x_big, y_big, x_small, y_small, pad=3):
 
 ds.close()
 
-if args.start_type == 'continuation':
+if args.start_type in ['continuation','perfect']:
     pad = 20
 elif args.start_type == 'new':
     pad = 0
@@ -117,9 +117,13 @@ def get_bounds(x_big, y_big, x_small, y_small, pad=3):
         vv[mm[tag]==1] = vtrim[mtrim[tag]==1][xyT[tag].query(xynew[tag], workers=-1)[1]]
         # note that "workers" has replaced "n_jobs"
         if vn == 'zeta':
+            # change sea level to match what was used in pgrid for z_offset.
+            z_offset = -1
+            vv = vv + z_offset
             # enforce a minimum depth
-            zmask = vv+hh <= 0.3
-            vv[zmask] = -hh[zmask] + 0.3
+            min_depth = 0.3
+            zmask = vv+hh <= min_depth
+            vv[zmask] = -hh[zmask] + min_depth
         data_dict[vn][:, :] = vv
     elif dm == 3:
         for nn in range(N):

forcing/ocnN/make_forcing_main.py

@@ -127,6 +127,17 @@ def print_info(fn):
     # print info about the files to the screen
     for fn in nc_list:
         print_info(out_dir / fn)
+    # check on min depth
+    print('')
+    ds = xr.open_dataset(out_dir / 'ocean_clm.nc', decode_times=False)
+    dsg = xr.open_dataset(grid_fn)
+    zz = ds.zeta[0,:,:].values.squeeze()
+    hh = dsg.h.values
+    print('Minimum Depth = %0.2f m' % (np.nanmin(zz+hh)))
+    print('Min zeta = %0.2f m' % (np.nanmin(zz)))
+    print('Max zeta = %0.2f m' % (np.nanmax(zz)))
+    ds.close()
+    dsg.close()
 result_dict['result'] = 'success'
 for fn in nc_list:
     if (out_dir / fn).is_file():

lo_tools/lo_tools/plotting_functions.py

@@ -146,6 +146,11 @@ def add_map_field(ax, ds, vn, vlims_dict, slev=-1, cmap='rainbow', fac=1,
         v = ds[vn][0, slev,:,:].values
     v_scaled = fac*v
 
+    # account for WET_DRY
+    if (tag=='rho') and ('wetdry_mask_rho' in ds.data_vars):
+        mwd = ds.wetdry_mask_rho[0,:,:].values.squeeze()
+        v_scaled[mwd==0] = np.nan
+
     # SETTING COLOR LIMITS
     # First see if they are already set. If so then we are done.
     vlims = vlims_dict[vn]
@@ -166,7 +171,9 @@ def add_map_field(ax, ds, vn, vlims_dict, slev=-1, cmap='rainbow', fac=1,
         # dicts have essentially global scope, so setting it here sets it everywhere
 
     if do_mask_edges:
-        v_scaled = mask_edges(v_scaled, x, y)  
+        v_scaled = mask_edges(v_scaled, x, y)
+
+
 
     cs = ax.pcolormesh(px, py, v_scaled, vmin=vlims[0], vmax=vlims[1], cmap=cmap, alpha=alpha)
     return cs

pgrid/gfun_user.py

@@ -16,7 +16,7 @@
 import gfun
 
 # This is the name of the grid that you are working on.
-gridname = 'wgh0'
+gridname = 'wgh1'
 
 # default s-coordinate info (could override below)
 s_dict = {'THETA_S': 4, 'THETA_B': 2, 'TCLINE': 10, 'N': 30,
@@ -186,6 +186,41 @@ def make_initial_info(gridname=gridname):
         if dch['use_z_offset']:
             z = z + dch['z_offset']
 
+    elif gridname == 'wgh1':
+        # Willapa Bay and Grays Harbor nest
+        dch = gfun.default_choices()
+        aa = [-124.4,-123.7,46.35,47.1]
+        res = 200 # target resolution (m)
+        Lon_vec, Lat_vec = gfu.simple_grid(aa, res)
+        dch['t_list'] = [dch['t_dir'] / 'nw_pacific' / 'nw_pacific.nc']
+        dch['z_offset'] = -1
+        # The docs for nw_pacific say the vertical datum is "sea level" so to match
+        # this we would use z_offset = 0, but the intention here is to make the z=0
+        # level be higher up, so that we catch more of the intertidal when using
+        # WET_DRY. This must be matched by a similar intervention to zeta in ocnN.
+        dch['nudging_edges'] = ['north', 'south', 'west']
+        dch['nudging_days'] = (0.1, 1.0)
+
+        # by setting a small min_depth were are planning to use
+        # WET_DRY in ROMS, but maintaining positive depth
+        # for all water cells
+        dch['min_depth'] = 0.2 # meters (positive down)
+
+        # Make the rho grid.
+        lon, lat = np.meshgrid(Lon_vec, Lat_vec)
+        # initialize the bathymetry array
+        z = np.nan * lon
+        # add bathymetry automatically from files
+        for t_fn in dch['t_list']:
+            print('\nOPENING BATHY FILE: ' + t_fn.name)
+            tlon_vec, tlat_vec, tz = gfu.load_bathy_nc(t_fn)
+            tlon, tlat = np.meshgrid(tlon_vec, tlat_vec)
+            z_part = zfun.interp2(lon, lat, tlon, tlat, tz)
+            # put good values of z_part in z
+            z[~np.isnan(z_part)] = z_part[~np.isnan(z_part)]
+        if dch['use_z_offset']:
+            z = z + dch['z_offset']
+
     elif gridname == 'ae0':
         # analytical model estuary
         dch = gfun.default_choices()

pgrid/plot_grid.py

@@ -68,7 +68,7 @@
 # bathymetry
 fig = plt.figure()
 ax = fig.add_subplot(111)
-cs = ax.pcolormesh(plon, plat, zm, vmin=-20, vmax=0, cmap='Spectral_r')
+cs = ax.pcolormesh(plon, plat, zm, vmin=-5, vmax=0, cmap='Spectral_r')
 # cs = ax.pcolormesh(plon, plat, zm, vmin=-120, vmax=-100, cmap='Spectral_r')
 fig.colorbar(cs, ax=ax)
 if dch['analytical'] == True: