# Example of usage of the Chesapeake Bay atlas (https://doi.org/10.17882/99441)
# MACAN webinar of 2024-09-05, pst-laurent@vims.edu
#
# Use package 'seacarb' to derive carbonate chemistry variables such as
#   saturation states and pCO2 from the alkalinity and DIC available in the
#   atlas.
 
# Import the different libraries that we will need.
# If a library/package isn't already installed on your computer, type:
#   install.packages( 'name_of_library' )

  library( ncdf4        ) # To read NetCDF files.
  library( seacarb      ) # To model the inorganic carbon system.
 
# Read from the atlas the different variables that we will need. We specifically
#   select bottom fields since this is where bivalves live.

  atla <- nc_open( 'atlas_chesbay_v20240319.nc' )

  long <- ncvar_get( atla, 'longitude'          ) # Dims: 348
  lati <- ncvar_get( atla, 'latitude'           ) # Dims:     567
  salt <- ncvar_get( atla, 'salinity_bottom'    ) # Dims: 348x567x12 (PSU    ).
  temp <- ncvar_get( atla, 'temperature_bottom' ) # Dims: 348x567x12 (deg.C  ).
  dic  <- ncvar_get( atla, 'DIC_bottom'         ) # Dims: 348x567x12 (umol/kg).
  ta   <- ncvar_get( atla, 'TA_bottom'          ) # Dims: 348x567x12 (umol/kg).

  nc_close( atla )                                # Close file once done.
  rm(       atla )                                # No longer necessary.

# We will focus on Mobjack Bay, and more specifically the monitoring station
#   WE4.1 which is roughly at the mouth of that Bay in water of depth ~6meters.

  lo_s <- - 76.34634 # LOngitude of Station.
  la_s <-   37.31181 # LAtitude  of Station.

# Spatially interpolate the atlas' fields at this particular station.
#   Begin by noting the DIMensions of the atlas' Fields (fdim), then find the
#   four points that surround the station, and finally, conduct a bilinear
#   interpolation, for each field of interest, and each MONTh of the year.

  fdim <- dim( salt ) # 348, 567, 12.

# Four neighbors surrounding the station: (i1,j1),(i2,j1),(i1,j2),(i2,j2).
#   (0<x<1, 0<y<1) is the station's position relative to the 4 neighbors.

  i    <- approx( long, seq( from = 1, to = fdim[1] ), lo_s )
  j    <- approx( lati, seq( from = 1, to = fdim[2] ), la_s )
  i1   <- floor( i$y )
  i2   <- i1 + 1
  j1   <- floor( j$y )
  j2   <- j1 + 1
  x    <- (lo_s - long[ i1 ]) / (long[ i2 ] - long[ i1 ])
  y    <- (la_s - lati[ j1 ]) / (lati[ j2 ] - lati[ j1 ])
  rm( fdim, i, j, long, lati ) # No longer necessary.

  timeseries_salt <- rep( NA, times = 12 ) # Empty containers that we will fill.
  timeseries_temp <- rep( NA, times = 12 )
  timeseries_dic  <- rep( NA, times = 12 )
  timeseries_ta   <- rep( NA, times = 12 )

# Loop over the months of year and do the bilinear interpolation:

  for ( mont in 1 : 12 ) {
    timeseries_salt[ mont ] = salt[ i1, j1, mont ] * (1. - x) * (1. - y) +
                              salt[ i1, j2, mont ] * (1. - x) *       y  +
                              salt[ i2, j1, mont ] *       x  * (1. - y) +
                              salt[ i2, j2, mont ] *       x  *       y
    timeseries_temp[ mont ] = temp[ i1, j1, mont ] * (1. - x) * (1. - y) +
                              temp[ i1, j2, mont ] * (1. - x) *       y  +
                              temp[ i2, j1, mont ] *       x  * (1. - y) +
                              temp[ i2, j2, mont ] *       x  *       y
    timeseries_dic [ mont ] = dic [ i1, j1, mont ] * (1. - x) * (1. - y) +
                              dic [ i1, j2, mont ] * (1. - x) *       y  +
                              dic [ i2, j1, mont ] *       x  * (1. - y) +
                              dic [ i2, j2, mont ] *       x  *       y
    timeseries_ta  [ mont ] = ta  [ i1, j1, mont ] * (1. - x) * (1. - y) +
                              ta  [ i1, j2, mont ] * (1. - x) *       y  +
                              ta  [ i2, j1, mont ] *       x  * (1. - y) +
                              ta  [ i2, j2, mont ] *       x  *       y
  }
  rm( salt, temp, dic, ta, i1, i2, j1, j2, x, y ) # No longer necessary.

# Use 'seacarb' to compute the full seawater Carbonate SYStem, using TA, DIC,
#   salinity and temperature as the four inputs.

  csys <- carb( flag    = 15,                     # We are providing TA and DIC (see doc).
                var1    = timeseries_ta  * 1.e-6, # TA,  converted from umol/kg to mol/kg.
                var2    = timeseries_dic * 1.e-6, # DIC, converted from umol/kg to mol/kg.
                S       = timeseries_salt,        # Seawater practical salinity.
                T       = timeseries_temp,        # Seawater potential temperature (deg.C).
                Patm    = 1.,                     # Atmospheric pressure (default = 1atm).
                P       = 0.6,                    # Bottom of WE4.1 = 6meters = 0.6bar.
                Pt      = 0.,                     # Total phosphate concentration (mol/kg).
                Sit     = 0.,                     # Total silicate  concentration (mol/kg).
                k1k2    = 'cw',                   # Diss.constants from Cai & Wang 1998.
                kf      = 'x',                    # Default value for Kf.
                ks      = 'd',                    # Default value for Ks (Dickson 1990).
                pHscale = 'T',                    # Default value is Total Scale.
                b       = 'u74',                  # Default value for boron concentration.
                gas     = 'potential',            # (Irrelevant unless pCO2 is an input.)
                warn    = 'y',                    # Yes, show me warnings (if any).
                eos     = 'eos80',                # When providing potential T and practical S.
                long    = lo_s,                   # Longitude (irrelevant when using eos80).
                lat     = la_s )                  # Latitude  (irrelevant when using eos80).

