ocean.jl (15284B)
1 import Pkg 2 using Test 3 using LinearAlgebra 4 5 hasNetCDF = false 6 if VERSION < v"0.7.0-alpha" 7 if typeof(Pkg.installed("NetCDF")) == VersionNumber 8 import NetCDF 9 hasNetCDF = true 10 end 11 else 12 import Pkg 13 if haskey(Pkg.installed(), "NetCDF") 14 import NetCDF 15 hasNetCDF = true 16 end 17 end 18 if !hasNetCDF 19 @warn "Package NetCDF not found. " * 20 "Ocean/atmosphere grid read not supported. " * 21 "If required, install NetCDF and its " * 22 "requirements with `Pkg.add(\"NetCDF\")`." 23 end 24 25 export createEmptyOcean 26 "Returns empty ocean type for initialization purposes." 27 function createEmptyOcean() 28 return Ocean(false, 29 30 zeros(1), 31 32 zeros(1,1), 33 zeros(1,1), 34 zeros(1,1), 35 zeros(1,1), 36 37 zeros(1), 38 zeros(1), 39 40 zeros(1,1,1,1), 41 zeros(1,1,1,1), 42 zeros(1,1,1,1), 43 zeros(1,1,1,1), 44 Array{Array{Int, 1}}(undef, 1, 1), 45 zeros(1,1), 46 1, 1, 1, 1, 47 false, [0.,0.,0.], [1.,1.,1.], [1,1,1], [1.,1.,1.]) 48 end 49 50 export readOceanNetCDF 51 """ 52 Read ocean NetCDF files generated by MOM6 from disk and return as `Ocean` data 53 structure. 54 55 # Arguments 56 * `velocity_file::String`: path to NetCDF file containing ocean velocities, 57 etc., (e.g. `prog__####_###.nc`). 58 * `grid_file::String`: path to NetCDF file containing ocean super-grid 59 information (typically `INPUT/ocean_hgrid.nc`). 60 * `regular_grid::Bool=false`: `true` if the grid is regular (all cells 61 equal and grid is Cartesian) or `false` (default). 62 """ 63 function readOceanNetCDF(velocity_file::String, grid_file::String; 64 regular_grid::Bool=false) 65 66 if !hasNetCDF 67 @warn "Package NetCDF not found. " * 68 "Ocean/atmosphere grid read not supported. " * 69 "Please install NetCDF and its " * 70 "requirements with `Pkg.add(\"NetCDF\")`." 71 else 72 73 time, u, v, h, e, zl, zi = readOceanStateNetCDF(velocity_file) 74 xh, yh, xq, yq = readOceanGridNetCDF(grid_file) 75 76 if size(u[:,:,1,1]) != size(xq) || size(v[:,:,1,1]) != size(xq) || 77 size(xq) != size(yq) 78 error("size mismatch between velocities and grid 79 (u: $(size(u[:,:,1,1])), v: $(size(v[:,:,1,1])), 80 xq: $(size(xq)), yq: $(size(yq)))") 81 end 82 83 ocean = Ocean([grid_file, velocity_file], 84 85 time, 86 87 xq, 88 yq, 89 xh, 90 yh, 91 92 zl, 93 zi, 94 95 u, 96 v, 97 h, 98 e, 99 Array{Array{Int, 1}}(undef, size(xh, 1), size(xh, 2)), 100 zeros(size(xh)), 101 1, 1, 1, 1, 102 103 false, [0.,0.,0.], [1.,1.,1.], [1,1,1], [1.,1.,1.] 104 ) 105 return ocean 106 end 107 end 108 109 export readOceanStateNetCDF 110 """ 111 Read NetCDF file with ocean state generated by MOM6 (e.g. `prog__####_###.nc` 112 or `########.ocean_month.nc`) from disk and return time stamps, velocity fields, 113 layer thicknesses, interface heights, and vertical coordinates. 114 115 # Returns 116 * `time::Vector{Float64}`: Time [s] 117 * `u::Array{Float64, 2}`: Cell corner zonal velocity [m/s], 118 dimensions correspond to placement in `[xq, yq, zl, time]` 119 * `v::Array{Float64, 2}`: Cell corner meridional velocity [m/s], 120 dimensions correspond to placement in `[xq, yq, zl, time]` 121 * `h::Array{Float64, 2}`: layer thickness [m], dimensions correspond to 122 placement in `[xh, yh, zl, time]` 123 * `e::Array{Float64, 2}`: interface height relative to mean sea level [m], 124 dimensions correspond to placement in `[xh, yh, zi, time]` 125 * `zl::Vector{Float64}`: layer target potential density [kg m^-3] 126 * `zi::Vector{Float64}`: interface target potential density [kg m^-3] 127 """ 128 function readOceanStateNetCDF(filename::String) 129 130 if !hasNetCDF 131 @warn "Package NetCDF not found. " * 132 "Ocean/atmosphere grid read not supported. " * 133 "Please install NetCDF and its " * 134 "requirements with `Pkg.add(\"NetCDF\")`." 135 else 136 137 if !isfile(filename) 138 error("$(filename) could not be opened") 139 end 140 141 u_staggered = convert(Array{Float64, 4}, NetCDF.ncread(filename, "u")) 142 v_staggered = convert(Array{Float64, 4}, NetCDF.ncread(filename, "v")) 143 u, v = interpolateOceanVelocitiesToCorners(u_staggered, v_staggered) 144 145 time = convert(Vector{Float64}, 146 NetCDF.ncread(filename, "time") .* 24. * 60. * 60.) 147 h = convert(Array{Float64, 4}, NetCDF.ncread(filename, "h")) 148 e = convert(Array{Float64, 4}, NetCDF.ncread(filename, "e")) 149 150 zl = convert(Vector{Float64}, NetCDF.ncread(filename, "zl")) 151 zi = convert(Vector{Float64}, NetCDF.ncread(filename, "zi")) 152 153 return time, u, v, h, e, zl, zi 154 end 155 end 156 157 export readOceanGridNetCDF 158 """ 159 Read NetCDF file with ocean *supergrid* information generated by MOM6 (e.g. 160 `ocean_hrid.nc`) from disk and return as `Ocean` data structure. This file is 161 located in the simulation `INPUT/` subdirectory. 162 163 # Returns 164 * `xh::Array{Float64, 2}`: Longitude for cell centers [deg] 165 * `yh::Array{Float64, 2}`: Latitude for cell centers [deg] 166 * `xq::Array{Float64, 2}`: Longitude for cell corners [deg] 167 * `yq::Array{Float64, 2}`: Latitude for cell corners [deg] 168 """ 169 function readOceanGridNetCDF(filename::String) 170 171 if !hasNetCDF 172 @warn "Package NetCDF not found. " * 173 "Ocean/atmosphere grid read not supported. " * 174 "Please install NetCDF and its " * 175 "requirements with `Pkg.add(\"NetCDF\")`." 176 else 177 178 if !isfile(filename) 179 error("$(filename) could not be opened") 180 end 181 x = convert(Array{Float64, 2}, NetCDF.ncread(filename, "x")) 182 y = convert(Array{Float64, 2}, NetCDF.ncread(filename, "y")) 183 184 xh = x[2:2:end, 2:2:end] 185 yh = y[2:2:end, 2:2:end] 186 187 xq = x[1:2:end, 1:2:end] 188 yq = y[1:2:end, 1:2:end] 189 190 return xh, yh, xq, yq 191 end 192 end 193 194 export interpolateOceanVelocitiesToCorners 195 """ 196 Convert gridded data from Arakawa-C type (decomposed velocities at faces) to 197 Arakawa-B type (velocities at corners) through interpolation. 198 """ 199 function interpolateOceanVelocitiesToCorners(u_in::Array{Float64, 4}, 200 v_in::Array{Float64, 4}) 201 202 if size(u_in) != size(v_in) 203 error("size of u_in ($(size(u_in))) must match v_in ($(size(v_in)))") 204 end 205 206 nx, ny, nz, nt = size(u_in) 207 u = zeros(nx+1, ny+1, nz, nt) 208 v = zeros(nx+1, ny+1, nz, nt) 209 for i=1:nx 210 for j=1:ny 211 if j < ny - 1 212 u[i, j, :, :] = (u_in[i, j, :, :] + u_in[i, j+1, :, :])/2. 213 else 214 u[i, j, :, :] = u_in[i, j, :, :] 215 end 216 if i < nx - 1 217 v[i, j, :, :] = (v_in[i, j, :, :] + v_in[i+1, j, :, :])/2. 218 else 219 v[i, j, :, :] = v_in[i, j, :, :] 220 end 221 end 222 end 223 return u, v 224 end 225 226 export interpolateOceanState 227 """ 228 Ocean data is containted in `Ocean` type at discrete times (`Ocean.time`). This 229 function performs linear interpolation between time steps to get the approximate 230 ocean state at any point in time. If the `Ocean` data set only contains a 231 single time step, values from that time are returned. 232 """ 233 function interpolateOceanState(ocean::Ocean, t::Float64) 234 if length(ocean.time) == 1 235 return ocean.u, ocean.v, ocean.h, ocean.e 236 elseif t < ocean.time[1] || t > ocean.time[end] 237 error("selected time (t = $(t)) is outside the range of time steps " * 238 "in the ocean data") 239 end 240 241 i = 1 242 rel_time = 0. 243 while i < length(ocean.time) 244 if ocean.time[i+1] < t 245 i += 1 246 continue 247 end 248 249 dt = ocean.time[i+1] - ocean.time[i] 250 rel_time = (t - ocean.time[i])/dt 251 if rel_time < 0. || rel_time > 1. 252 error("time bounds error") 253 end 254 break 255 end 256 257 return ocean.u[:,:,:,i]*(1. - rel_time) + ocean.u[:,:,:,i+1]*rel_time, 258 ocean.v[:,:,:,i]*(1. - rel_time) + ocean.v[:,:,:,i+1]*rel_time, 259 ocean.h[:,:,:,i]*(1. - rel_time) + ocean.h[:,:,:,i+1]*rel_time, 260 ocean.e[:,:,:,i]*(1. - rel_time) + ocean.e[:,:,:,i+1]*rel_time 261 end 262 263 export createRegularOceanGrid 264 """ 265 266 createRegularOceanGrid(n, L[, origo, time, name, 267 bc_west, bc_south, bc_east, bc_north]) 268 269 Initialize and return a regular, Cartesian `Ocean` grid with `n[1]` by `n[2]` 270 cells in the horizontal dimension, and `n[3]` vertical cells. The cell corner 271 and center coordinates will be set according to the grid spatial dimensions 272 `L[1]`, `L[2]`, and `L[3]`. The grid `u`, `v`, `h`, and `e` fields will contain 273 one 4-th dimension matrix per `time` step. Sea surface will be at `z=0.` with 274 the ocean spanning `z<0.`. Vertical indexing starts with `k=0` at the sea 275 surface, and increases downwards. 276 277 # Arguments 278 * `n::Vector{Int}`: number of cells along each dimension [-]. 279 * `L::Vector{Float64}`: domain length along each dimension [m]. 280 * `origo::Vector{Float64}`: domain offset in each dimension [m] (default = 281 `[0.0, 0.0]`). 282 * `time::Vector{Float64}`: vector of time stamps for the grid [s]. 283 * `name::String`: grid name (default = `"unnamed"`). 284 * `bc_west::Integer`: grid boundary condition for the grains. 285 * `bc_south::Integer`: grid boundary condition for the grains. 286 * `bc_east::Integer`: grid boundary condition for the grains. 287 * `bc_north::Integer`: grid boundary condition for the grains. 288 """ 289 function createRegularOceanGrid(n::Vector{Int}, 290 L::Vector{Float64}; 291 origo::Vector{Float64} = zeros(2), 292 time::Vector{Float64} = zeros(1), 293 name::String = "unnamed", 294 bc_west::Integer = 1, 295 bc_south::Integer = 1, 296 bc_east::Integer = 1, 297 bc_north::Integer = 1) 298 299 xq = repeat(range(origo[1], stop=origo[1] + L[1], 300 length=n[1] + 1), 301 outer=[1, n[2] + 1]) 302 yq = repeat(range(origo[2], stop=origo[2] + L[2], 303 length=n[2] + 1)', 304 outer=[n[1] + 1, 1]) 305 306 dx = L./n 307 xh = repeat(range(origo[1] + .5*dx[1], 308 stop=origo[1] + L[1] - .5*dx[1], 309 length=n[1]), 310 outer=[1, n[2]]) 311 yh = repeat(range(origo[2] + .5*dx[2], 312 stop=origo[2] + L[2] - .5*dx[2], 313 length=n[2])', 314 outer=[n[1], 1]) 315 316 zl = -range(.5*dx[3], stop=L[3] - .5*dx[3], length=n[3]) 317 zi = -range(0., stop=L[3], length=n[3] + 1) 318 319 u = zeros(n[1] + 1, n[2] + 1, n[3], length(time)) 320 v = zeros(n[1] + 1, n[2] + 1, n[3], length(time)) 321 h = zeros(n[1] + 1, n[2] + 1, n[3], length(time)) 322 e = zeros(n[1] + 1, n[2] + 1, n[3], length(time)) 323 324 return Ocean(name, 325 time, 326 xq, yq, 327 xh, yh, 328 zl, zi, 329 u, v, h, e, 330 Array{Array{Int, 1}}(undef, size(xh, 1), size(xh, 2)), 331 zeros(size(xh)), 332 bc_west, bc_south, bc_east, bc_north, 333 true, origo, L, n, dx) 334 end 335 336 export addOceanDrag! 337 """ 338 Add drag from linear and angular velocity difference between ocean and all ice 339 floes. 340 """ 341 function addOceanDrag!(simulation::Simulation) 342 if typeof(simulation.ocean.input_file) == Bool 343 error("no ocean data read") 344 end 345 346 u, v, h, e = interpolateOceanState(simulation.ocean, simulation.time) 347 uv_interp = Vector{Float64}(undef, 2) 348 sw = Vector{Float64}(undef, 2) 349 se = Vector{Float64}(undef, 2) 350 ne = Vector{Float64}(undef, 2) 351 nw = Vector{Float64}(undef, 2) 352 353 for grain in simulation.grains 354 355 if !grain.enabled 356 continue 357 end 358 359 i, j = grain.ocean_grid_pos 360 k = 1 361 362 x_tilde, y_tilde = getNonDimensionalCellCoordinates(simulation.ocean, 363 i, j, 364 grain.lin_pos) 365 x_tilde = clamp(x_tilde, 0., 1.) 366 y_tilde = clamp(y_tilde, 0., 1.) 367 368 bilinearInterpolation!(uv_interp, u, v, x_tilde, y_tilde, i, j, k, 1) 369 applyOceanDragToGrain!(grain, uv_interp[1], uv_interp[2]) 370 applyOceanVorticityToGrain!(grain, 371 curl(simulation.ocean, x_tilde, y_tilde, 372 i, j, k, 1, sw, se, ne, nw)) 373 end 374 nothing 375 end 376 377 export applyOceanDragToGrain! 378 """ 379 Add Stokes-type drag from velocity difference between ocean and a single ice 380 floe. 381 """ 382 function applyOceanDragToGrain!(grain::GrainCylindrical, 383 u::Float64, v::Float64) 384 freeboard = .1*grain.thickness # height above water 385 ρ_o = 1000. # ocean density 386 draft = grain.thickness - freeboard # height of submerged thickness 387 388 drag_force = ρ_o * π * 389 (2.0*grain.ocean_drag_coeff_vert*grain.areal_radius*draft + 390 grain.ocean_drag_coeff_horiz*grain.areal_radius^2.0) * 391 ([u, v] - grain.lin_vel[1:2])*norm([u, v] - grain.lin_vel[1:2]) 392 393 grain.force += vecTo3d(drag_force) 394 grain.ocean_stress = vecTo3d(drag_force/grain.horizontal_surface_area) 395 nothing 396 end 397 398 export applyOceanVorticityToGrain! 399 """ 400 Add Stokes-type torque from angular velocity difference between ocean and a 401 single grain. See Eq. 9.28 in "Introduction to Fluid Mechanics" by Nakayama 402 and Boucher, 1999. 403 """ 404 function applyOceanVorticityToGrain!(grain::GrainCylindrical, 405 ocean_curl::Float64) 406 freeboard = .1*grain.thickness # height above water 407 ρ_o = 1000. # ocean density 408 draft = grain.thickness - freeboard # height of submerged thickness 409 410 grain.torque[3] += 411 π * grain.areal_radius^4. * ρ_o * 412 (grain.areal_radius/5. * grain.ocean_drag_coeff_horiz + 413 draft * grain.ocean_drag_coeff_vert) * 414 norm(.5 * ocean_curl - grain.ang_vel[3]) * 415 (.5 * ocean_curl - grain.ang_vel[3]) 416 nothing 417 end 418 419 export compareOceans 420 """ 421 compareOceans(ocean1::Ocean, ocean2::Ocean) 422 423 Compare values of two `Ocean` objects using the `Base.Test` framework. 424 """ 425 function compareOceans(ocean1::Ocean, ocean2::Ocean) 426 427 @test ocean1.input_file == ocean2.input_file 428 @test ocean1.time ≈ ocean2.time 429 430 @test ocean1.xq ≈ ocean2.xq 431 @test ocean1.yq ≈ ocean2.yq 432 433 @test ocean1.xh ≈ ocean2.xh 434 @test ocean1.yh ≈ ocean2.yh 435 436 @test ocean1.zl ≈ ocean2.zl 437 @test ocean1.zi ≈ ocean2.zi 438 439 @test ocean1.u ≈ ocean2.u 440 @test ocean1.v ≈ ocean2.v 441 @test ocean1.h ≈ ocean2.h 442 @test ocean1.e ≈ ocean2.e 443 444 if isassigned(ocean1.grain_list, 1) 445 @test ocean1.grain_list == ocean2.grain_list 446 end 447 nothing 448 end