Add option to warp surface by upperSurface

Author: matthewhoffman

landice/mesh_tools_li/convert_cell_coordinates_to_sphere.py

@@ -24,6 +24,15 @@
                     help="file to convert to a spherical mesh convention")
 parser.add_argument("-r", dest="radius", type=float, default=6371220.,
                     help="sphere radius to use (m)")
+parser.add_argument("-w", dest="warp_surface", action='store_true',
+                    help="warp surface by upperSurface variable. " \
+                         "Uses scale factor from -s" \
+                         "If upperSurface is not in file, will attempt " \
+                         "to calculate it from thickness and bedTopography. " \
+                         "Uses geometry from first time level only; it is not " \
+                         "possible to have scaling evolve in time with this method.")
+parser.add_argument("-s", dest="scale", type=float, default=100.,
+                    help="scale factor for warping surface")
 args = parser.parse_args()
 
 
@@ -45,7 +54,19 @@ def compute_xyz(latCell_deg, lonCell_deg, radius):
 latCell = np.degrees(ds.variables["latCell"][:])
 lonCell = np.degrees(ds.variables["lonCell"][:])
 
-xCell, yCell, zCell = compute_xyz(latCell, lonCell, args.radius)
+if args.warp_surface:
+    if 'upperSurface' in ds.variables:
+        sfc = ds.variables['upperSurface'][0,:]
+    else:
+        thickness = ds.variables['thickness'][0,:]
+        bedTopography = ds.variables['bedTopography'][0,:]
+        sfc = np.maximum(bedTopography + thickness,
+                         (1.0 - 910.0 / 1028.0) * thickness)
+        #sfc = thickness + np.minimum(bedTopography, 0.0) * 1028.0 / 910.0
+    radius = args.radius + np.maximum(sfc, 0.0) * args.scale
+else:
+    radius = args.radius
+xCell, yCell, zCell = compute_xyz(latCell, lonCell, radius)
 
 # Add to NetCDF file (assumes variables already exist)
 ds.variables["xCell"][:] = xCell