tadd read of MOM6 netcdf files - Granular.jl - Julia package for granular dynamics simulation
(HTM) git clone git://src.adamsgaard.dk/Granular.jl
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) LICENSE
---
(DIR) commit ed06106295c2a593d7de8aba8e897e2843767423
(DIR) parent d17a415e6abf470072f60d210616caa78722aa31
(HTM) Author: Anders Damsgaard <andersd@riseup.net>
Date: Fri, 21 Apr 2017 16:02:07 -0400
add read of MOM6 netcdf files
Diffstat:
M src/datatypes.jl | 63 +++++++++++++++++++++++++++++++
M src/io.jl | 30 +++++++++++++++++++++++++++++-
A test/netcdf.jl | 20 ++++++++++++++++++++
A test/prog__0001_006.nc | 0
M test/runtests.jl | 1 +
5 files changed, 113 insertions(+), 1 deletion(-)
---
(DIR) diff --git a/src/datatypes.jl b/src/datatypes.jl
t@@ -100,3 +100,66 @@ type Simulation
overlaps::Array{Array{Float64, 1}, 1}
end
+#=
+Type containing all relevant data from MOM6 NetCDF file. The ocean grid is a
+staggered of Arakawa-C type, with north-east convention centered on the
+h-points.
+
+ q(i-1, j) --- v( i, j) --- q( i, j)
+ | |
+ | |
+ | |
+ | |
+ u(i-1, j) h( i, j) u( i, j)
+ | |
+ | |
+ | |
+ | |
+ q(i-1,j-1) --- v( i,j-1) --- q( i,j-1)
+
+Source:
+https://mom6.readthedocs.io/en/latest/api/generated/pages/Horizontal_indexing.html
+
+# Fields
+* `input_file::String`: path to input NetCDF file
+* `time::Array{Float64, 1}`: time in days
+* `xq::Array{Float64, 1}`: nominal longitude of q-points [degrees_E]
+* `yq::Array{Float64, 1}`: nominal latitude of q-points [degrees_N]
+* `xh::Array{Float64, 1}`: nominal longitude of h-points [degrees_E]
+* `yh::Array{Float64, 1}`: nominal latitude of h-points [degrees_N]
+* `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]`.
+* `v::Array{Float64, Int}`: meridional velocity (positive towards north) [m/s],
+ dimensions correspond to placement in `[xh, 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]`.
+
+
+=#
+type Ocean
+ input_file::String
+
+ time::Array{Float64, 1}
+
+ # q-point (cell corner) positions
+ xq::Array{Float64, 1}
+ yq::Array{Float64, 1}
+
+ # h-point (cell center) positions
+ xh::Array{Float64, 1}
+ yh::Array{Float64, 1}
+
+ # Vertical positions
+ zl::Array{Float64, 1}
+ zi::Array{Float64, 1}
+
+ # Field values
+ u::Array{Float64, 4}
+ v::Array{Float64, 4}
+ h::Array{Float64, 4}
+ e::Array{Float64, 4}
+end
(DIR) diff --git a/src/io.jl b/src/io.jl
t@@ -1,4 +1,5 @@
-import WriteVTK # Install with Pkg.add("WriteVTK")
+import WriteVTK
+import NetCDF
## IO functions
t@@ -67,3 +68,30 @@ function writeVTK(simulation::Simulation;
return nothing
end
end
+
+"""
+Read NetCDF file generated by MOM6 (e.g. `prog__####_###.nc`) from disk and
+return as `Ocean` data structure.
+"""
+function readNetCDF(filename::String)
+
+ if !isfile(filename)
+ error("$(filename) could not be opened")
+ end
+
+ ocean = Ocean(filename,
+ NetCDF.ncread(filename, "Time"),
+
+ NetCDF.ncread(filename, "xq"),
+ NetCDF.ncread(filename, "yq"),
+ NetCDF.ncread(filename, "xh"),
+ NetCDF.ncread(filename, "yh"),
+ NetCDF.ncread(filename, "zl"),
+ NetCDF.ncread(filename, "zi"),
+
+ NetCDF.ncread(filename, "u"),
+ NetCDF.ncread(filename, "v"),
+ NetCDF.ncread(filename, "h"),
+ NetCDF.ncread(filename, "e"))
+ return ocean
+end
(DIR) diff --git a/test/netcdf.jl b/test/netcdf.jl
t@@ -0,0 +1,20 @@
+#!/usr/bin/env julia
+
+# Check if NetCDF files are read correctly from the disk.
+
+import Base.Test
+import SeaIce
+
+info("#### $(basename(@__FILE__)) ####")
+
+info("Testing dimensions of content read from prog__0001_006.nc")
+ocean = SeaIce.readNetCDF("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)
(DIR) diff --git a/test/prog__0001_006.nc b/test/prog__0001_006.nc
Binary files differ.
(DIR) diff --git a/test/runtests.jl b/test/runtests.jl
t@@ -3,3 +3,4 @@ using Base.Test
include("contact-search-and-geometry.jl")
include("collision-2floes-normal.jl")
+include("netcdf.jl")