Add two new tools: one automatically fixes non-advective cells found …

…in a mask file; the other combines two masks. Also modified the topog2mask.py utility to take treat cells with a fraction of ocean lower than a user-provided value as land.

combine_masks.py

@@ -0,0 +1,58 @@
+#!/usr/bin/env python
+
+import sys
+import os
+import numpy as np
+import argparse
+import netCDF4 as nc
+
+def fix_mask(mask1_filename, mask2_filename, grid_filename, lat, output_mask_filename):
+
+    with nc.Dataset(mask1_filename, 'r') as mf1, \
+         nc.Dataset(mask2_filename, 'r') as mf2, \
+         nc.Dataset(grid_filename, 'r') as gr, \
+         nc.Dataset(output_mask_filename, 'w') as mf_out:
+
+        # Sanity checks
+        num_lons = mf1.dimensions['nx'].size
+        num_lats = mf1.dimensions['ny'].size
+        if 2*num_lons != gr.dimensions['nx'].size or 2*num_lats != gr.dimensions['ny'].size:
+            print('Error: dimensions of {} and {} are not compatible'.format(mask1_filename, grid_filename), file=sys.stderr)
+            return 1
+
+        if num_lons != mf2.dimensions['nx'].size or num_lats != mf2.dimensions['ny'].size:
+            print('Error: dimensions of {} and {} are not compatible'.format(mask1_filename, mask2_filename), file=sys.stderr)
+            return 1
+
+        # Get the coordinates of the T points from the supergrid:
+        yt = gr.variables['y'][1::2,1::2]
+
+        # Combine the two masks, taking mask 1 north of 'lat' and mask 2 south of 'lat'.
+        combined = np.where(yt > lat, mf1.variables['mask'][:], mf2.variables['mask'][:])
+
+        # Output mask
+        mf_out.createDimension('nx', num_lons)
+        mf_out.createDimension('ny', num_lats)
+        mask_out = mf_out.createVariable('mask', 'f8', dimensions=('ny', 'nx'))
+        mask_out[:] = combined[:]
+
+def main():
+
+    parser = argparse.ArgumentParser()
+    parser.add_argument('mask1', help='First mask to combine.')
+    parser.add_argument('mask2', help='Second mask to combine.')
+    parser.add_argument('grid', help='Ocean super-grid.')
+    parser.add_argument('output_mask', help='The combined mask.')
+    parser.add_argument('latitude', type=float, help='Take mask1 north of this value and mask2 otherwise.')
+
+    args = parser.parse_args()
+
+    if not os.path.exists(args.mask1) or not os.path.exists(args.mask2) or not os.path.exists(args.grid):
+        print('Error: one of {}, {} or {} was not found'.format(args.topog, args.input_mask, args.grid), file=sys.stderr)
+        parser.print_help()
+        return 1
+
+    return fix_mask(args.mask1, args.mask2, args.grid, args.latitude, args.output_mask)
+
+if __name__ == '__main__':
+    sys.exit(main())

fix_nonadvective_mask.py

