Script used to load, model, and plot Sargassum trajectories in the North Atlantic for McGillicuddy et al.

Sargassum Trajectories.R

library(TheSource)
source('R/Particle Trajectories v2 functions.R')

stn = readRDS('stn.rds') # Load sampling site information

## Setup model
n = 500
particles = init.lagrangian.particles(lon = rep(stn$Lon, n), lat = rep(stn$Lat, n), time = rep(stn$Datetime, n))
particles$Station = rep(stn$Station, n)

## Load the wind vectors once:
u.wind = load.nc('A:/NCEP GDAS/wind/uwnd.sfc.2021.nc', verbose = F)
v.wind = load.nc('A:/NCEP GDAS/wind/vwnd.sfc.2021.nc', verbose = F)
u.wind$lon = seq(0, 357.5, length.out = 144)
u.wind$time = make.time(1800) + 86400 * u.wind$time
u.wind$grid = expand.grid(lon = u.wind$lon, lat = u.wind$lat)

windage = 0.02
model = init.lagrangian.model(t.start = make.time(2021, 1, 7), t.end = max(sargas$Datetime) + 86400, t.step = -3600*6, nparticles = nrow(particles), save.freq = 1)
  
## Load data products
advection = load.advection.oscar('Z:/Data/OSCAR/oscar_vel2021.nc', lats = c(0, 50), lons = c(-100, -20)+360)
  
## Apply Windage
index = index.grid(x = u.wind$grid$lon, y = u.wind$grid$lat, x.new = advection$grid$lon, y.new = advection$grid$lat)
  
for (k in 1:length(index)) {
  i = ((k - 1) %% length(advection$lon)) + 1
  j = floor((k - 1) / length(advection$lon)) + 1
    
  for (t in 1:length(advection$time)) {
    wind.t = which.min(abs(advection$time[t] - u.wind$time))
    wind.i = ((index[k] - 1) %% length(u.wind$lon)) + 1
    wind.j = floor((index[k]-1) / length(u.wind$lon)) + 1
      
    advection$u[i,j,t] = advection$u[i,j,t] + windage * u.wind$uwnd[wind.i,wind.j,wind.t]
    advection$v[i,j,t] = advection$v[i,j,t] + windage * v.wind$vwnd[wind.i,wind.j,wind.t]
  }
}


## Run model
model = run.advection(particles = particles, model = model, advection = advection, Kd = 4000)
saveRDS(model, file = paste0('_rdata/Model v3 (', windage, ' windage - OSCAR 4k).rds'))


## Save Figure
png(paste0('_figures/Back Traj for Sargassum Data v3 (', windage, ' windage - OSCAR 4k).png'), width = 2400, height = 1800, res = 250)
    
p = make.proj('aea', lon = 310, lat = 30)
map = make.map2('coastline3', lon = 300, lat = 25, scale = 2500, land.col = '#383838', dlon = 20, dlat = 20, p = p)
add.map.axis(map, sides = c(1:4), N = 1e4)

for (i in 1:nrow(stn)) {
      
  l = which(particles$Station == stn$Station[i])
  k = which(cumsum(model$hist[l[1], ,9]) * abs(model$meta$t.step) * model$meta$save.freq / 86400 < 60) # Select the length of time [days] to show the tails.
      
  for (j in l) {
    add.map.line(map, lon = model$hist[j,k,1], lat = model$hist[j,k,2], lwd = 2, greatCircle = F, col = '#22229902')
  }
  add.map.line(map, lon = apply(model$hist[l,k,1], 2, function(x){mean(x, na.rm = T)}), lat = apply(model$hist[l,k,2], 2, function(x){mean(x, na.rm = T)}), lwd = 2, greatCircle = F, col = '#2244bb90')
}
add.map.scale(pos = 1)
    
add.map.points(map, stn$Lon, stn$Lat, pch = 16, cex = 2, col = '#00000099')