Granular.jl

commit 89b95b4d0410273d568677bcb6242c49b1019ef1
parent 764bf1a96f0c1c5defe0afe0cc35c1db7f59eb07
Author: Anders Damsgaard <andersd@riseup.net>
Date:   Mon, 24 Apr 2017 21:18:50 -0400

interpolate ocean velocities to cell center

Diffstat:
M
src/ocean.jl
 | 
35
+++++++++++++++++++++++++++++++++--

1 file changed, 33 insertions(+), 2 deletions(-)
diff --git a/src/ocean.jl b/src/ocean.jl
@@ -23,6 +23,9 @@ function readOceanNetCDF(filename::String)
     if !isfile(filename)
         error("$(filename) could not be opened")
     end
+    u_staggered::Array{float, 4} = NetCDF.ncread(filename, "u")
+    v_staggered::Array{float, 4} = NetCDF.ncread(filename, "v")
+    u, v = convertToColocatedOceanGrid(u_staggered, v_staggered)
 
     ocean = Ocean(filename,
                   NetCDF.ncread(filename, "Time"),
@@ -34,9 +37,37 @@ function readOceanNetCDF(filename::String)
                   NetCDF.ncread(filename, "zl"),
                   NetCDF.ncread(filename, "zi"),
 
-                  NetCDF.ncread(filename, "u"),
-                  NetCDF.ncread(filename, "v"),
+                  u,
+                  v,
                   NetCDF.ncread(filename, "h"),
                   NetCDF.ncread(filename, "e"))
     return ocean
 end
+
+"""
+Convert gridded data from staggered (Arakawa-C) to collocated grid (Arakawa-A) 
+through interpolation.  The new data points are located in the centers of the 
+original staggered grid (spatial coordinates `xh` and `yh`).
+"""
+function convertToColocatedOceanGrid(u_in::Array{float, 4},
+                                     v_in::Array{float, 4})
+    u = Array{float}(size(u_in))
+    v = Array{float}(size(v_in))
+    nx = size(u_in)[1]
+    ny = size(u_in)[2]
+    for i=1:nx
+        for j=1:ny
+            if j < ny - 1
+                u[i, j, :, :] = (u_in[i, j, :, :] + u_in[i, j+1, :, :])/2.
+            else
+                u[i, j, :, :] = u_in[i, j, :, :]
+            end
+            if i < nx - 1
+                v[i, j, :, :] = (v_in[i, j, :, :] + v_in[i+1, j, :, :])/2.
+            else
+                v[i, j, :, :] = v_in[i, j, :, :]
+            end
+        end
+    end
+    return u, v
+end