Remaining issues with hydro files

[jessmt]

Need to:

1. Add vertical exchanges
2. Figure out why in the time series visualisation in Shiny we see elements of vertical exchange even when vertical exchange is excluded from plotting
3. Figure out why the magnitude of horizontal exchange is still relatively small (should be seeing exchange for some boxes in the order of 0.1 - 1.)

Have asked Javier to help look at the sediment particle tracking to validate exchanges in Atlantis.

@mdsumner @mmori1

[mdsumner]

When looking at the bottoming out via botz it changes boxes 1-6 to maximum level 7. (all others go to 10, including box 0)

And we get this tally of raw ROMS data (for a given time step) , different between the two cases

romstab %>% dplyr::select(-temp, -salt) %>% 
+   inner_join(face_index, c("cellindex" = "cell_")) %>% 
+   inner_join(box_index[c("cell_", "maxlayer")], c("cellindex" = "cell_")) %>% 
+   dplyr::mutate(layer = pmin(layer, maxlayer)) %>% count(layer)
# A tibble: 10 x 2
   layer     n
   <int> <int>
 1     1  4741
 2     2  2607
 3     3  2902
 4     4  4279
 5     5  3351
 6     6  2520
 7     7  7713
 8     8  3092
 9     9  9598
10    10 16640
> romstab %>% dplyr::select(-temp, -salt) %>% 
+   inner_join(face_index, c("cellindex" = "cell_")) %>% 
+   inner_join(box_index[c("cell_", "maxlayer")], c("cellindex" = "cell_")) %>% count(layer)
# A tibble: 10 x 2
   layer     n
   <int> <int>
 1     1  4741
 2     2  2607
 3     3  2902
 4     4  4279
 5     5  3351
 6     6  2520
 7     7  6348
 8     8  3572
 9     9 10438
10    10 16685

[mdsumner]

Here's some quick plotting (which immediately shows that verticalflux still needs sum/scaling ...). But, should be simple to modify this to compare directly in subsequent runs.

  fp <- "~/Git/EA-Atlantis-dev/hydroconstruct/hydroruns.6"

ncfiles <- list.files(fp, pattern  = "nc$", full.names = TRUE)
library(raster)
#> Loading required package: sp
## super quick plots by time
transport <- brick(grep("transport", ncfiles, value = TRUE), stopIfNotEqualSpaced = FALSE)
#> Loading required namespace: ncdf4
## 16 days throughout the year
plot(subset(transport, seq(1, nlayers(transport), length = 16)), asp = NA)

temp <- brick(grep("mass", ncfiles, value = TRUE), varname = "temperature", stopIfNotEqualSpaced = FALSE)
salt <- brick(grep("mass", ncfiles, value = TRUE), varname = "salinity", stopIfNotEqualSpaced = FALSE)
vert <- brick(grep("mass", ncfiles, value = TRUE), varname = "verticalflux", stopIfNotEqualSpaced = FALSE)
## 16 days throughout the year
plot(subset(temp, seq(1, nlayers(temp), length = 16)), asp = NA)

plot(subset(salt, seq(1, nlayers(salt), length = 16)), asp = NA)

plot(subset(vert, seq(1, nlayers(vert), length = 16)), asp = NA)

^{Created on 2019-01-24 by the reprex package (v0.2.1)}

[mdsumner]

run 7 (scale/sum vertical flux now)

  fp <- "~/Git/EA-Atlantis-dev/hydroconstruct/hydroruns.7"

ncfiles <- list.files(fp, pattern  = "nc$", full.names = TRUE)
library(raster)
#> Loading required package: sp
## super quick plots by time
transport <- brick(grep("transport", ncfiles, value = TRUE), stopIfNotEqualSpaced = FALSE)
#> Loading required namespace: ncdf4
## 16 days throughout the year
plot(subset(transport, seq(1, nlayers(transport), length = 16)), asp = NA)

temp <- brick(grep("mass", ncfiles, value = TRUE), varname = "temperature", stopIfNotEqualSpaced = FALSE)
salt <- brick(grep("mass", ncfiles, value = TRUE), varname = "salinity", stopIfNotEqualSpaced = FALSE)
vert <- brick(grep("mass", ncfiles, value = TRUE), varname = "verticalflux", stopIfNotEqualSpaced = FALSE)
## 16 days throughout the year
plot(subset(temp, seq(1, nlayers(temp), length = 16)), asp = NA)

plot(subset(salt, seq(1, nlayers(salt), length = 16)), asp = NA)

plot(subset(vert, seq(1, nlayers(vert), length = 16)), asp = NA)

^{Created on 2019-01-24 by the reprex package (v0.2.1)}

[jessmt]

Need to make sure that the velocities are multiplied by area so that flow is in units of m3/s (Jess is confirming with Javier that this is NOT done by hydroconstruct - I don't think it is but for some reason David's code stopped at average velocities across faces and didn't multiply by area). There seem to be two alternative ways to do this: (i) calculate the mean velocity across a face and then multiply by the area of the face (ii) sum the flows (velocity x area) for each ROMS layer that intersects (Jess has asked Javier if one of these is preferable). Next week need to check the magnitudes of the flows in the raw output versus the output from hydroconstruct to check we understand how hyperdiffusion is corrected for (from the wiki it's just a division by the area of each box) and the relative magnitudes of the two outputs.

[mdsumner]

I think (i) is relatively straightforward from current script, but I have to unpick the face/layer arrangement a little. (These are just notes for self, for now).

EDIT: repeat after me, faces are segments.

• add a segment id, we need the identity of every segment - not just face - and corresponds to each x0, y0, x1, y1
• delay creating a face_index and carrying that into the roms data
• within ROMS time step group by segment_id, layer and get the mean flow
• multiply each flow by segment length and then join on face attributes
• multiply each (flow * segment length) by layer depth, then group by face id and get overall mean

[mdsumner]

Forget that last post. My problem was the different bottoming out of layers, I now directly index the dz based on the maxlayer determined for each face. I was getting more than 900 groups (90 faces 10 layers) sometimes because of this layer/maxlayer decoupling and so the indexes didn't check. Now running hydroruns.8*

[mdsumner]

Hmm, numbers are quite a lot bigger now

transport: