Merge pull request #450 from hollyhan/fix_MPAS_mesh_tools

Fix mpas mesh tools

conda_package/mpas_tools/scrip/from_mpas.py

@@ -6,7 +6,7 @@
 import numpy as np
 
 from optparse import OptionParser
-
+from mpas_tools.cime.constants import constants
 
 def scrip_from_mpas(mpasFile, scripFile, useLandIceMask=False):
     """
@@ -47,10 +47,17 @@ def scrip_from_mpas(mpasFile, scripFile, useLandIceMask=False):
     sphereRadius = float(fin.sphere_radius)
     on_a_sphere = str(fin.on_a_sphere)
 
+    # check the longitude convention to use positive values [0 2pi]
+    if np.any(np.logical_or(lonCell < 0, lonCell > 2.0 * np.pi)):
+       raise ValueError("lonCell is not in the desired range (0, 2pi)")
+
+    if np.any(np.logical_or(lonVertex < 0, lonVertex > 2.0 * np.pi)):
+       raise ValueError("lonVertex is not in the desired range (0, 2pi)")
+
     if sphereRadius <= 0:
-        print(" -- WARNING: conservative remapping is NOT possible when "
-              "'sphereRadius' <= 0 because 'grid_area' field will be infinite "
-              "(from division by 0)")
+       sphereRadius =  constants['SHR_CONST_REARTH']
+       print(f" -- WARNING: sphereRadius<0 so setting sphereRadius = "
+             f"{constants['SHR_CONST_REARTH']}")
 
     if on_a_sphere == "NO":
         print(" -- WARNING: 'on_a_sphere' attribute is 'NO', which means that "

mesh_tools/planar_grid_transformations/set_lat_lon_fields_in_planar_grid.py

@@ -6,6 +6,7 @@
 from __future__ import absolute_import, division, print_function, unicode_literals
 
 import sys
+import numpy as np
 import netCDF4
 import pyproj
 from optparse import OptionParser
@@ -68,24 +69,36 @@
 xEdge = f.variables['xEdge']
 yEdge = f.variables['yEdge']
 
-latCell = f.variables['latCell']
-lonCell = f.variables['lonCell']
-latVertex = f.variables['latVertex']
-lonVertex = f.variables['lonVertex']
-latEdge = f.variables['latEdge']
-lonEdge = f.variables['lonEdge']
+latCellVar = f.variables['latCell']
+lonCellVar = f.variables['lonCell']
+latVertexVar = f.variables['latVertex']
+lonVertexVar = f.variables['lonVertex']
+latEdgeVar = f.variables['latEdge']
+lonEdgeVar = f.variables['lonEdge']
 
 print("Input file xCell min/max values:", xCell[:].min(), xCell[:].max())
 print("Input file yCell min/max values:", yCell[:].min(), yCell[:].max())
 
 # populate x,y fields
 # MPAS uses lat/lon in radians, so have pyproj return fields in radians.
-lonCell[:], latCell[:] = pyproj.transform(projections[options.projection], projections['latlon'], xCell[:], yCell[:], radians=True)
-lonVertex[:], latVertex[:] = pyproj.transform(projections[options.projection], projections['latlon'], xVertex[:], yVertex[:], radians=True)
-lonEdge[:], latEdge[:] = pyproj.transform(projections[options.projection], projections['latlon'], xEdge[:], yEdge[:], radians=True)
-
-print("Calculated latCell min/max values (radians):", latCell[:].min(), latCell[:].max())
-print("Calculated lonCell min/max values (radians):", lonCell[:].min(), lonCell[:].max())
+lonCell, latCell = pyproj.transform(projections[options.projection], projections['latlon'], xCell[:], yCell[:], radians=True)
+lonVertex, latVertex = pyproj.transform(projections[options.projection], projections['latlon'], xVertex[:], yVertex[:], radians=True)
+lonEdge, latEdge = pyproj.transform(projections[options.projection], projections['latlon'], xEdge[:], yEdge[:], radians=True)
+
+# change the longitude convention to use positive values [0 2pi]
+lonCell = np.mod(lonCell, 2.0*np.pi)
+lonVertex = np.mod(lonVertex, 2.0*np.pi)
+lonEdge = np.mod(lonEdge, 2.0*np.pi)
+
+print("Calculated latCell min/max values (radians):", latCell.min(), latCell.max())
+print("Calculated lonCell min/max values (radians):", lonCell.min(), lonCell.max())
+
+latCellVar[:] = latCell
+lonCellVar[:] = lonCell
+latVertexVar[:] = latVertex
+lonVertexVar[:] = lonVertex
+latEdgeVar[:] = latEdge
+lonEdgeVar[:] = lonEdge
 
 # Update history attribute of netCDF file
 thiscommand = datetime.now().strftime("%a %b %d %H:%M:%S %Y") + ": " + " ".join(sys.argv[:])