commit 289da8a97695063b56b4ceb0f34e8e57ebafc8c3
parent 57ade89c259cdb4ce92992bab7a0606b68e33f77
Author: Anders Damsgaard <andersd@riseup.net>
Date: Tue, 25 Apr 2017 15:02:10 -0400
use NetCDF files from Baltic test case, interpolate to B grid
Diffstat:
4 files changed, 141 insertions(+), 53 deletions(-)
diff --git a/src/datatypes.jl b/src/datatypes.jl
@@ -91,21 +91,22 @@ type IceFloeArrays
end
#=
-Type containing all relevant data from MOM6 NetCDF file. The ocean grid is a
+Type containing all relevant data from MOM6 NetCDF files. The ocean grid is a
staggered of Arakawa-C type, with north-east convention centered on the
-h-points.
+ h-points. During read, the velocities are interpolated to the cell corners
+ (q-points).
- q(i-1, j) --- v( i, j) --- q( i, j)
+ q(i-1, j) ------------------ q( i, j)
| |
| |
| |
| |
- u(i-1, j) h( i, j) u( i, j)
+ | h( i, j) |
| |
| |
| |
| |
- q(i-1,j-1) --- v( i,j-1) --- q( i,j-1)
+ q(i-1,j-1) ------------------ q( i,j-1)
Source:
https://mom6.readthedocs.io/en/latest/api/generated/pages/Horizontal_indexing.html
@@ -120,13 +121,13 @@ https://mom6.readthedocs.io/en/latest/api/generated/pages/Horizontal_indexing.ht
* `zl::Array{Float64, 1}`: layer target potential density [kg m^-3]
* `zi::Array{Float64, 1}`: interface target potential density [kg m^-3]
* `u::Array{Float64, Int}`: zonal velocity (positive towards west) [m/s],
- dimensions correspond to placement in `[xq, yh, zl, time]`.
+ dimensions correspond to placement in `[xq, yq, zl, time]`.
* `v::Array{Float64, Int}`: meridional velocity (positive towards north) [m/s],
- dimensions correspond to placement in `[xh, yq, zl, time]`.
+ dimensions correspond to placement in `[xq, yq, zl, time]`.
* `h::Array{Float64, Int}`: layer thickness [m], dimensions correspond to
placement in `[xh, yh, zl, time]`.
* `e::Array{Float64, Int}`: interface height relative to mean sea level [m],
- dimensions correspond to placement in `[xh, yq, zi, time]`.
+ dimensions correspond to placement in `[xh, yh, zi, time]`.
=#
type Ocean
input_file::Any
@@ -134,12 +135,12 @@ type Ocean
time::Array{Float64, 1}
# q-point (cell corner) positions
- xq::Array{Float64, 1}
- yq::Array{Float64, 1}
+ xq::Array{Float64, 2}
+ yq::Array{Float64, 2}
# h-point (cell center) positions
- xh::Array{Float64, 1}
- yh::Array{Float64, 1}
+ xh::Array{Float64, 2}
+ yh::Array{Float64, 2}
# Vertical positions
zl::Array{Float64, 1}
diff --git a/src/ocean.jl b/src/ocean.jl
@@ -1,13 +1,17 @@
"Returns empty ocean type for initialization purposes."
function createEmptyOcean()
return Ocean(false,
+
zeros(1),
+
+ zeros(1,1),
+ zeros(1,1),
+ zeros(1,1),
+ zeros(1,1),
+
zeros(1),
zeros(1),
- zeros(1),
- zeros(1),
- zeros(1),
- zeros(1),
+
zeros(1,1,1,1),
zeros(1,1,1,1),
zeros(1,1,1,1),
@@ -15,46 +19,128 @@ function createEmptyOcean()
end
"""
-Read NetCDF file generated by MOM6 (e.g. `prog__####_###.nc`) from disk and
-return as `Ocean` data structure.
+Read ocean NetCDF files generated by MOM6 from disk and return as `Ocean` data
+structure.
+
+# Arguments
+* `velocity_file::String`: Path to NetCDF file containing ocean velocities,
+ etc., (e.g. `prog__####_###.nc`).
+* `grid_file::String`: Path to NetCDF file containing ocean super-grid
+ information (typically `INPUT/ocean_hgrid.nc`).
+"""
+function readOceanNetCDF(velocity_file::String, grid_file::String)
+
+ time, u, v, h, e, zl, zi = readOceanStateNetCDF(velocity_file)
+ xh, yh, xq, yq = readOceanGridNetCDF(grid_file)
+
+ if size(u[:,:,1,1]) != size(xq) || size(v[:,:,1,1]) != size(xq) ||
+ size(xq) != size(yq)
+ error("size mismatch between velocities and grid
+ (u: $(size(u[:,:,1,1])), v: $(size(v[:,:,1,1])),
+ xq: $(size(xq)), yq: $(size(yq)))")
+ end
+
+ ocean = Ocean([grid_file, velocity_file],
+
+ time,
+
+ xq,
+ yq,
+ xh,
+ yh,
+
+ zl,
+ zi,
+
+ u,
+ v,
+ h,
+ e)
+ return ocean
+end
+
+"""
+Read NetCDF file with ocean state generated by MOM6 (e.g. `prog__####_###.nc`
+or `########.ocean_month.nc`) from disk and return time stamps, velocity fields,
+layer thicknesses, interface heights, and vertical coordinates.
+
+# Returns
+* `time::Array{Float, 1}`: Time [s]
+* `u::Array{Float, 2}`: Cell corner zonal velocity [m/s],
+ dimensions correspond to placement in `[xq, yq, zl, time]`
+* `v::Array{Float, 2}`: Cell corner meridional velocity [m/s],
+ dimensions correspond to placement in `[xq, yq, zl, time]`
+* `h::Array{Float64, 2}`: layer thickness [m], dimensions correspond to
+ placement in `[xh, yh, zl, time]`
+* `e::Array{Float64, 2}`: interface height relative to mean sea level [m],
+ dimensions correspond to placement in `[xh, yh, zi, time]`
+* `zl::Array{Float64, 1}`: layer target potential density [kg m^-3]
+* `zi::Array{Float64, 1}`: interface target potential density [kg m^-3]
"""
-function readOceanNetCDF(filename::String)
+function readOceanStateNetCDF(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)
+ u, v = interpolateOceanVelocitiesToCorners(u_staggered, v_staggered)
- ocean = Ocean(filename,
- NetCDF.ncread(filename, "Time"),
+ time = NetCDF.ncread(filename, "time")*24.*60.*60.
+ h = NetCDF.ncread(filename, "h")
+ e = NetCDF.ncread(filename, "e")
- NetCDF.ncread(filename, "xq"),
- NetCDF.ncread(filename, "yq"),
- NetCDF.ncread(filename, "xh"),
- NetCDF.ncread(filename, "yh"),
- NetCDF.ncread(filename, "zl"),
- NetCDF.ncread(filename, "zi"),
+ zl = NetCDF.ncread(filename, "zl")
+ zi = NetCDF.ncread(filename, "zi")
- u,
- v,
- NetCDF.ncread(filename, "h"),
- NetCDF.ncread(filename, "e"))
- return ocean
+ return time, u, v, h, e, zl, zi
+end
+
+"""
+Read NetCDF file with ocean *supergrid* information generated by MOM6 (e.g.
+`ocean_hrid.nc`) from disk and return as `Ocean` data structure. This file is
+located in the simulation `INPUT/` subdirectory.
+
+# Returns
+* `xh::Array{Float, 2}`: Longitude for cell centers [deg]
+* `yh::Array{Float, 2}`: Latitude for cell centers [deg]
+* `xq::Array{Float, 2}`: Longitude for cell corners [deg]
+* `yq::Array{Float, 2}`: Latitude for cell corners [deg]
+"""
+function readOceanGridNetCDF(filename::String)
+
+ if !isfile(filename)
+ error("$(filename) could not be opened")
+ end
+ x::Array{float, 2} = NetCDF.ncread(filename, "x")
+ y::Array{float, 2} = NetCDF.ncread(filename, "y")
+
+ xh = x[2:2:end, 2:2:end]
+ yh = y[2:2:end, 2:2:end]
+
+ xq = x[1:2:end, 1:2:end]
+ yq = y[1:2:end, 1:2:end]
+
+ return xh, yh, xq, yq
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`).
+Convert gridded data from Arakawa-C type (decomposed velocities at faces) to
+Arakawa-B type (velocities at corners) through interpolation.
"""
-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]
+function interpolateOceanVelocitiesToCorners(u_in::Array{float, 4},
+ v_in::Array{float, 4})
+
+ if size(u_in) != size(v_in)
+ error("size of u_in ($(size(u_in))) must match v_in ($(size(v_in)))")
+ end
+
+ nx, ny, nz, nt = size(u_in)
+ #u = Array{float}(nx+1, ny+1, nz, nt)
+ #v = Array{float}(nx+1, ny+1, nz, nt)
+ u = zeros(nx+1, ny+1, nz, nt)
+ v = zeros(nx+1, ny+1, nz, nt)
for i=1:nx
for j=1:ny
if j < ny - 1
diff --git a/test/netcdf.jl b/test/netcdf.jl
@@ -4,14 +4,15 @@
info("#### $(basename(@__FILE__)) ####")
-info("Testing dimensions of content read from prog__0001_006.nc")
-ocean = SeaIce.readOceanNetCDF("prog__0001_006.nc")
-@test length(ocean.xq) == 44
-@test length(ocean.xh) == 44
-@test length(ocean.yq) == 40
-@test length(ocean.yh) == 40
-@test ocean.time ≈ [5., 10.]
-@test size(ocean.u) == (44, 40, 2, 2)
-@test size(ocean.v) == (44, 40, 2, 2)
-@test size(ocean.h) == (44, 40, 2, 2)
-@test size(ocean.e) == (44, 40, 3, 2)
+info("Testing dimensions of content read from Baltic test case")
+ocean = SeaIce.readOceanNetCDF("Baltic/00010101.ocean_month.nc",
+ "Baltic/ocean_hgrid.nc")
+@test ocean.time/(24.*60.*60.) ≈ [.5, 1.5, 2.5, 3.5, 4.5]
+@test size(ocean.xq) == (24, 15)
+@test size(ocean.yq) == (24, 15)
+@test size(ocean.xh) == (23, 14)
+@test size(ocean.yh) == (23, 14)
+@test size(ocean.u) == (24, 15, 63, 5)
+@test size(ocean.v) == (24, 15, 63, 5)
+@test size(ocean.h) == (23, 14, 63, 5)
+@test size(ocean.e) == (23, 14, 64, 5)
diff --git a/test/prog__0001_006.nc b/test/prog__0001_006.nc
Binary files differ.
