Adding Alistair's grid conversion program.

- Make a MOM6 grid from a ROMS grid.

README.md

@@ -136,14 +136,14 @@ $ make DEVELOP=1 uninstall
 
 ## Running Pyroms
 
-We have a gridid.txt file that's pointed to by the PYROMS_GRIDID_FILE
+We have a gridid.txt file that's pointed to by the PYROMS\_GRIDID\_FILE
 environment variable. If you are operating on files containing
 sufficient grid information already, you won't need to use this.
 An example is provided in the examples directory.
 
 
 ## Doxygen
 
-Running "doxygen .doxygen" in any of pyroms, pyroms_toolbox or
-bathy_smoother will generate doxygen files. Edit the .doxygen files to
+Running "doxygen .doxygen" in any of pyroms, pyroms\_toolbox or
+bathy\_smoother will generate doxygen files. Edit the .doxygen files to
 specify html vs. some other output format.

examples/grid_MOM6/README

@@ -0,0 +1,8 @@
+These scripts came from Mehmet Ilicak who says the conversion script came from Alistair Adcroft.
+
+When running this on my pan-Arctic grid, I got:
+
+convert_ROMS_grid_to_MOM6.py:230: RuntimeWarning: divide by zero encountered in true_divide
+  mom6_grid['supergrid']['angle'][:,1:-1] = numpy.arctan( (lat[:,2:] - lat[:,:-2]) / ((lon[:,2:] - lon[:,:-2]) * cos_lat[:,1:-1]) )
+
+It seemed to work anyway.

examples/grid_MOM6/Spherical.py

@@ -0,0 +1,51 @@
+# -*- coding: utf-8 -*-
+# vim: set fileencoding=utf-8 :
+
+# Code to compute distances and angles on a sphere.
+# Based on code written by Alistair Adcroft and Matthew Harrison of GFDL
+
+import numpy
+
+def angle_through_center(p1, p2):
+    """Angle at center of sphere between two points on the surface of the sphere.
+    Positions are given as (latitude,longitude) tuples measured in degrees."""
+    phi1 = numpy.deg2rad( p1[0] )
+    phi2 = numpy.deg2rad( p2[0] )
+    dphi_2 = 0.5 * ( phi2 - phi1 )
+    dlambda_2 = 0.5 * numpy.deg2rad( p2[1] - p1[1] )
+    a = numpy.sin( dphi_2 )**2 + numpy.cos( phi1 ) * numpy.cos( phi2 ) * ( numpy.sin( dlambda_2 )**2 )
+    c = 2. * numpy.arctan2( numpy.sqrt(a), numpy.sqrt( 1. - a ) )
+    return c
+
+def angle_between(v1, v2, v3):
+    """Returns angle v2-v1-v3 i.e betweeen v1-v2 and v1-v3."""
+    # vector product between v1 and v2
+    px = v1[1] * v2[2] - v1[2] * v2[1]
+    py = v1[2] * v2[0] - v1[0] * v2[2]
+    pz = v1[0] * v2[1] - v1[1] * v2[0]
+    # vector product between v1 and v3
+    qx = v1[1] * v3[2] - v1[2] * v3[1]
+    qy = v1[2] * v3[0] - v1[0] * v3[2]
+    qz = v1[0] * v3[1] - v1[1] * v3[0]
+
+    ddd = (px * px + py * py + pz * pz) * (qx * qx + qy * qy + qz * qz)
+    ddd = (px * qx + py * qy + pz * qz) / numpy.sqrt(ddd)
+    angle = numpy.arccos( ddd )
+    return angle
+
+def quad_area(lat, lon):
+    """Returns area of spherical quad (bounded by great arcs)."""
+    # x,y,z are 3D coordinates
+    d2r = numpy.deg2rad(1.)
+    x = numpy.cos(d2r * lat) * numpy.cos(d2r * lon)
+    y = numpy.cos(d2r * lat) * numpy.sin(d2r * lon)
+    z = numpy.sin(d2r * lat)
+    c0 = (x[ :-1, :-1], y[ :-1, :-1], z[ :-1, :-1])
+    c1 = (x[ :-1,1:  ], y[ :-1,1:  ], z[ :-1,1:  ])
+    c2 = (x[1:  ,1:  ], y[1:  ,1:  ], z[1:  ,1:  ])
+    c3 = (x[1:  , :-1], y[1:  , :-1], z[1:  , :-1])
+    a0 = angle_between(c1, c0, c2)
+    a1 = angle_between(c2, c1, c3)
+    a2 = angle_between(c3, c2, c0)
+    a3 = angle_between(c0, c3, c1)
+    return a0 + a1 + a2 + a3 - 2. * numpy.pi