@@ -0,0 +1,72 @@
+#!/usr/bin/env python
+
+import sys
+import os
+import numpy as np
+import argparse
+import netCDF4 as nc
+
+def fix_mask(input_mask_filename, output_mask_filename):
+
+    with nc.Dataset(input_mask_filename, 'r') as mf_in, \
+         nc.Dataset(output_mask_filename, 'w') as mf_out:
+
+        mask = mf_in.variables['mask'][:,:]
+        num_lons = mf_in.dimensions['nx'].size
+        num_lats = mf_in.dimensions['ny'].size
+        print('mask dimensions {:11d} {:11d}'.format(num_lons, num_lats))
+
+        # Create mask with halo (includes copied points along LH edge and RH edge and points along tripole seam)
+        mask_halo = np.zeros((num_lats+1,num_lons+2), dtype=int)
+        mask_halo[:-1,1:-1] = mask[:,:]
+
+        mask_halo[:,-1] = mask_halo[:,1]
+        mask_halo[:,0] = mask_halo[:,-2]
+        mask_halo[-1, :] = mask_halo[-2, ::-1]
+
+        # Find non-advective cells
+        done = False
+        n_total = 0
+        print('non-advective cells:')
+        while (not done):
+            sw = (mask_halo[1:-1,  :-2] == 0) | (mask_halo[ :-2,  :-2] == 0) | (mask_halo[ :-2, 1:-1] == 0)
+            se = (mask_halo[ :-2, 1:-1] == 0) | (mask_halo[ :-2, 2:  ] == 0) | (mask_halo[1:-1, 2:  ] == 0)
+            ne = (mask_halo[1:-1, 2:  ] == 0) | (mask_halo[2:  , 2:  ] == 0) | (mask_halo[2:  , 1:-1] == 0)
+            nw = (mask_halo[2:  , 1:-1] == 0) | (mask_halo[2:  ,  :-2] == 0) | (mask_halo[1:-1,  :-2] == 0)
+            advective_cells = sw & se & ne & nw & (mask_halo[1:-1,1:-1] == 1)
+
+            jj,ii = np.nonzero(advective_cells)
+            for ij in list(zip(ii,jj+1)):
+                print('{:11d} {:11d}'.format(ij[0], ij[1]))
+
+            mask_halo[jj+1,ii+1] = 0
+
+            n_new = np.count_nonzero(advective_cells)
+            done = n_new == 0
+            n_total += n_new
+
+        print('Found {} non-advective cells'.format(n_total))
+
+        mf_out.createDimension('nx', num_lons)
+        mf_out.createDimension('ny', num_lats)
+
+        mask_out = mf_out.createVariable('mask', 'f8', dimensions=('ny', 'nx'))
+        mask_out[:,:] = mask_halo[:-1,1:-1]
+
+def main():
+
+    parser = argparse.ArgumentParser()
+    parser.add_argument('input_mask',  help='Mask to fix.')
+    parser.add_argument('output_mask', help='Output mask.')
+
+    args = parser.parse_args()
+
+    if not os.path.exists(args.input_mask):
+        print('Error: {} was not found'.format(args.input_mask), file=sys.stderr)
+        parser.print_help()
+        return 1
+
+    fix_mask(args.input_mask, args.output_mask)
+
+if __name__ == '__main__':
+    sys.exit(main())

topog2mask.py

@@ -8,7 +8,7 @@
 import argparse
 import netCDF4 as nc
 
-def make_mask(topog_filename, mask_filename, ice=False):
+def make_mask(topog_filename, mask_filename, frac, ice=False):
 
     with nc.Dataset(topog_filename, 'r') as tf, \
             nc.Dataset(mask_filename, 'w') as mf:
@@ -26,10 +26,7 @@ def make_mask(topog_filename, mask_filename, ice=False):
         else:
             mask = mf.createVariable('mask', 'f8', dimensions=('ny', 'nx'))
         # CICE and MOM use 0 as masked
-        m = np.zeros_like(mask)
-        m[:] = 1.0
-        m[np.where(tf.variables['depth'][:] <= 0.0)] = 0.0
-        mask[:] = m[:]
+        mask[:] = np.where((tf.variables['frac'][:] < frac) | (tf.variables['depth'][:] <= 0.0), 0, 1)
 
 
 def main():
@@ -38,7 +35,7 @@ def main():
     parser.add_argument('topog', help='The topog file.')
     parser.add_argument('cice_mask',  help='The new CICE mask file.')
     parser.add_argument('mom_mask',  help='The new MOM mask file.')
-
+    parser.add_argument('fraction', type=float, nargs='?', default=0.0, help='Cells with a fraction of ocean lower than this value will be treated as land')
     args = parser.parse_args()
 
     if not os.path.exists(args.topog):
@@ -53,8 +50,14 @@ def main():
         parser.print_help()
         return 1
 
-    make_mask(args.topog, args.cice_mask, ice=True)
-    make_mask(args.topog, args.mom_mask)
+    if args.fraction > 1.0 or args.fraction < 0.0:
+        print('Error: fraction must be between 0 and 1, but it is {}'.format(args.fraction),
+              file=sys.stderr)
+        parser.print_help()
+        return 1
+
+    make_mask(args.topog, args.cice_mask, args.fraction, ice=True)
+    make_mask(args.topog, args.mom_mask, args.fraction)
 
 if __name__ == '__main__':
     sys.exit(main())