# The first figure will show TA and DIC, with TA on the left axis and DIC on
#   the right axis. Modify the default plot margins to accomodate the extra
#   labels on the right axis:

  par( mar = c( 5.25, 4.25, 4.25, 4.25 ) )

# TA as a function of the months of the year; discard the default labeling of
#   the x axis (xaxt='n', xlab='') as you will enter your own (see below). Set
#   the range of the left vertical axis to 1720-1800 umol/kg:

  plot( seq( from = 1, to = 12 ), timeseries_ta, type = 'l', col = 'blue',
    lwd = 5, xlab = '', ylab = 'TA (umol kg-1)', col.lab = 'blue', xaxt = 'n',
    ylim = c( 1720, 1800 ) )

# Set your own labels for the horizontal axis (12 months of the year):

  axis( 1, at = seq( from = 1, to = 12, by = 1 ),
    labels = c('Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec') )

# Add DIC on the right axis of the same figure, colored in red:

  par( new = TRUE )

  plot( seq( from = 1, to = 12 ), timeseries_dic, type = 'l', col = 'red',
    lwd = 5, axes = FALSE, bty = 'n', xlab = '', ylab = '', ylim = c(1590, 1670) )

  axis( 4, at = seq( from = 1590, to = 1670, by = 20 ) )

  mtext( 'DIC (umol kg-1)', side = 4, line = 3, col = 'red' )

# Begin the second figure, showing the saturation state (left axis) and pCO2
#   (right axis):

  plot( seq( from = 1, to = 12 ), csys$OmegaCalcite, type = 'l', col = 'blue',
    lwd = 5, xlab = '', ylab = 'Calcite saturation state', col.lab = 'blue',
    xaxt = 'n', ylim = c( 0, 3 ) )

  axis( 1, at = seq( from = 1, to = 12, by = 1 ), # Customized labels on horizontal axis.
    labels = c('Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec') )

# Add pCO2 on the right axis of the same figure, colored in red. For reference,
#   the average atmospheric CO2 over the period 1985-2023 (i.e. the period the
#   atlas is representative of) was ~375 uatm. pCO2 values above/below that
#   value are associated with outgassing and ingassing, respectively.

  par( new = TRUE )

  plot( seq( from = 1, to = 12 ), csys$pCO2, type = 'l', col = 'red',
    lwd = 5, axes = FALSE, bty = 'n', xlab = '', ylab = '', ylim = c(0, 600) )

  axis( 4, at = seq( from = 0, to = 600, by = 50 ) )

  mtext( 'pCO2 (uatm)', side = 4, line = 3, col = 'red' )

# Now that we have the two figures created, we can have a bit of fun with
#   seacarb. Let's ask the question: how much of the seasonality that we see
#   in the saturation state and in pCO2 is due to the *direct* influence of
#   temperature on these two variables? To answer this, replace the temperature
#   timeseries by its annual average value, and repeat the same steps as above.

  timeseries_temp <- rep( mean( timeseries_temp ), times = 12 )

  csys <- carb( flag    = 15,                     # We are providing TA and DIC (see doc).
                var1    = timeseries_ta  * 1.e-6, # TA,  converted from umol/kg to mol/kg.
                var2    = timeseries_dic * 1.e-6, # DIC, converted from umol/kg to mol/kg.
                S       = timeseries_salt,        # Seawater practical salinity.
                T       = timeseries_temp,        # Seawater potential temperature (deg.C).
                Patm    = 1.,                     # Atmospheric pressure (default = 1atm).
                P       = 0.6,                    # Bottom of WE4.1 = 6meters = 0.6bar.
                Pt      = 0.,                     # Total phosphate concentration (mol/kg).
                Sit     = 0.,                     # Total silicate  concentration (mol/kg).
                k1k2    = 'cw',                   # Diss.constants from Cai & Wang 1998.
                kf      = 'x',                    # Default value for Kf.
                ks      = 'd',                    # Default value for Ks (Dickson 1990).
                pHscale = 'T',                    # Default value is Total Scale.
                b       = 'u74',                  # Default value for boron concentration.
                gas     = 'potential',            # (Irrelevant unless pCO2 is an input.)
                warn    = 'y',                    # Yes, show me warnings (if any).
                eos     = 'eos80',                # When providing potential T and practical S.
                long    = lo_s,                   # Longitude (irrelevant when using eos80).
                lat     = la_s )                  # Latitude  (irrelevant when using eos80).

  plot( seq( from = 1, to = 12 ), csys$OmegaCalcite, type = 'l', col = 'blue',
    lwd = 5, xlab = '', ylab = 'Calcite saturation state with constant T', col.lab = 'blue',
    xaxt = 'n', ylim = c( 0, 3 ) )

  axis( 1, at = seq( from = 1, to = 12, by = 1 ), # Customized labels on horizontal axis.
    labels = c('Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec') )

# Add pCO2 on the right axis of the same figure, colored in red:

  par( new = TRUE )

  plot( seq( from = 1, to = 12 ), csys$pCO2, type = 'l', col = 'red',
    lwd = 5, axes = FALSE, bty = 'n', xlab = '', ylab = '', ylim = c(0, 600) )

  axis( 4, at = seq( from = 0, to = 600, by = 50 ) )

  mtext( 'pCO2 (uatm) with constant T', side = 4, line = 3, col = 'red' )