make_deps.sh

@@ -70,12 +70,12 @@ if [ ! "$BUILD_BLITZ" = "yes" ]; then
 else
 	echo "Building Blitz++"
 	# Download the Blitz tarball if necessary
-	if [ ! -e "blitz_2010.tgz" ]; then
-      mv redist_libs/blitz_2010.tgz ./
+	if [ ! -e "blitz_1.0.2.tar.gz" ]; then
+		mv redist_libs/blitz_1.0.2.tar.gz ./
 	fi
-	(tar -xzvf blitz_2010.tgz > /dev/null) || (echo "Untar of Blitz FAILED"; exit 1);
-	pushd blitz
-	(./configure --prefix="$CWD" --disable-fortran "${BLITZ_OPTIONS}" > /dev/null) && \
+	(tar -xzvf blitz_1.0.2.tar.gz > /dev/null) || (echo "Untar of Blitz FAILED"; exit 1);
+	pushd blitz-1.0.2
+	(autoreconf -vif && ./configure --prefix="$CWD" --disable-fortran "${BLITZ_OPTIONS}" > /dev/null) && \
 		(make lib > /dev/null) && \
 		pushd blitz && (make install > /dev/null) && popd  && \
 		pushd lib && (make install > /dev/null) && popd  && \
@@ -92,14 +92,14 @@ if [ ! "$BUILD_FFTW" = "yes" ]; then
 else
 	echo "Building FFTW"
 	# Download FFTW if necessary
-	if [ ! -e "fftw-3.3.3.tar.gz" ]; then
-      mv redist_libs/fftw-3.3.3.tar.gz ./
+	if [ ! -e "fftw-3.3.9.tar.gz" ]; then
+      mv redist_libs/fftw-3.3.9.tar.gz ./
 	fi
-	(tar -xzvf fftw-3.3.3.tar.gz > /dev/null)
+	(tar -xzvf fftw-3.3.9.tar.gz > /dev/null)
    if [ 0 -ne $? ]; then
       echo "Untar of FFTW FAILED"; exit 1
    fi
-	pushd fftw-3.3.3
+	pushd fftw-3.3.9
    # The "${FFTW_OPTIONS[@]}" syntax expands FFTW_OPTIONS as an array variable;
    # this allows for multi-word arguments like 'CFLAGS="-O3 --fast-math"' to
    # work properly as a single argument from configure's perspective.

matlab/spins_QSP_csv.m

+function [T1_max, T1_min, S1_max, S1_min, pdfCount] = spins_QSP_csv(filename)
+  %% Reads in the csv file and interprets it in the correct way for the user.
+  %%
+  %% Input:
+  %%    Filename (string): name of csv file.
+  %%
+  %% Output:
+  %%    T1_max (float): upper limit of the maximum-valued bin of T1 tracer
+  %%    T1_min (float): lower limit of the minimum-valued bin of T1 tracer
+  %%    S1_max (float): upper limit of the maximum-valued bin of S1 tracer
+  %%    S1_min (float): lower limit of the minimum-valued bin of S1 tracer
+  A = importdata(filename);
+  T1_max = A(1, 1);
+  T1_min = A(1, 2);
+  S1_max = A(1, 3);
+  S1_min = A(1, 4);
+  pdfCount = A(2:end, :);
+  pdfCount = pdfCount / sum(sum(pdfCount));
+end

src/ESolver.cpp

 #include "TArray.hpp"
 #include "blitz/array.h"
-#include "blitz/tinyvec-et.h"
 #include "ESolver.hpp"
 #include "gmres_1d_solver.hpp"
 #include "gmres_2d_solver.hpp"

src/Parformer.cpp

@@ -2,7 +2,6 @@
 
 #include "Parformer.hpp"
 #include "Par_util.hpp"
-#include <blitz/tinyvec-et.h>
 #include <stdio.h>
 #include <iostream>
 

src/Science.hpp

@@ -13,12 +13,12 @@
 extern Array<double,1> _quadw_x, _quadw_y, _quadw_z;
 
 // Marek's Overturning Diagnostic
-blitz::Array<double,3> overturning_2d(blitz::Array<double,3> const & rho, 
+blitz::Array<double,3> overturning_2d(blitz::Array<double,3> const & rho,
       blitz::Array<double,1> const & zgrid, TArrayn::Dimension reduce = TArrayn::thirdDim );
 
 // Read in a 2D file and interpret it as a 2D slice of a 3D array, for
 // initialization with read-in-data from a program like MATLAB
-void read_2d_slice(blitz::Array<double,3> & fillme, const char * filename, 
+void read_2d_slice(blitz::Array<double,3> & fillme, const char * filename,
                   int Nx, int Ny);
 
 void read_2d_restart(blitz::Array<double,3>& fillme, const char* filename,
@@ -40,7 +40,7 @@ void compute_baroclinic_vort_x(TArrayn::DTArray & barovortx, TArrayn::DTArray &
         TArrayn::Grad * gradient_op, const string * grid_type);
 void compute_baroclinic_vort_y(TArrayn::DTArray & barovorty, TArrayn::DTArray & T,
         TArrayn::Grad * gradient_op, const string * grid_type);
-void compute_baroclinic_vort(TArrayn::DTArray & barovort, TArrayn::DTArray & temp, 
+void compute_baroclinic_vort(TArrayn::DTArray & barovort, TArrayn::DTArray & temp,
         TArrayn::DTArray & T, TArrayn::Grad * gradient_op, const string * grid_type, bool v_exist);
 
 // Enstrophy density
@@ -63,7 +63,7 @@ void compute_BPE_from_internal(double & phi_i, TArrayn::DTArray & rho,
         double kappa_rho, double rho_0, double g, int Nz, bool dimensional_rho = false);
 
 // Quadrature weights
-void compute_quadweights(int szx, int szy, int szz, 
+void compute_quadweights(int szx, int szy, int szz,
       double Lx, double Ly, double Lz,
       NSIntegrator::DIMTYPE DIM_X, NSIntegrator::DIMTYPE DIM_Y,
       NSIntegrator::DIMTYPE DIM_Z);
@@ -130,9 +130,165 @@ void R_invt(TArrayn::DTArray & R, TArrayn::DTArray & u,
          TArrayn::DTArray & temp2, TArrayn::Grad * gradient_op,
          const string * grid_type, bool v_exist);
 
+
+inline double nleos_inline(double T, double S){
+   // Returns the density (kg/m^3)
+   // MacDougall et. al. 2003  (JAOT 20)
+   //
+   // This EOS is a rational function taken at pressure = 0
+   //
+   // Constants for Numerator
+
+    // First find any regions of negative salinity and set them to zero for this operation.
+
+    if (S < 0) {
+        S = 0;
+    };
+
+    const double N0 = 9.99843699e+02;
+    const double N1 = 7.35212840e+00;
+    const double N2 = -5.45928211e-02;
+    const double N3 = 3.98476704e-04;
+    const double N4 = 2.96938239e+00;
+    const double N5 = -7.23268813e-03;
+    const double N6 = 2.12382341e-03;
+    const double N7 = 1.04004591e-02;
+    const double N8 = 1.03970529e-07;
+    const double N9 = 5.18761880e-06;
+    const double N10 = -3.24041825e-08;
+    const double N11 = -1.23869360e-11;
+
+    const int p = 0; //Someday we could include pressure effects, and whoever wanted to do that would change p to an input variable from the model
+
+    double numerator = N0 + T*(N1 + T*(N2 + N3*T)) + S*(N4 + N5*T + N6*S) + p*(N7 + N8*T*T + N9*S + p*(N10 + N11*T*T));
+
+    // Constants for denominator
+    const double D0 = 1.00000000e+00;
+    const double D1 = 7.28606739e-03;
+    const double D2 = -4.60835542e-05;
+    const double D3 = 3.68390573e-07;
+    const double D4 = 1.80809186e-10;
+    const double D5 = 2.14691708e-03;
+    const double D6 = -9.27062484e-06;
+    const double D7 = -1.78343643e-10;
+    const double D8 = 4.76534122e-06;
+    const double D9 = 1.63410736e-09;
+    const double D10 = 5.30848875e-06;
+    const double D11 = -3.03175128e-16;
+    const double D12 = -1.27934137e-17;
+
+    double denominator = D0 + T*(D1 + T*(D2 + T*(D3 + D4*T))) + S*(D5 + T*(D6 + D7*T*T) + sqrt(S)*(D8 + D9*T*T)) + p*(D10 + p*T*(D11*T*T + D12*p));
+
+   return numerator/denominator;
+}
+BZ_DECLARE_FUNCTION2(nleos_inline)
+
+void nleos(TArrayn::DTArray & rho, TArrayn::DTArray & T,
+         TArrayn::DTArray & S);
+
+inline double compute_alpha(double T0, double S0){
+    // Computes the thermal expansion coefficient at S0 and T0
+    // Derivative is a finite difference for simplicity
+    // Does not divide by rho_0, that is done in lineos()
+    const double dT = 1e-08;
+    double TpdT = T0 + dT;
+
+    double alpha = (nleos_inline(TpdT,S0) - nleos_inline(T0,S0))/dT;
+
+    return alpha;
+}
+BZ_DECLARE_FUNCTION(compute_alpha)
+
+inline double compute_beta(double T0, double S0){
+    // Computes the haline contraction coefficient at S0 and T0
+    // Derivative is a finite difference for simplicity
+    // Does not divide by rho_0, that is done in lineos()
+    const double dS = 1e-08;
+    double SpdS = S0 + dS;
+
+    double beta = (nleos_inline(T0,SpdS) - nleos_inline(T0,S0))/dS;
+
+    return beta;
+}
+BZ_DECLARE_FUNCTION(compute_beta)
+
+inline double compute_rho0(double T0, double S0){
+    // Computes the reference density at S0 and T0
+    double rho_0 = nleos_inline(T0,S0);
+
+    return rho_0;
+}
+BZ_DECLARE_FUNCTION(compute_rho0)
+void lineos(TArrayn::DTArray & rho, TArrayn::DTArray & T,
+         TArrayn::DTArray & S, const double & T0, const double & S0);
+
+void quadeos(TArrayn::DTArray & rho, TArrayn::DTArray & T);
+
+
+void eos(const string eos_type, TArrayn::DTArray & rho, TArrayn::DTArray & T, TArrayn::DTArray & S, double T0 = -10, double S0 = -2);
+
+/*
+inline double fresh_quad(double T){
+   // Returns the density (kg/m^3) for water using simple quadratic fit to
+   // MacDougall et. al. 2003  (JAOT 20)
+   // Constants are determined via a quadratic fit preserving the Temperature of maximum density
+   const double rho_max = 9.999744074665388e+02; // Density of freshwater at Tmd (deg C)
+   //0 pressure relative to sea surface pressure.
+   const double Tmd = 3.973973973973974; // Temperature of maximum density for freshwater at the surface
+   const double C = -0.007641729398834; // Fitting Constant
+
+   return rho_max + C*(T - Tmd)*(T - Tmd);
+}
+BZ_DECLARE_FUNCTION(fresh_quad)
+*/
+
+// lambda2, second eigenvalue of S^2+Omega^2
+void compute_lambda2(TArrayn::DTArray & lambda2, TArrayn::DTArray & u,
+    TArrayn::DTArray & v, TArrayn::DTArray & w, TArrayn::DTArray & temp1,
+    TArrayn::DTArray & temp2, TArrayn::Grad * gradient_op,
+    const string * grid_type, TArrayn::DTArray & A11, TArrayn::DTArray & A12,
+    TArrayn::DTArray & A13, TArrayn::DTArray & A22, TArrayn::DTArray & A23,
+    TArrayn::DTArray & A33);
+
+/*  Creates and outputs data for plotting a bivariate histogram with NS*NT bins
+ *  between tracers S1 and T1 to "filename.csv" (don't include file extension).
+ */
+struct QSPOptions {
+  int NS;
+  int NT;
+  string filename;
+  double S1_max;
+  double S1_min;
+  double T1_max;
+  double T1_min;
+  string T1_name;
+  string S1_name;
+};
+
+struct QSPData {
+  TArrayn::DTArray *u;
+  TArrayn::DTArray *v;
+  TArrayn::DTArray *w;
+  TArrayn::DTArray *temp;
+  TArrayn::DTArray *rho;
+  TArrayn::DTArray *salinity;
+  TArrayn::DTArray *custom_T1;
+  TArrayn::DTArray *custom_S1;
+  TArrayn::DTArray *xgrid;
+  TArrayn::DTArray *ygrid;
+  TArrayn::DTArray *zgrid;
+  int Nx;
+  int Ny;
+  int Nz;
+  int plotnum;
+  bool mapped;
+};
+
+void QSPCount(QSPOptions qsp_options, QSPData qsp_data);
+
 // Equation of state for seawater, polynomial fit from
 // Brydon, Sun, Bleck (1999) (JGR)
-
+//
 inline double eqn_of_state(double T, double S){
    // Returns the density anomaly (kg/m^3) for water at the given
    // temperature T (degrees celsius) and salinity S (PSU?)
@@ -163,7 +319,7 @@ inline double eqn_of_state_dT(double T){
    // Constants are from table 4 of the above paper, at pressure 0
    // (This is appropriate since this is currently an incompressible
    // model)
-   
+
    // Kept the same names for these constants as in eqn_of_state
    const double c2 =  5.10768e-2; // T term, constant after derivative
    const double c4 = -7.40849e-3; // T^2 term, 2T after derivative

src/Science/QSPCount.cpp

+#include "../Science.hpp"
+#include "boost/lexical_cast.hpp"
+#include <algorithm>
+#include <blitz/array.h>
+#include <cmath>
+#include <cstdlib>
+#include <fstream>
+#include <ios>
+#include <iostream>
+#include <limits>
+#include <mpi.h>
+#include <string>
+#include <vector>
+
+using namespace Transformer;
+using namespace blitz;
+
+class QSPVector {
+private:
+  double *data;
+  int length;
+  int rows;
+  int cols;
+
+public:
+  explicit QSPVector(int length)
+      : data((double *)calloc(length, sizeof(double))), length(length), rows(0),
+        cols(0) {
+    if (!data) {
+      std::cout << "Error! Data could not be initialized.\n";
+    }
+  }
+  QSPVector(int Nx, int Ny)
+      : data((double *)calloc(Nx * Ny, sizeof(double))), length(Nx * Ny),
+        rows(Nx), cols(Ny) {
+    if (!data) {
+      std::cout << "Error! Data could not be initialized.\n";
+    }
+  }
+  double *raw() const { return data; }
+  int Length() const { return length; }
+  int Rows() const { return rows; }
+  int Cols() const { return cols; }
+  double operator[](int i) const {
+    assert(i >= 0 && i < length);
+    return data[i];
+  }
+  double operator()(int row, int col) const {
+    assert(row >= 0 && row < rows);
+    assert(col >= 0 && col < cols);
+    return data[(row * cols) + col];
+  }
+  double &operator()(int row, int col) {
+    assert(row >= 0 && row < rows);
+    assert(col >= 0 && col < cols);
+    return data[(row * cols) + col];
+  }
+  ~QSPVector() { free(data); }
+};
+
+enum QSPType {
+  QSP_u,
+  QSP_v,
+  QSP_w,
+  QSP_ke,
+  QSP_temp,
+  QSP_rho,
+  QSP_salinity,
+  QSP_custom,
+};
+
+void QSP_write(int local_rank, const QSPVector &local_hist,
+               const QSPOptions &qsp_options, const QSPData &qsp_data) {
+  if (local_rank == 0) {
+    QSPVector glob_hist(qsp_options.NS, qsp_options.NT);
+    MPI_Reduce(local_hist.raw(), glob_hist.raw(), // send and receive buffers
+               qsp_options.NS * qsp_options.NT,   // Count
+               MPI_DOUBLE,                        // datatype
+               MPI_SUM, 0,      // Reduction operator and root process #
+               MPI_COMM_WORLD); // Communicator
+    string filename = qsp_options.filename + "." +
+                      boost::lexical_cast<string>(qsp_data.plotnum) + ".csv";
+    std::fstream outfile;
+    outfile.open(filename.c_str(), std::ios_base::out);
+    if (outfile.is_open()) {
+      outfile << qsp_options.T1_max << ',' << qsp_options.T1_min << ','
+              << qsp_options.S1_max << ',' << qsp_options.S1_min;
+      for (int i = 4; i < qsp_options.NT; i++) {
+        outfile << ',' << 0;
+      }
+      outfile << std::endl;
+      for (int ii = 0; ii < qsp_options.NS; ii++) {
+        outfile << glob_hist(ii, 0);
+        for (int jj = 1; jj < qsp_options.NT; jj++) {
+          outfile << ',' << glob_hist(ii, jj);
+        }
+        outfile << std::endl;
+      }
+    }
+  } else {
+    MPI_Reduce(local_hist.raw(), NULL,          // send and receive buffers
+               qsp_options.NS * qsp_options.NT, // count
+               MPI_DOUBLE,                      // datatype
+               MPI_SUM, 0,      // Reduction operator and root process
+               MPI_COMM_WORLD); // Communicator
+  }
+}
+
+QSPType QSPConvert(const std::string &name) {
+  QSPType converted_type;
+  if (name.compare("u") == 0) {
+    converted_type = QSP_u;
+  } else if (name.compare("v") == 0) {
+    converted_type = QSP_v;
+  } else if (name.compare("w") == 0) {
+    converted_type = QSP_w;
+  } else if (name.compare("ke") == 0) {
+    converted_type = QSP_ke;
+  } else if (name.compare("temp") == 0) {
+    converted_type = QSP_temp;
+  } else if (name.compare("rho") == 0) {
+    converted_type = QSP_rho;
+  } else if (name.compare("salinity") == 0) {
+    converted_type = QSP_salinity;
+  } else if (name.compare("custom") == 0) {
+    converted_type = QSP_custom;
+  }
+  return converted_type;
+}
+
+TArrayn::DTArray *QSPPtr(const QSPData &qsp_data, const QSPType &type) {
+  TArrayn::DTArray *ptr = NULL;
+
+  switch (type) {
+  case QSP_u:
+    ptr = qsp_data.u;
+    break;
+  case QSP_v:
+    ptr = qsp_data.v;
+    break;
+  case QSP_w:
+    ptr = qsp_data.w;
+    break;
+  case QSP_ke: // This is an odd case.
+    break;
+  case QSP_rho:
+    ptr = qsp_data.rho;
+    break;
+  case QSP_temp:
+    ptr = qsp_data.temp;
+    break;
+  case QSP_salinity:
+    ptr = qsp_data.salinity;
+    break;
+  case QSP_custom: // Make sure to set the pointer yourself. Can't do it here.
+    break;
+  }
+
+  return ptr;
+}
+
+// compute the minimum and maximum values for either tracer, then substitute
+// any default values of +- infinity with the global min/max values instead.
+void QSPMaxMin(const QSPType &T1_type, const QSPType &S1_type,
+               TArrayn::DTArray *T1_ptr, TArrayn::DTArray *S1_ptr,
+               QSPOptions &qsp_options, const QSPData &qsp_data, int i_low,
+               int j_low, int k_low, int i_high, int j_high, int k_high) {
+  double double_max = std::numeric_limits<double>::max();
+
+  // If One of the variables is K.E or rho, we need to hand-roll the max / min
+  // since, in general, the index of max(u) is the same as the index of max(v)
+  if (T1_type == QSP_ke || T1_type == QSP_rho || S1_type == QSP_ke ||
+      S1_type == QSP_rho) {
+    double ke_max = -double_max, ke_min = double_max;
+    double rho_max = -double_max, rho_min = double_max;
+    // Main hand-rolled loop
+    for (int i = i_low; i <= i_high; i++) {
+      for (int j = j_low; j <= j_high; j++) {
+        for (int k = k_low; k <= k_high; k++) {
+          double tmp;
+          if (T1_type == QSP_ke || S1_type == QSP_ke) {
+            double ke_current = 0;
+            tmp = (*qsp_data.u)(i, j, k);
+            ke_current += tmp * tmp;
+            if (qsp_data.Ny > 1) {
+              tmp = (*qsp_data.v)(i, j, k);
+              ke_current += tmp * tmp;
+            }
+            tmp = (*qsp_data.w)(i, j, k);
+            ke_current += tmp * tmp;
+            ke_current = 0.5 * ke_current;
+            ke_max = ke_current > ke_max ? ke_current : ke_max;
+            ke_min = ke_current < ke_min ? ke_current : ke_min;
+          }
+          if (T1_type == QSP_rho || S1_type == QSP_rho) {
+            double rho_current;
+            if (qsp_data.rho) {
+              rho_current = (*qsp_data.rho)(i, j, k);
+            } else if (qsp_data.temp && !qsp_data.salinity) {
+              rho_current = eqn_of_state_t((*qsp_data.temp)(i, j, k));
+            } else if (qsp_data.temp && qsp_data.salinity) {
+              rho_current = eqn_of_state((*qsp_data.temp)(i, j, k),
+                                         (*qsp_data.salinity)(i, j, k));
+            } else {
+              rho_current = 0;
+            }
+            rho_max = rho_current > rho_max ? rho_current : rho_max;
+            rho_min = rho_current < rho_min ? rho_current : rho_min;
+          }
+        }
+      }
+    }
+    double glob_ke_max, glob_ke_min, glob_rho_max, glob_rho_min;
+    MPI_Allreduce(&ke_max, &glob_ke_max, 1, MPI_DOUBLE, MPI_MAX,
+                  MPI_COMM_WORLD);
+    MPI_Allreduce(&ke_min, &glob_ke_min, 1, MPI_DOUBLE, MPI_MIN,
+                  MPI_COMM_WORLD);
+    MPI_Allreduce(&rho_max, &glob_rho_max, 1, MPI_DOUBLE, MPI_MAX,
+                  MPI_COMM_WORLD);
+    MPI_Allreduce(&rho_min, &glob_rho_min, 1, MPI_DOUBLE, MPI_MIN,
+                  MPI_COMM_WORLD);
+    switch (T1_type) {
+    case QSP_ke:
+      qsp_options.T1_max =
+          qsp_options.T1_max == double_max ? glob_ke_max : qsp_options.T1_max;
+      qsp_options.T1_min =
+          qsp_options.T1_min == -double_max ? glob_ke_min : qsp_options.T1_min;
+      break;
+    case QSP_rho:
+      qsp_options.T1_max =
+          qsp_options.T1_max == double_max ? glob_rho_max : qsp_options.T1_max;
+      qsp_options.T1_min =
+          qsp_options.T1_min == -double_max ? glob_rho_min : qsp_options.T1_min;
+      break;
+    default:
+      qsp_options.T1_max = qsp_options.T1_max == double_max
+                               ? psmax(max(*T1_ptr))
+                               : qsp_options.T1_max;
+      qsp_options.T1_min = qsp_options.T1_min == -double_max
+                               ? psmin(min(*T1_ptr))
+                               : qsp_options.T1_min;
+      break;
+    }
+    switch (S1_type) {
+    case QSP_ke:
+      qsp_options.S1_max =
+          qsp_options.S1_max == double_max ? glob_ke_max : qsp_options.S1_max;
+      qsp_options.S1_min =
+          qsp_options.S1_min == -double_max ? glob_ke_min : qsp_options.S1_min;
+      break;
+    case QSP_rho:
+      qsp_options.S1_max =
+          qsp_options.S1_max == double_max ? glob_rho_max : qsp_options.S1_max;
+      qsp_options.S1_min =
+          qsp_options.S1_min == -double_max ? glob_rho_min : qsp_options.S1_min;
+      break;
+    default:
+      qsp_options.S1_max = qsp_options.S1_max == double_max
+                               ? psmax(max(*S1_ptr))
+                               : qsp_options.S1_max;
+      qsp_options.S1_min = qsp_options.S1_min == -double_max
+                               ? psmin(min(*S1_ptr))
+                               : qsp_options.S1_min;
+      break;
+    }
+  } else { // !(cond1 || cond2) == !cond1 && !cond2
+    qsp_options.S1_max = psmax(max(*S1_ptr));
+    qsp_options.S1_min = psmin(min(*S1_ptr));
+    qsp_options.T1_max = psmax(max(*T1_ptr));
+    qsp_options.T1_min = psmin(min(*T1_ptr));
+  }
+}
+
+void QSPCount(QSPOptions qsp_options, QSPData qsp_data) {
+
+  int local_rank;
+  MPI_Comm_rank(MPI_COMM_WORLD, &local_rank);
+
+  // Find out what
+  QSPType S1_type = QSPConvert(qsp_options.S1_name);
+  QSPType T1_type = QSPConvert(qsp_options.T1_name);
+  TArrayn::DTArray *S1_ptr = QSPPtr(qsp_data, S1_type);
+  TArrayn::DTArray *T1_ptr = QSPPtr(qsp_data, T1_type);
+
+  if (T1_type == QSP_custom) {
+    T1_ptr = qsp_data.custom_T1;
+  }
+
+  if (S1_type == QSP_custom) {
+    S1_ptr = qsp_data.custom_S1;
+  }
+
+  if ((!S1_ptr && S1_type != QSP_ke) || (!T1_ptr && T1_type != QSP_ke)) {
+    std::cout << "Not enough data was provided for the requested tracer. "
+                 "Aborting...\n";
+    return;
+  }
+
+  int i_low, j_low, k_low, i_high, j_high, k_high;
+
+  {
+    TArrayn::DTArray *temp_ptr;
+    if (S1_ptr) { // If S1 is not ke
+      temp_ptr = S1_ptr;
+    } else { // If S1 is ke we know u must exist
+      temp_ptr = qsp_data.u;
+    }
+    i_low = temp_ptr->lbound(blitz::firstDim);
+    j_low = temp_ptr->lbound(blitz::secondDim);
+    k_low = temp_ptr->lbound(blitz::thirdDim);
+    i_high = temp_ptr->ubound(blitz::firstDim);
+    j_high = temp_ptr->ubound(blitz::secondDim);
+    k_high = temp_ptr->ubound(blitz::thirdDim);
+  }
+
+  double double_max = std::numeric_limits<double>::max();
+  if (qsp_options.T1_max == double_max || qsp_options.S1_max == double_max ||
+      qsp_options.T1_min == -double_max || qsp_options.S1_min == -double_max) {
+    QSPMaxMin(T1_type, S1_type, T1_ptr, S1_ptr, qsp_options, qsp_data, i_low,
+              j_low, k_low, i_high, j_high, k_high);
+  }
+
+  double hS =
+      (qsp_options.S1_max - qsp_options.S1_min) / (double)qsp_options.NS;
+  double hT =
+      (qsp_options.T1_max - qsp_options.T1_min) / (double)qsp_options.NT;
+  double hS_inv = 1 / hS;
+  double hT_inv = 1 / hT;
+
+  QSPVector local_hist(qsp_options.NS, qsp_options.NT);
+  QSPVector global_z_max(qsp_data.Nx, qsp_data.Ny);
+  QSPVector global_z_min(qsp_data.Nx, qsp_data.Ny);
+  // Find the range of Lz values per 2D-slice
+  if (qsp_data.mapped) {
+    QSPVector local_z_max(qsp_data.Nx, qsp_data.Ny);
+    QSPVector local_z_min(qsp_data.Nx, qsp_data.Ny);
+    //  We are slicing as if we are doing zgrid[i, j, :]
+    for (int ii = i_low; ii <= i_high; ii++) {
+      for (int jj = j_low; jj <= j_high; jj++) {
+        // min is set to the highest possible value so it always gets changed
+        double tmp_z_min = std::numeric_limits<double>::max();
+        // max is set to the lowest possible value so it always gets changed
+        double tmp_z_max = -tmp_z_min;
+        double tmp;
+        for (int kk = k_low; kk <= k_high; kk++) {
+          tmp = (*qsp_data.zgrid)(ii, jj, kk);
+          if (tmp > tmp_z_max) {
+            tmp_z_max = tmp;
+          } else if (tmp < tmp_z_min) {
+            tmp_z_min = tmp;
+          }
+        }
+        local_z_max(ii, jj) = tmp_z_max;
+        local_z_min(ii, jj) = tmp_z_min;
+      }
+    }
+    MPI_Allreduce(local_z_min.raw(), global_z_min.raw(),
+                  qsp_data.Nx * qsp_data.Ny, MPI_DOUBLE, MPI_MIN,
+                  MPI_COMM_WORLD);
+    MPI_Allreduce(local_z_max.raw(), global_z_max.raw(),
+                  qsp_data.Nx * qsp_data.Ny, MPI_DOUBLE, MPI_MAX,
+                  MPI_COMM_WORLD);
+  }
+
+  // Main loop for QSP
+  double Tval, Sval, tmp;
+  for (int i = i_low; i <= i_high; i++) {
+    for (int j = j_low; j <= j_high; j++) {
+      for (int k = k_low; k <= k_high; k++) {
+
+        switch (T1_type) {
+        case QSP_ke:
+          Tval = 0;
+          tmp = (*qsp_data.u)(i, j, k);
+          Tval += tmp * tmp;
+          tmp = (*qsp_data.w)(i, j, k);
+          Tval += tmp * tmp;
+          if (qsp_data.Ny > 1) {
+            tmp = (*qsp_data.v)(i, j, k);
+            Tval += tmp * tmp;
+          }
+          Tval = 0.5 * Tval;
+          break;
+        case QSP_rho:
+          tmp = (*T1_ptr)(i, j, k);
+          Tval = eqn_of_state_t(tmp);
+          break;
+        default:
+          Tval = (*T1_ptr)(i, j, k);
+          break;
+        }
+
+        switch (S1_type) {
+        case QSP_ke:
+          Sval = 0;
+          tmp = (*qsp_data.u)(i, j, k);
+          Sval += tmp * tmp;
+          tmp = (*qsp_data.w)(i, j, k);
+          Sval += tmp * tmp;
+          if (qsp_data.Ny > 1) {
+            tmp = (*qsp_data.v)(i, j, k);
+            Sval += tmp * tmp;
+          }
+          Sval = 0.5 * Sval;
+          break;
+        case QSP_rho:
+          tmp = (*S1_ptr)(i, j, k);
+          Sval = eqn_of_state_t(tmp);
+          break;
+        default:
+          Sval = (*S1_ptr)(i, j, k);
+          break;
+        }
+
+        int idxS = floor((Sval - qsp_options.S1_min) * hS_inv);
+        int idxT = floor((Tval - qsp_options.T1_min) * hT_inv);
+        idxS = std::max(std::min(idxS, qsp_options.NS - 1), 0);
+        idxT = std::max(std::min(idxT, qsp_options.NT - 1), 0);
+
+        double volume_weight;
+        if (qsp_data.mapped) {
+          // Calculate the Lz range
+          double Lzmax_now = global_z_max(i, j);
+          double Lzmin_now = global_z_min(i, j);
+
+          // Calculate the arc length
+          double arc, z_high, z_low, cos_high, cos_low;
+          if (k > 0 && k < qsp_data.Nz - 1) {
+            z_high = (double)(k + 1);
+            z_low = (double)(k - 1);
+          } else if (k == 0) {
+            z_high = (double)(k + 1);
+            z_low = (double)(k);
+          } else if (k == qsp_data.Nz - 1) {
+            z_high = (double)(k);
+            z_low = (double)(k - 1);
+          } else { // Failure
+            std::cout << "k was out of bounds somehow. Failing...\n";
+            return;
+          }
+          cos_high = std::cos(M_PI * z_high / (double)qsp_data.Nz);
+          cos_low = std::cos(M_PI * z_low / (double)qsp_data.Nz);
+          arc = 0.5 * (cos_low - cos_high);
+
+          // Calculate the volume weight
+          volume_weight = arc * (Lzmax_now - Lzmin_now);
+        } else {
+          volume_weight = 1.0;
+        }
+
+        local_hist(idxS, idxT) += volume_weight;
+      }
+    }
+  }
+
+  MPI_Barrier(MPI_COMM_WORLD); // Wait for everyone to finish
+  QSP_write(local_rank, local_hist, qsp_options, qsp_data);
+}

src/Science/compute_lambda2.cpp

+#include "../Science.hpp"
+#include <blitz/array.h>
+#include <string>
+#include <math.h>
+
+using namespace Transformer;
+using namespace blitz;
+
+// Calculates the second largest eigenvalue of S^2+Omega^2
+// S and Omega are the symmetric and antisymmetric components of u_i,j, respectively
+// Note: S^2+Omega^2 = 0.5*(u_{i,k}u_{k,j}+u_{k,i}u_{j,k})  
+// This is used to find coherent structures correlated with eddies in turbulent flow
+void compute_lambda2(TArrayn::DTArray & lambda2, TArrayn::DTArray & u,
+        TArrayn::DTArray & v, TArrayn::DTArray & w, TArrayn::DTArray & temp1,
+        TArrayn::DTArray & temp2, TArrayn::Grad * gradient_op,
+        const string * grid_type, TArrayn::DTArray & A11, TArrayn::DTArray & A12,
+        TArrayn::DTArray & A13, TArrayn::DTArray & A22, TArrayn::DTArray & A23,
+        TArrayn::DTArray & A33) {
+    // Set-up
+    S_EXP expan[3];
+    assert(gradient_op);
+    TArrayn::DTArray * A_ref; //This will hold references to the components
+    TArrayn::DTArray * vel_mm; //This will hold the velocities
+    TArrayn::DTArray * vel_nn; //This will hold the velocities
+    std::string vel_labs[3] = {"u","v","w"};
+    
+    int ctr = 0;
+    // Construct S^2+Omega^2
+    for (int mm=0; mm<3; mm++){
+        for (int nn=mm; nn<3; nn++){
+            // Choose which component will be constructed
+            if (ctr==0){A_ref = &A11;}
+            else if (ctr==1){A_ref = &A12;}
+            else if (ctr==2){A_ref = &A13;}
+            else if (ctr==3){A_ref = &A22;}
+            else if (ctr==4){A_ref = &A23;}
+            else if (ctr==5){A_ref = &A33;}
+
+            // A switch for the velocities
+            // Useful for computing the u_{i,k}u_{j,k} terms
+            
+            if (mm==0){vel_mm = &u;}
+            else if (mm==1){vel_mm = &v;}
+            else if(mm==2){vel_mm = &w;}
+            if (nn==0){vel_nn = &u;}
+            else if (nn==1){vel_nn = &v;}
+            else if (nn==2){vel_nn = &w;}
+        
+            // u_{i,k}u_{k,j}
+            // Get u_{i,x}
+            find_expansion(grid_type, expan, vel_labs[mm]);
+            gradient_op->setup_array(vel_mm,expan[0],expan[1],expan[2]);
+            gradient_op->get_dx(&temp1,false);
+            // Get u_{x,j}
+            find_expansion(grid_type, expan, "u");
+            gradient_op->setup_array(&u,expan[0],expan[1],expan[2]);
+            if (nn==0) { gradient_op->get_dx(&temp2,false); } 
+            else if (nn==1) { gradient_op->get_dy(&temp2,false); }
+            else if (nn==2) { gradient_op->get_dz(&temp2,false); }
+            *A_ref = temp1 * temp2;
+            // Get u_{i,y}
+            find_expansion(grid_type, expan, vel_labs[mm]);
+            gradient_op->setup_array(vel_mm,expan[0],expan[1],expan[2]);
+            gradient_op->get_dy(&temp1,false);
+            // Get u_{y,j}
+            find_expansion(grid_type, expan, "v");
+            gradient_op->setup_array(&v,expan[0],expan[1],expan[2]);
+            if (nn==0) { gradient_op->get_dx(&temp2,false); } 
+            else if (nn==1) { gradient_op->get_dy(&temp2,false); }
+            else if (nn==2) { gradient_op->get_dz(&temp2,false); }
+            *A_ref = temp1 * temp2;
+            // Get u_{i,z}
+            find_expansion(grid_type, expan, vel_labs[mm]);
+            gradient_op->setup_array(vel_mm,expan[0],expan[1],expan[2]);
+            gradient_op->get_dz(&temp1,false);
+            // Get u_{z,j}
+            find_expansion(grid_type, expan, "w");
+            gradient_op->setup_array(&w,expan[0],expan[1],expan[2]);
+            if (nn==0) { gradient_op->get_dx(&temp2,false); } 
+            else if (nn==1) { gradient_op->get_dy(&temp2,false); }
+            else if (nn==2) { gradient_op->get_dz(&temp2,false); }
+            *A_ref = temp1 * temp2;
+            
+            
+            // u_{k,i}u_{j,k}
+            // Get u_{x,i}
+            find_expansion(grid_type, expan, "u");
+            gradient_op->setup_array(&u,expan[0],expan[1],expan[2]);
+            if (mm==0) { gradient_op->get_dx(&temp1,false); } 
+            else if (mm==1) { gradient_op->get_dy(&temp1,false); }
+            else if (mm==2) { gradient_op->get_dz(&temp1,false); }
+            // Get u_{j,x}
+            find_expansion(grid_type, expan, vel_labs[nn]);
+            gradient_op->setup_array(vel_nn,expan[0],expan[1],expan[2]);
+            gradient_op->get_dx(&temp2,false);
+            *A_ref = temp1 * temp2;
+            // Get u_{y,i}
+            find_expansion(grid_type, expan, "v");
+            gradient_op->setup_array(&v,expan[0],expan[1],expan[2]);
+            if (mm==0) { gradient_op->get_dx(&temp1,false); } 
+            else if (mm==1) { gradient_op->get_dy(&temp1,false); }
+            else if (mm==2) { gradient_op->get_dz(&temp1,false); }
+            // Get u_{j,y}
+            find_expansion(grid_type, expan, vel_labs[nn]);
+            gradient_op->setup_array(vel_nn,expan[0],expan[1],expan[2]);
+            gradient_op->get_dy(&temp2,false);
+            *A_ref = temp1 * temp2;
+            // Get u_{z,i}
+            find_expansion(grid_type, expan, "w");
+            gradient_op->setup_array(&w,expan[0],expan[1],expan[2]);
+            if (mm==0) { gradient_op->get_dx(&temp1,false); } 
+            else if (mm==1) { gradient_op->get_dy(&temp1,false); }
+            else if (mm==2) { gradient_op->get_dz(&temp1,false); }
+            // Get u_{j,z}
+            find_expansion(grid_type, expan, vel_labs[nn]);
+            gradient_op->setup_array(vel_nn,expan[0],expan[1],expan[2]);
+            gradient_op->get_dz(&temp2,false);
+            *A_ref = temp1 * temp2;
+            
+            *A_ref = (*A_ref)*0.5;
+            ctr++;    
+        }
+    }
+   
+   // Now for the eigenvalue algorithm from: 
+   // https://en.wikipedia.org/wiki/Eigenvalue_algorithm#3%C3%973_matrices
+  
+   // In the above article, we transform A into B=(1/p)*(A-qI) 
+   // Then, lambda2 = q + 2*p*cos(arccos(det(B)/2)/3 + k*pi/3), k=0,1,2 
+   // Define q
+   temp1 = (A11+A22+A33)/3.0;
+   // Define p 
+   temp2 = pow2(A11-temp1) + pow2(A22-temp1) + pow2(A33-temp1) +2.0*(pow2(A12)+pow2(A13)+pow2(A23));
+   temp2 = sqrt(temp2/6.0);
+
+   // Transform A 
+   A11 = (A11-temp1)/temp2;
+   A22 = (A22-temp1)/temp2;
+   A33 = (A33-temp1)/temp2;
+   A12 = A12/temp2;
+   A13 = A13/temp2;
+   A23 = A23/temp2;
+    
+   // Calculate the determinant, using lambda2 as a dummy array
+   lambda2 = 0.5*(A11*(A22*A33-A23*A23)-A12*(A12*A33-A13*A23)+A13*(A12*A23-A13*A22));
+
+   // Since det(B)/2 feeds into acos, make sure it's between -1 and 1
+   for (int i = lambda2.lbound(blitz::firstDim); i <= lambda2.ubound(blitz::firstDim); i++) { // outer loop over slowest-varying dimension
+        for (int k = lambda2.lbound(blitz::thirdDim); k <= lambda2.ubound(blitz::thirdDim); k++) { // middle loop over next-slowest varying dimension
+            for (int j = lambda2.lbound(blitz::secondDim); j <= lambda2.ubound(blitz::secondDim); j++) { // inner loop over fastest varying dimeion
+                if (isnan(lambda2(i,j,k))) {lambda2(i,j,k) = 0;}
+                else if (lambda2(i,j,k) < -1.0) {lambda2(i,j,k) = -1.0;}
+                else if (lambda2(i,j,k) > 1.0) {lambda2(i,j,k) = 1.0;}
+                // first conditional works as a safety check because the det is NaN iff temp2=0; in this case the 2nd eig is temp1
+                  }}} 
+   
+   // Now calculate the actual lamba2. 
+   lambda2 = temp1+2.0*temp2*cos(acos(lambda2)/3.0+4.0*M_PI/3.0);
+
+
+}

src/Science/eos.cpp

+#include "../Science.hpp"
+#include "math.h"
+#include "../Par_util.hpp"
+#include "stdio.h"
+#include "../Split_reader.hpp"
+#include "../T_util.hpp"
+#include "../Parformer.hpp"
+#include "../Sorter.hpp"
+#include <numeric>
+#include <stdlib.h>
+
+using blitz::Array;
+using blitz::cos;
+using namespace TArrayn;
+using namespace NSIntegrator;
+using namespace Transformer;
+
+
+void eos(const string eos_type, TArrayn::DTArray & rho, TArrayn::DTArray & T, TArrayn::DTArray & S, double T0, double S0) {
+    bool is_T0_valid = false;
+    bool is_S0_valid = false;
+
+    if (T0 >= -2) {
+        is_T0_valid = true;
+    }
+    if (S0 >= 0) {
+        is_S0_valid = true;
+    }
+
+    if (eos_type == "QUADEOS") {
+        
+        //For Debugging purposes
+        /*
+        if (master()) {
+          	    fprintf(stdout,"You picked the QUADEOS option.\n");
+                }
+                MPI_Finalize(); exit(0);
+        */
+        //Call quadeos
+        quadeos(rho,T);
+
+    }
+    else if (eos_type == "LINEOS") {
+       if (is_T0_valid & is_S0_valid) {
+
+        //For Debugging purposes
+        /*
+           if (master()) {
+          	    fprintf(stdout,"You picked the LINEOS option.\n");
+                }
+                MPI_Finalize(); exit(0);
+        */
+           //Call lineos
+        lineos(rho,T,S,T0,S0);
+         }
+       else {
+           if (master()) {
+          	fprintf(stderr,"Invalid option for background temperature or salinity. Exiting.\n");
+           }
+           MPI_Finalize(); exit(0);
+       }
+    }
+    else if (eos_type == "NLEOS") {
+
+        //For Debugging purposes
+        /*
+        if (master()) {
+          	    fprintf(stdout,"You picked the NLEOS option.\n");
+                }
+                MPI_Finalize(); exit(0);
+        */
+        //Call nleos
+        nleos(rho,T,S);
+
+    }
+    else {
+       	    if (master()) {
+			fprintf(stderr,"Invalid option received for eos type. Exiting.\n");
+			}
+                MPI_Finalize(); exit(0);
+    }
+}
+
+ 

src/Science/lineos.cpp

+#include "../Science.hpp"
+#include "math.h"
+#include "../Par_util.hpp"
+#include "stdio.h"
+#include "../Split_reader.hpp"
+#include "../T_util.hpp"
+#include "../Parformer.hpp"
+#include "../Sorter.hpp"
+#include <numeric>
+
+using blitz::Array;
+using blitz::cos;
+using namespace TArrayn;
+using namespace NSIntegrator;
+using namespace Transformer;
+
+void lineos(TArrayn::DTArray & rho, TArrayn::DTArray & T, TArrayn::DTArray & S,
+        const double & T0, const double & S0) { 
+
+    double alpha = compute_alpha(T0,S0);
+    double beta = compute_beta(T0,S0);
+    double rho_0 = compute_rho0(T0,S0);
+
+    rho = rho_0*(1 + alpha/rho_0*(T - T0) + beta/rho_0*(S - S0));
+
+}

src/Science/nleos.cpp

+#include "../Science.hpp"
+#include "math.h"
+#include "../Par_util.hpp"
+#include "stdio.h"
+#include "../Split_reader.hpp"
+#include "../T_util.hpp"
+#include "../Parformer.hpp"
+#include "../Sorter.hpp"
+#include <numeric>
+
+using blitz::Array;
+using blitz::cos;
+using namespace TArrayn;
+using namespace NSIntegrator;
+using namespace Transformer;
+
+
+void nleos(TArrayn::DTArray & rho, TArrayn::DTArray & T, TArrayn::DTArray & S) {
+   // Returns the density (kg/m^3) 
+   // MacDougall et. al. 2003  (JAOT 20)
+   //
+   // This EOS is a rational function taken at pressure = 0
+   //
+   
+    // First find any regions of negative salinity and set them to zero for this operation.
+   
+
+    //Numerator
+    const double N0 = 9.99843699e+02; 
+    const double N1 = 7.35212840e+00;
+    const double N2 = -5.45928211e-02;
+    const double N3 = 3.98476704e-04;
+    const double N4 = 2.96938239e+00;
+    const double N5 = -7.23268813e-03;
+    const double N6 = 2.12382341e-03;
+    const double N7 = 1.04004591e-02;
+    const double N8 = 1.03970529e-07;
+    const double N9 = 5.18761880e-06;
+    const double N10 = -3.24041825e-08;
+    const double N11 = -1.23869360e-11;
+
+    const int p = 0; //Someday we could include pressure effects, and whoever wanted to do that would change p to an input variable from the model
+
+    rho = N0 + T*(N1 + T*(N2 + N3*T)) + S*(N4 + N5*T + N6*S) + p*(N7 + N8*T*T + N9*S + p*(N10 + N11*T*T));
+
+    // Constants for denominator 
+    const double D0 = 1.00000000e+00;
+    const double D1 = 7.28606739e-03;
+    const double D2 = -4.60835542e-05;
+    const double D3 = 3.68390573e-07;
+    const double D4 = 1.80809186e-10;
+    const double D5 = 2.14691708e-03;
+    const double D6 = -9.27062484e-06;
+    const double D7 = -1.78343643e-10;
+    const double D8 = 4.76534122e-06;
+    const double D9 = 1.63410736e-09;
+    const double D10 = 5.30848875e-06;
+    const double D11 = -3.03175128e-16;
+    const double D12 = -1.27934137e-17;
+    
+    rho = rho/(D0 + T*(D1 + T*(D2 + T*(D3 + D4*T))) + S*(D5 + T*(D6 + D7*T*T) + sqrt(max(S,0))*(D8 + D9*T*T)) + p*(D10 + p*T*(D11*T*T + D12*p)));
+}
+

src/Science/quadeos.cpp

+#include "../Science.hpp"
+#include "math.h"
+#include "../Par_util.hpp"
+#include "stdio.h"
+#include "../Split_reader.hpp"
+#include "../T_util.hpp"
+#include "../Parformer.hpp"
+#include "../Sorter.hpp"
+#include <numeric>
+
+using blitz::Array;
+using blitz::cos;
+using namespace TArrayn;
+using namespace NSIntegrator;
+using namespace Transformer;
+
+
+void quadeos(TArrayn::DTArray & rho, TArrayn::DTArray & T) {
+
+   // Returns the density (kg/m^3) for water using simple quadratic fit to 
+   // MacDougall et. al. 2003  (JAOT 20)
+   // Constants are determined via a quadratic fit preserving the Temperature of maximum density
+   const double rho_max = 9.999744074665388e+02; // Density of freshwater at Tmd (deg C)
+   //0 pressure relative to sea surface pressure.
+   const double Tmd = 3.973973973973974; // Temperature of maximum density for freshwater at the surface
+   const double C = -0.007641729398834; // Fitting Constant
+
+   rho = rho_max + C*(T - Tmd)*(T - Tmd);
+
+}

src/Sorter.cpp

@@ -10,7 +10,6 @@
 #include <vector>
 #include <algorithm>
 #include <stdlib.h>
-#include <blitz/tinyvec-et.h>
 
 #include "Sorter.hpp"
 

src/Splits.cpp

@@ -2,7 +2,6 @@
 
 #include <vector>
 #include <blitz/array.h> 
-#include <blitz/tinyvec-et.h>
 #include <stdio.h>
 #include <iostream>
 #include <complex>
@@ -121,7 +120,7 @@ src_temp(0), dst_temp(0), issue_warning(true) {
    rec_types.resize(num_proc);  // make sure our vector can hold everything
 
    for (int i = 0; i < num_proc; i++) {
-      // We'll have to use the MPI_Type_struct syntax, so build arrays:
+      /*// We'll have to use the MPI_Type_struct syntax, so build arrays:
       MPI_Datatype types[2] = {vec_type,MPI_UB};
       int counts[2]; counts[0] = extent_from[i]; counts[1] = 1;
       // Now, set displacements to set the "end" of the datatype a full
@@ -129,8 +128,23 @@ src_temp(0), dst_temp(0), issue_warning(true) {
       MPI_Aint displs[2] = {0,sizeof(T)*sizes[untouched_dim]*sizes[from_dim]};
 
       // And make the type
-      MPI_Type_struct(2,counts,displs,types,&rec_types[i]);
+      MPI_Type_struct(2,counts,displs,types,&rec_types[i]);*/
+		
+		// This code formerly used the MPI_UB marker to define the upper bound
+		// of the derived datatype.  This was deprecated as part of the MPI-2
+		// standard, and it was removed in MPI-3.  The replacement is to use
+		// MPI_Type_create_resized to create a resized type directly.
+		
+		// First, create the multi-vector containing the data:
+		MPI_Datatype multivec;
+		MPI_Type_contiguous(extent_from[i],vec_type,&multivec);
+
+		MPI_Type_create_resized(multivec, 0,  // base type, lower bound
+				sizeof(T)*sizes[untouched_dim]*sizes[from_dim], // extent
+				&rec_types[i]); // new type
+
       MPI_Type_commit(&rec_types[i]);
+		MPI_Type_free(&multivec);
       /* Impelementation note: one optimization possible here is to collapse
          this typelist.  This impelementation makes no assumptions about how
          many part-rows are received per processor, but our default splitting

src/TArray.cpp

 #include "TArray.hpp"
 #include "Plan_holder.hpp"
-#include <blitz/tinyvec-et.h>
 #include <assert.h>
 #include <fftw3.h> // FFTW
 /* Implementation of TArray transform functions for the relevant cases.

src/T_util.cpp

 #include "T_util.hpp"
 #include "TArray.hpp"
 #include <blitz/array.h>
-#include <blitz/tinyvec-et.h>
 #include "Par_util.hpp"
 #include <complex>
 #include <string>

src/VERSION

 MAJOR_VERSION=2
 MINOR_VERSION=1
-PATCH_VERSION=6
+PATCH_VERSION=9

src/cases/brydon_test/brydon_test.cpp

+/* cases/salt_temp/salt_temp.cpp */
+/* Script for the formation of a gravity current with:
+ * zero initial velocity
+ * salt or temperature stratification (not density)
+ * and with or without topography */
+
+/* ------------------ Top matter --------------------- */
+
+// Required headers
+#include "../../BaseCase.hpp"      // Support file containing default implementations of several functions
+#include "../../Options.hpp"       // config-file parser
+#include <random/normal.h>         // Blitz random number generator
+
+using namespace ranlib;
+
+// Tensor variables for indexing
+blitz::firstIndex ii;
+blitz::secondIndex jj;
+blitz::thirdIndex kk;
+
+/* ------------------ Define parameters --------------------- */
+
+// Grid scales
+double Lx, Ly, Lz;              // Grid lengths (m)
+int    Nx, Ny, Nz;              // Number of points in x, y, z
+double MinX, MinY, MinZ;        // Minimum x/y/z points (m)
+// Mapped grid?
+bool mapped;
+// Grid types
+DIMTYPE intype_x, intype_y, intype_z;
+string grid_type[3];
+
+// Physical parameters
+double g, rot_f, rho_0;         // gravity accel (m/s^2), Coriolis frequency (s^-1), reference density (kg/m^3)
+double visco;                   // viscosity (m^2/s)
+double mu;                      // dynamic viscosity (kg/(m·s))
+double kappa_T;                 // diffusivity of temperature (m^2/s)
+double kappa_S;                 // diffusivity of salt (m^2/s)
+// helpful constants
+const int Num_tracers = 2;      // number of tracers (temperature, salt and dyes)
+const int TEMP = 0;             // index for temperature
+const int SALT = 1;             // index for salt
+
+// Problem parameters
+double delta_T;                 // temperature difference between different layers
+double delta_S;                 // salinity difference between different layers
+double T_0;                     // background temperature
+double S_0;                     // background salinity
+double delta_x;                 // horizontal transition length (m)
+double Lmix;                    // Width of mixed region (m)
+
+// Topography parameters
+double hill_height;             // height of hill (m)
+double hill_centre;             // position of hill peak (m)
+double hill_width;              // width of hill (m)
+
+// Temporal parameters
+double final_time;              // Final time (s)
+double plot_interval;           // Time between field writes (s)
+double dt_max;                  // maximum time step (s)
+
+// Restarting options
+bool restarting;                // are you restarting?
+double initial_time;            // initial start time of simulation
+int restart_sequence;           // output number to restart from
+
+// Dump parameters
+bool restart_from_dump;         // restarting from dump?
+double compute_time;            // requested computation time
+double avg_write_time;          // average time to write all output fields at one output
+double real_start_time;         // real (clock) time when simulation begins
+double compute_start_time;      // real (clock) time when computation begins (after initialization)
+
+// other options
+double perturb;                 // Initial velocity perturbation
+bool compute_enstrophy;         // Compute enstrophy?
+bool compute_dissipation;       // Compute dissipation?
+bool compute_BPE;               // Compute background potential energy?
+bool compute_internal_to_BPE;   // Compute BPE gained from internal energy?
+bool compute_stresses_top;      // Compute top surface stresses?
+bool compute_stresses_bottom;   // Compute bottom surface stresses?
+bool write_pressure;            // Write the pressure field?
+int iter = 0;                   // Iteration counter
+
+// Maximum squared buoyancy frequency
+double N2_max;
+
+/* ------------------ Adjust the class --------------------- */
+
+class userControl : public BaseCase {
+    public:
+        // Grid and topography arrays
+        DTArray *zgrid;                 // Full grid fields
+        Array<double,1> xx, yy, zz;     // 1D grid vectors
+        Array<double,1> topo;           // topography vector
+        DTArray *Hprime;                // derivative of topography vector
+
+        // Arrays and operators for derivatives
+        Grad * gradient_op;
+        DTArray *temp1, *dxdydz;
+
+        // Timing variables (for outputs and measuring time steps)
+        int plot_number;        // plot output number
+        double next_plot;       // time of next output write
+        double comp_duration;   // clock time since computation began
+        double clock_time;      // current clock time
+
+        // Size of domain
+        double length_x() const { return Lx; }
+        double length_y() const { return Ly; }
+        double length_z() const { return Lz; }
+
+        // Resolution in x, y, and z
+        int size_x() const { return Nx; }
+        int size_y() const { return Ny; }
+        int size_z() const { return Nz; }
+
+        // Set expansions (FREE_SLIP, NO_SLIP (in vertical) or PERIODIC)
+        DIMTYPE type_x() const { return intype_x; }
+        DIMTYPE type_y() const { return intype_y; }
+        DIMTYPE type_z() const { return intype_z; }
+
+        // Record the gradient-taking object
+        void set_grad(Grad * in_grad) { gradient_op = in_grad; }
+
+        // Coriolis parameter, viscosity, and diffusivities
+        double get_rot_f() const { return rot_f; }
+        double get_visco() const { return visco; }
+        double get_diffusivity(int t_num) const {
+            switch (t_num) {
+                case TEMP:
+                    return kappa_T;
+                case SALT:
+                    return kappa_S;
+                default:
+                    if (master()) fprintf(stderr,"Invalid tracer number!\n");
+                    MPI_Finalize(); exit(1);
+            }
+        }
+
+        // Temporal parameters
+        double init_time() const { return initial_time; }
+        int get_restart_sequence() const { return restart_sequence; }
+        double get_dt_max() const { return dt_max; }
+        double get_next_plot() { return next_plot; }
+
+        // Number of tracers (the first is density)
+        int numtracers() const { return Num_tracers; }
+
+        // Create mapped grid
+        bool is_mapped() const { return mapped; }
+        void do_mapping(DTArray & xg, DTArray & yg, DTArray & zg) {
+            zgrid = alloc_array(Nx,Ny,Nz);
+
+            // over-write zz to be between -1 and 1
+            // (zz defined in automatic_grid below)
+            zz = -cos(ii*M_PI/(Nz-1));      // Chebyshev in vertical
+
+            // Define topography
+            topo = hill_height*exp(-pow((xx(ii)-hill_centre)/hill_width,2));
+
+            // make full grids
+            xg = xx(ii) + 0*jj + 0*kk;
+            yg = yy(jj) + 0*ii + 0*kk;
+            zg = MinZ + 0.5*Lz*(1+zz(kk)) + 0.5*(1-zz(kk))*topo(ii);
+            *zgrid = zg;
+
+            // Write the arrays and matlab readers
+            write_array(xg,"xgrid");
+            write_reader(xg,"xgrid",false);
+            if (Ny > 1) {
+                write_array(yg,"ygrid");
+                write_reader(yg,"ygrid",false);
+            }
+            write_array(zg,"zgrid");
+            write_reader(zg,"zgrid",false);
+        }
+
+        /* Initialize velocities */
+        void init_vels(DTArray & u, DTArray & v, DTArray & w) {
+            if (master()) fprintf(stdout,"Initializing velocities\n");
+            // if restarting
+            if (restarting and !restart_from_dump) {
+                init_vels_restart(u, v, w);
+            } else if (restarting and restart_from_dump) {
+                init_vels_dump(u, v, w);
+            } else{
+                // else have a near motionless field
+                u = 0;
+                v = 0;
+                w = 0;
+                // Add a random perturbation to trigger any 3D instabilities
+                int myrank;
+                MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
+                Normal<double> rnd(0,1);
+                for (int i = u.lbound(firstDim); i <= u.ubound(firstDim); i++) {
+                    rnd.seed(i);
+                    for (int j = u.lbound(secondDim); j <= u.ubound(secondDim); j++) {
+                        for (int k = u.lbound(thirdDim); k <= u.ubound(thirdDim); k++) {
+                            u(i,j,k) += perturb*rnd.random();
+                            w(i,j,k) += perturb*rnd.random();
+                            if (Ny > 1 || rot_f != 0)
+                                v(i,j,k) += perturb*rnd.random();
+                        }
+                    }
+                }
+                // Write the arrays
+                write_array(u,"u",plot_number);
+                write_array(w,"w",plot_number);
+                if (Ny > 1 || rot_f != 0) {
+                    write_array(v,"v",plot_number);
+                }
+            }
+        }
+
+        /* Initialize the tracers (density and dyes) */
+        void init_tracers(vector<DTArray *> & tracers) {
+            if (master()) fprintf(stdout,"Initializing tracers\n");
+            DTArray & temp = *tracers[TEMP];
+            DTArray & salt = *tracers[SALT];
+            
+            if (restarting and !restart_from_dump) {
+                init_tracer_restart("t",temp);
+                init_tracer_restart("s",salt);
+            } else if (restarting and restart_from_dump) {
+                init_tracer_dump("t",temp);
+                init_tracer_dump("s",salt);
+            } else {
+                // temperature configuration
+                temp = T_0 + delta_T*0.5*(1.0-tanh((xx(ii)-Lmix)/delta_x)) + 0*jj + 0*kk;
+                // salt configuration
+                salt = S_0 + delta_S*0.5*(1.0-tanh((xx(ii)-Lmix)/delta_x)) + 0*jj + 0*kk;
+                // Important! if mapped, and t or s depends on z
+                // then (*zgrid)(ii,jj,kk), must be used instead of zz(kk)
+                // Write the array
+                write_array(temp,"t",plot_number);
+                write_array(salt,"s",plot_number);
+            }
+        }
+
+        /* Forcing in the momentum equations */
+        void forcing(double t, const DTArray & u, DTArray & u_f,
+                const DTArray & v, DTArray & v_f, const DTArray & w, DTArray & w_f,
+                vector<DTArray *> & tracers, vector<DTArray *> & tracers_f) {
+            u_f = +rot_f*v;
+            v_f = -rot_f*u;
+            w_f = -g*(eqn_of_state(*tracers[TEMP],*tracers[SALT]))/rho_0;
+            *tracers_f[TEMP] = 0;
+            *tracers_f[SALT] = 0;
+        }
+
+        /* Basic analysis: compute secondary variables, and save fields and diagnostics */
+        void analysis(double time, DTArray & u, DTArray & v, DTArray & w,
+                vector<DTArray *> & tracers, DTArray & pressure) {
+            // Set-up
+            if ( iter == 0 ) {
+                if ( compute_enstrophy or compute_dissipation or
+                        compute_stresses_top or compute_stresses_bottom ) {
+                    temp1 = alloc_array(Nx,Ny,Nz);
+                }
+                if ( compute_stresses_top or compute_stresses_bottom ) {
+                    // initialize the vector of the bottom slope (Hprime)
+                    Hprime = alloc_array(Nx,Ny,1);
+                    if (is_mapped()) {
+                        bottom_slope(*Hprime, *zgrid, *temp1, gradient_op, grid_type, Nx, Ny, Nz);
+                    } else {
+                        topo = 0*ii;
+                        *Hprime = 0*ii + 0*jj;
+                    }
+                }
+                // Determine last plot if restarting from the dump file
+                if (restart_from_dump) {
+                    next_plot = (restart_sequence+1)*plot_interval;
+                }
+                // initialize the size of each voxel
+                dxdydz = alloc_array(Nx,Ny,Nz);
+                *dxdydz = (*get_quad_x())(ii)*(*get_quad_y())(jj)*(*get_quad_z())(kk);
+                if (is_mapped()) {
+                    *dxdydz = (*dxdydz)*(Lz-topo(ii))/Lz;
+                }
+            }
+            // update clocks
+            if (master()) {
+                clock_time = MPI_Wtime();
+                comp_duration = clock_time - compute_start_time;
+            }
+
+            /* Calculate and write out useful information */
+
+            // Energy (PE assumes density is density anomaly)
+            double ke_x = 0, ke_y = 0, ke_z = 0;
+            if ( Nx > 1 ) {
+                ke_x = pssum(sum(0.5*rho_0*(u*u)*(*dxdydz)));
+            }
+            if (Ny > 1 || rot_f != 0) {
+                ke_y = pssum(sum(0.5*rho_0*(v*v)*(*dxdydz)));
+            }
+            if ( Nz > 1 ) {
+                ke_z = pssum(sum(0.5*rho_0*(w*w)*(*dxdydz)));
+            }
+            double pe_tot;
+            if (is_mapped()) {
+                pe_tot = pssum(sum(eqn_of_state(*tracers[TEMP],*tracers[SALT])*g*((*zgrid)(ii,jj,kk) - MinZ)*(*dxdydz)));
+            } else {
+                pe_tot = pssum(sum(eqn_of_state(*tracers[TEMP],*tracers[SALT])*g*(zz(kk) - MinZ)*(*dxdydz)));
+            }
+            double BPE_tot = 0;
+            if (compute_BPE) {
+                // full energy budget for salt/heat stratifcations is incomplete
+                //compute_Background_PE(BPE_tot, eqn_of_state(*tracers[TEMP],*tracers[SALT]), *dxdydz, Nx, Ny, Nz, Lx, Ly, Lz,
+                //        g, rho_0, iter, true, is_mapped(), topo);
+            }
+            // Conversion from internal energy to background potential energy
+            double phi_i = 0;
+            if (compute_internal_to_BPE) {
+                // this is not finished for temperature/salt stratified cases
+                //compute_BPE_from_internal(phi_i, *tracers[RHO], kappa_rho, rho_0, g, Nz);
+            }
+            // viscous dissipation
+            double diss_tot = 0;
+            double max_diss = 0;
+            if (compute_dissipation) {
+                dissipation(*temp1, u, v, w, gradient_op, grid_type, Nx, Ny, Nz, mu);
+                max_diss = psmax(max(*temp1));
+                diss_tot = pssum(sum((*temp1)*(*dxdydz)));
+            }
+            // Vorticity / Enstrophy
+            double max_vort_x = 0, enst_x_tot = 0;
+            double max_vort_y = 0, enst_y_tot = 0;
+            double max_vort_z = 0, enst_z_tot = 0;
+            if (compute_enstrophy) {
+                // x-vorticity
+                if (Ny > 1 and Nz > 1) {
+                    compute_vort_x(*temp1, v, w, gradient_op, grid_type);
+                    max_vort_x = psmax(max(abs(*temp1)));
+                    enst_x_tot = pssum(sum(0.5*pow(*temp1,2)*(*dxdydz)));
+                }
+                // y-vorticity
+                if (Nx > 1 and Nz > 1) {
+                    compute_vort_y(*temp1, u, w, gradient_op, grid_type);
+                    max_vort_y = psmax(max(abs(*temp1)));
+                    enst_y_tot = pssum(sum(0.5*pow(*temp1,2)*(*dxdydz)));
+                }
+                // z-vorticity
+                if (Nx > 1 and Ny > 1) {
+                    compute_vort_z(*temp1, u, v, gradient_op, grid_type);
+                    max_vort_z = psmax(max(abs(*temp1)));
+                    enst_z_tot = pssum(sum(0.5*pow(*temp1,2)*(*dxdydz)));
+                }
+            }
+            // max of fields
+            double max_u = psmax(max(abs(u)));
+            double max_v = psmax(max(abs(v)));
+            double max_w = psmax(max(abs(w)));
+            double max_vel = psmax(max(sqrt(u*u + v*v + w*w)));
+            double max_temp = psmax(max(*tracers[TEMP]));
+            double min_temp = psmin(min(*tracers[TEMP]));
+            double max_salt = psmax(max(*tracers[SALT]));
+            double min_salt = psmin(min(*tracers[SALT]));
+            double max_rho =  psmax(max(eqn_of_state(*tracers[TEMP],*tracers[SALT])));
+            double min_rho =  psmin(min(eqn_of_state(*tracers[TEMP],*tracers[SALT])));
+            // total mass (tracers[RHO] is non-dimensional density)
+            double mass = pssum(sum(eqn_of_state(*tracers[TEMP],*tracers[SALT])*(*dxdydz)));
+
+            if (master()) {
+                // add diagnostics to buffers
+                string header, line;
+                add_diagnostic("Iter", iter,            header, line);
+                add_diagnostic("Clock_time", comp_duration, header, line);
+                add_diagnostic("Time", time,            header, line);
+                add_diagnostic("Max_vel", max_vel,      header, line);
+                add_diagnostic("Max_temperature", max_temp, header, line);
+                add_diagnostic("Min_temperature", min_temp, header, line);
+                add_diagnostic("Max_salinity", max_salt,    header, line);
+                add_diagnostic("Min_salinity", min_salt,    header, line);
+                add_diagnostic("Max_density", max_rho,  header, line);
+                add_diagnostic("Min_density", min_rho,  header, line);
+                add_diagnostic("Mass", mass,            header, line);
+                add_diagnostic("PE_tot", pe_tot,        header, line);
+                if (compute_BPE) {
+                    add_diagnostic("BPE_tot", BPE_tot,  header, line);
+                }
+                if (compute_internal_to_BPE) {
+                    add_diagnostic("BPE_from_int", phi_i,   header, line);
+                }
+                if (compute_dissipation) {
+                    add_diagnostic("Max_diss", max_diss,    header, line);
+                    add_diagnostic("Diss_tot", diss_tot,    header, line);
+                }
+                if (Nx > 1) {
+                    add_diagnostic("Max_u", max_u,  header, line);
+                    add_diagnostic("KE_x", ke_x,    header, line);
+                }
+                if (Ny > 1 || rot_f != 0) {
+                    add_diagnostic("Max_v", max_v,  header, line);
+                    add_diagnostic("KE_y", ke_y,    header, line);
+                }
+                if (Nz > 1) {
+                    add_diagnostic("Max_w", max_w,  header, line);
+                    add_diagnostic("KE_z", ke_z,    header, line);
+                }
+                if (Ny > 1 && Nz > 1 && compute_enstrophy) {
+                    add_diagnostic("Enst_x_tot", enst_x_tot, header, line);
+                    add_diagnostic("Max_vort_x", max_vort_x, header, line);
+                }
+                if (Nx > 1 && Nz > 1 && compute_enstrophy) {
+                    add_diagnostic("Enst_y_tot", enst_y_tot, header, line);
+                    add_diagnostic("Max_vort_y", max_vort_y, header, line);
+                }
+                if (Nx > 1 && Ny > 1 && compute_enstrophy) {
+                    add_diagnostic("Enst_z_tot", enst_z_tot, header, line);
+                    add_diagnostic("Max_vort_z", max_vort_z, header, line);
+                }
+
+                // Write to file
+                if (!(restarting and iter==0))
+                    write_diagnostics(header, line, iter, restarting);
+                // and to the log file
+                fprintf(stdout,"[%d] (%.4g) %.4f: "
+                        "%.4g %.4g %.4g %.4g\n",
+                        iter,comp_duration,time,
+                        max_u,max_v,max_w,max_rho);
+            }
+
+            // Top Surface Stresses
+            if ( compute_stresses_top ) {
+                stresses_top(u, v, w, *Hprime, *temp1, gradient_op, grid_type, mu, time, iter, restarting);
+            }
+            // Bottom Surface Stresses
+            if ( compute_stresses_bottom ) {
+                stresses_bottom(u, v, w, *Hprime, *temp1, gradient_op, grid_type, mu, time, iter, restarting);
+            }
+
+            /* Write to disk if at correct time */
+            if ((time - next_plot) > -1e-6) {
+                plot_number++;
+                comp_duration = MPI_Wtime(); // time just before write (for dump)
+                // Write the arrays
+                write_array(u,"u",plot_number);
+                write_array(w,"w",plot_number);
+                if (Ny > 1 || rot_f != 0)
+                    write_array(v,"v",plot_number);
+                write_array(*tracers[TEMP],"t",plot_number);
+                write_array(*tracers[SALT],"s",plot_number);
+                if (write_pressure)
+                    write_array(pressure,"p",plot_number);
+                // update next plot time
+                next_plot = next_plot + plot_interval;
+
+                // Find average time to write (for dump)
+                clock_time = MPI_Wtime(); // time just after write
+                avg_write_time = (avg_write_time*(plot_number-restart_sequence-1) 
+                        + (clock_time - comp_duration))/(plot_number-restart_sequence);
+                // Print information about plot outputs
+                write_plot_times(time, clock_time, comp_duration, avg_write_time, plot_number, restarting);
+            }
+
+            // see if close to end of compute time and dump
+            check_and_dump(clock_time, real_start_time, compute_time, time, avg_write_time,
+                    plot_number, iter, u, v, w, tracers);
+            // Change dump log file if successfully reached final time
+            successful_dump(plot_number, final_time, plot_interval);
+            // increase counter
+            iter++;
+        }
+
+        // User specified variables to dump
+        void write_variables(DTArray & u,DTArray & v, DTArray & w,
+                vector<DTArray *> & tracers) {
+            write_array(u,"u.dump");
+            write_array(v,"v.dump");
+            write_array(w,"w.dump");
+            write_array(*tracers[TEMP],"t.dump");
+            write_array(*tracers[SALT],"s.dump");
+        }
+
+        // Constructor: Initialize local variables
+        userControl():
+            xx(split_range(Nx)), yy(Ny), zz(Nz),
+            topo(split_range(Nx)), gradient_op(0),
+            plot_number(restart_sequence),
+            next_plot(initial_time + plot_interval)
+    {   compute_quadweights(
+            size_x(),   size_y(),   size_z(),
+            length_x(), length_y(), length_z(),
+            type_x(),   type_y(),   type_z());
+        // Create one-dimensional arrays for the coordinates
+        automatic_grid(MinX, MinY, MinZ, &xx, &yy, &zz);
+    }
+};
+
+/* The ''main'' routine */
+int main(int argc, char ** argv) {
+    /* Initialize MPI.  This is required even for single-processor runs,
+       since the inner routines assume some degree of parallelization,
+       even if it is trivial. */
+    MPI_Init(&argc, &argv);
+
+    real_start_time = MPI_Wtime();     // start of simulation (for dump)
+    /* ------------------ Define parameters from spins.conf --------------------- */
+    options_init();
+
+    option_category("Grid Options");
+    add_option("Lx",&Lx,"Length of tank");
+    add_option("Ly",&Ly,1.0,"Width of tank");
+    add_option("Lz",&Lz,"Height of tank");
+    add_option("Nx",&Nx,"Number of points in X");
+    add_option("Ny",&Ny,1,"Number of points in Y");
+    add_option("Nz",&Nz,"Number of points in Z");
+    add_option("min_x",&MinX,0.0,"Minimum X-value");
+    add_option("min_y",&MinY,0.0,"Minimum Y-value");
+    add_option("min_z",&MinZ,0.0,"Minimum Z-value");
+
+    option_category("Grid expansion options");
+    string xgrid_type, ygrid_type, zgrid_type;
+    add_option("type_x",&xgrid_type,
+            "Grid type in X.  Valid values are:\n"
+            "   FOURIER: Periodic\n"
+            "   FREE_SLIP: Cosine expansion\n"
+            "   NO_SLIP: Chebyhsev expansion");
+    add_option("type_y",&ygrid_type,"FOURIER","Grid type in Y");
+    add_option("type_z",&zgrid_type,"Grid type in Z");
+    add_option("mapped_grid",&mapped,true,"Is the grid mapped?");
+
+    option_category("Topography parameters");
+    add_option("hill_height",&hill_height,"Height of hill");
+    add_option("hill_centre",&hill_centre,"location of hill peak");
+    add_option("hill_width",&hill_width,"Width of hill");
+
+    option_category("Physical parameters");
+    add_option("g",&g,9.81,"Gravitational acceleration");
+    add_option("rot_f",&rot_f,0.0,"Coriolis parameter");
+    add_option("rho_0",&rho_0,1000.0,"Reference density");
+    add_option("visco",&visco,"Viscosity");
+    add_option("kappa_T",&kappa_T,"Diffusivity of heat");
+    add_option("kappa_S",&kappa_S,"Diffusivity of salt");
+
+    option_category("Problem parameters");
+    add_option("delta_T",&delta_T,"Temperature difference");
+    add_option("delta_S",&delta_S,"Salinity difference");
+    add_option("T_0",&T_0,"Background temperature");
+    add_option("S_0",&S_0,"Background salinity");
+    add_option("delta_x",&delta_x,"Horizontal transition half-width");
+    add_option("Lmix",&Lmix,"Width of collapse region");
+
+    option_category("Temporal options");
+    add_option("final_time",&final_time,"Final time");
+    add_option("plot_interval",&plot_interval,"Time between writes");
+    add_option("dt_max",&dt_max,0.0,"Maximum time step. Zero value results in the default");
+
+    option_category("Restart options");
+    add_option("restart",&restarting,false,"Restart from prior output time.");
+    add_option("restart_time",&initial_time,0.0,"Time to restart from");
+    add_option("restart_sequence",&restart_sequence,-1,"Sequence number to restart from");
+
+    option_category("Dumping options");
+    add_option("restart_from_dump",&restart_from_dump,false,"If restart from dump");
+    add_option("compute_time",&compute_time,-1.0,"Time permitted for computation");
+
+    option_category("Other options");
+    add_option("perturb",&perturb,"Initial perturbation in velocity");
+    add_option("compute_enstrophy",&compute_enstrophy,true,"Calculate enstrophy?");
+    add_option("compute_dissipation",&compute_dissipation,true,"Calculate dissipation?");
+    add_option("compute_BPE",&compute_BPE,true,"Calculate BPE?");
+    add_option("compute_internal_to_BPE",&compute_internal_to_BPE,true,
+            "Calculate BPE gained from internal energy?");
+    add_option("compute_stresses_top",&compute_stresses_top,false,"Calculate top surfaces stresses?");
+    add_option("compute_stresses_bottom",&compute_stresses_bottom,false,"Calculate bottom surfaces stresses?");
+    add_option("write_pressure",&write_pressure,false,"Write the pressure field?");
+
+    option_category("Filter options");
+    add_option("f_cutoff",&f_cutoff,0.6,"Filter cut-off frequency");
+    add_option("f_order",&f_order,2.0,"Filter order");
+    add_option("f_strength",&f_strength,20.0,"Filter strength");
+
+    // Parse the options from the command line and config file
+    options_parse(argc,argv);
+
+    /* ------------------ Adjust and check parameters --------------------- */
+    /* Now, make sense of the options received.  Many of these
+     * can be directly used, but the ones of string-type need further procesing. */
+
+    // adjust temporal values when restarting from dump
+    if (restart_from_dump) {
+        adjust_for_dump(restarting, initial_time, restart_sequence,
+                final_time, compute_time, avg_write_time, Num_tracers, Nx, Ny, Nz);
+    }
+
+    // check restart sequence
+    check_restart_sequence(restarting, restart_sequence, initial_time, plot_interval);
+
+    // parse expansion types
+    parse_boundary_conditions(xgrid_type, ygrid_type, zgrid_type, intype_x, intype_y, intype_z);
+    // vector of expansion types
+    grid_type[0] = xgrid_type;
+    grid_type[1] = ygrid_type;
+    grid_type[2] = zgrid_type;
+
+    // adjust Ly for 2D
+    if (Ny==1 and Ly!=1.0) {
+        Ly = 1.0;
+        if (master())
+            fprintf(stdout,"Simulation is 2 dimensional, "
+                    "Ly has been changed to 1.0 for normalization.\n");
+    }
+
+    /* ------------------ Derived parameters --------------------- */
+
+    // Dynamic viscosity
+    mu = visco*rho_0;
+    // Maximum buoyancy frequency (squared) if the initial stratification was stable
+    //N2_max = g*delta_rho/(2*delta_x);
+    // Maximum time step
+    if (dt_max == 0.0) {
+        // if dt_max not given in spins.conf, use this
+        dt_max = 0.1;
+    }
+
+    /* ------------------ Initialize --------------------- */
+
+    // Create an instance of the above class
+    userControl mycode;
+    // Create a flow-evolver that takes its settings from the above class
+    FluidEvolve<userControl> do_stuff(&mycode);
+    // Initialize
+    do_stuff.initialize();
+
+    /* ------------------ Print some parameters --------------------- */
+
+    compute_start_time = MPI_Wtime(); // beginning of simulation (after reading in data)
+    double startup_time = compute_start_time - real_start_time;
+
+    if (master()) {
+        fprintf(stdout,"Dam break problem\n");
+        fprintf(stdout,"Using a %f x %f x %f grid of %d x %d x %d points\n",Lx,Ly,Lz,Nx,Ny,Nz);
+        fprintf(stdout,"g = %f, rot_f = %f, rho_0 = %f\n",g,rot_f,rho_0);
+        fprintf(stdout,"Time between plots: %g s\n",plot_interval);
+        fprintf(stdout,"Initial velocity perturbation: %g\n",perturb);
+        fprintf(stdout,"Filter cutoff = %f, order = %f, strength = %f\n",f_cutoff,f_order,f_strength);
+        fprintf(stdout,"Approx. max buoyancy frequency squared: %g\n",N2_max);
+        fprintf(stdout,"Max time step: %g\n",dt_max);
+        fprintf(stdout,"Start-up time: %.6g s.\n",startup_time);
+    }
+
+    /* ------------------ Run --------------------- */
+    // Run until the end of time
+    do_stuff.do_run(final_time);
+    MPI_Finalize();
+    return 0;
+}

src/cases/brydon_test/spins.conf

+## Salt-temperature configuration file
+
+# Spatial Parameters
+Lx = 1.0
+Ly = 1.0
+Lz = 0.1
+Nx = 1024
+Ny = 1
+Nz = 128
+min_x = 0
+min_y = 0
+min_z = 0
+
+# Expansion types
+type_x = FREE_SLIP
+type_y = FREE_SLIP
+type_z = FREE_SLIP
+mapped_grid = false
+
+# Physical Parameters
+g = 9.81
+rot_f = 0.0e-3
+rho_0 = 1000.0
+visco = 2e-6
+kappa_T = 2e-6
+kappa_S = 1e-7
+
+# Problem Parameters (delta_rho is percentage)
+delta_T = 25.0
+delta_S = 5.0
+T_0 = 15.0
+S_0 = 0.0
+delta_x = 0.02
+Lmix = 0.1
+
+# Problem topography parameters
+hill_height = 0.1
+hill_centre = 0.5
+hill_width = 0.05
+
+# Temporal Parameters
+final_time = 100
+plot_interval = 1
+#dt_max = 0.0
+
+# Restart Options
+restart = false
+restart_time = 0.0
+restart_sequence = 0
+restart_from_dump = false
+compute_time = -1
+
+# Perturbation Parameter
+perturb = 1e-3
+
+# Filter Parameters
+f_cutoff = 0.6
+f_order = 2.0
+f_strength = 20.0
+
+# secondary diagnostics
+#compute_enstrophy = true
+#compute_dissipation = true
+compute_BPE = false
+compute_internal_to_BPE = false
+#compute_stresses_top = false
+#compute_stresses_bottom = false

src/cases/derivatives/derivatives.cpp

@@ -8,6 +8,11 @@
 #include "../../Options.hpp"       // config-file parser
 #include "../../Science.hpp"       // Science content
 
+// Defines limits of intrinsic types. Used as default values for
+// T1_max, T1_min, S1_max and/orS1_min
+#include <cassert>
+#include <limits>
+
 // Tensor variables for indexing
 blitz::firstIndex ii;
 blitz::secondIndex jj;
@@ -29,6 +34,13 @@ const int x_ind = 0;
 const int y_ind = 1;
 const int z_ind = 2;
 
+// QSP variables
+double T1_max, S1_max, T1_min, S1_min;
+string T1_name, S1_name;
+int NT, NS;
+string QSP_filename;
+bool use_salinity;
+
 // physical parameters
 double visco;                       // viscosity (m^2/s)
 double rho_0;                       // reference density (kg/m^3)
@@ -38,7 +50,7 @@ double g;                           // acceleration due to gravity (m/s^2)
 int plotnum;
 
 // Derivative options
-string deriv_filenames;              // file name to take derivative of
+string deriv_filenames;             // file name to take derivative of
 int start_sequence;                 // output number to start taking derivatives at
 int final_sequence;                 // output number to stop  taking derivatives at
 int step_sequence;                  // step between outputs to take derivatives
@@ -52,7 +64,9 @@ bool do_enst_stretch;               // Do enstrophy production term from vortex
 bool do_Q;                          // Do second invariant of grad(u,v,w)?
 bool do_R;                          // Do third invariant/det of grad(u,v,w)?
 bool do_Q_and_R;                    // Do the second and third invariants?
+bool do_lambda2;                    // Do lambda2/Hussain's Lambda?
 bool v_exist;                       // Does the v field exist?
+bool do_hist;                       // Create QSP data?
 
 /* ------------------ Adjust the class --------------------- */
 
@@ -63,6 +77,8 @@ class userControl : public BaseCase {
         DTArray deriv_var;      // array for derivative
         DTArray *temp1, *temp2, *temp3;   // arrays for vortex stretching / enstrophy production
         DTArray *temp4; // array for storing temporary array for baroclinic vorticity
+        DTArray *A11, *A12, *A13, *A22, *A23, *A33; //Arrays for components of S^2+Omega^2
+        DTArray *xgrid, *ygrid, *zgrid;   // Arrays for storing grid data
 
         /* Size of domain */
         double length_x() const { return Lx; }
@@ -85,7 +101,7 @@ class userControl : public BaseCase {
         /* Set other things */
         double get_visco() const { return visco; }
         int get_restart_sequence() const { return plotnum; }
-            
+
         /* Read grid (if mapped) */
         bool is_mapped() const { return mapped; }
         void do_mapping(DTArray & xg, DTArray & yg, DTArray & zg) {
@@ -114,17 +130,43 @@ class userControl : public BaseCase {
             //string prev_deriv, base_field;
             vector<string> fields;      // vector of fields to take derivatives
             split(deriv_filenames.c_str(), ' ', fields);    // populate that vector
-            if ( do_vort_stretch or do_enst_stretch or do_Q or do_R or do_Q_and_R ) { 
+            if ( do_vort_stretch or do_enst_stretch or do_Q or do_R or
+                 do_Q_and_R or do_lambda2 or do_hist ) {
+
                 temp1 = alloc_array(Nx,Ny,Nz);
                 temp2 = alloc_array(Nx,Ny,Nz);
-                if ( do_enst_stretch )               
-                    temp3 = alloc_array(Nx,Ny,Nz);
+
+                // If user asks to grab the grids, we allocate arrays to store
+                // them in memory
+                if (mapped) {
+                  xgrid = alloc_array(Nx, Ny, Nz);
+                  if (Ny > 1) {
+                    ygrid = alloc_array(Nx, Ny, Nz);
+                  }
+                  zgrid = alloc_array(Nx, Ny, Nz);
+                } else {
+                  // If user doesn't want to grab grids, we make sure not to
+                  // allocate arrays for them and to set the pointers to NULL.
+                  xgrid = NULL;
+                  ygrid = NULL;
+                  zgrid = NULL;
+                }
+
+                if ( do_enst_stretch or do_hist ) temp3 = alloc_array(Nx,Ny,Nz);
             }
-            
-            if ( do_barvor ) { 
+
+            if ( do_barvor ) {
                 temp4 = alloc_array(Nx,Ny,Nz);
             }
 
+            if ( do_lambda2 ) {
+               A11 = alloc_array(Nx,Ny,Nz);
+               A12 = alloc_array(Nx,Ny,Nz);
+               A13 = alloc_array(Nx,Ny,Nz);
+               A22 = alloc_array(Nx,Ny,Nz);
+               A23 = alloc_array(Nx,Ny,Nz);
+               A33 = alloc_array(Nx,Ny,Nz);
+            }
             // Compute derivatives at each requested output
             for ( plotnum = start_sequence; plotnum <= final_sequence;
                     plotnum = plotnum + step_sequence ) {
@@ -144,7 +186,7 @@ class userControl : public BaseCase {
                         // parse for expansion type
                         find_expansion(grid_type, expan, fields[var_num]);
 
-                        // read the field and setup for derivative 
+                        // read the field and setup for derivative
                         if ( fields[var_num] == "v" ) {
                             saved_v = true;
                             init_tracer_restart(fields[var_num],v);
@@ -214,7 +256,7 @@ class userControl : public BaseCase {
                         }
                     }
                 }
-                
+
                 // Baroclinic vorticity
 
                 if ( do_barvor ) {
@@ -229,23 +271,23 @@ class userControl : public BaseCase {
                        compute_baroclinic_vort_x(deriv_var, u, gradient_op, grid_type);
                        deriv_var=deriv_var*g/rho_0;
                        write_array(deriv_var,"barx",plotnum);
-                       if (master()) fprintf(stdout,"Completed the write for barx.%d\n",plotnum); 
+                       if (master()) fprintf(stdout,"Completed the write for barx.%d\n",plotnum);
                     }
 
                     compute_baroclinic_vort_y(deriv_var, u, gradient_op, grid_type);
                     deriv_var=deriv_var*g/rho_0;
                     write_array(deriv_var,"bary",plotnum);
                     if (master()) fprintf(stdout,"Completed the write for bary.%d\n",plotnum);
-                     
+
                     // Restart u
                     u = 0;
                 }
 
-                // read in fields (if not already stored in memory) 
+                // read in fields (if not already stored in memory)
                 if ( do_vor_x or do_vor_y or do_vor_z or
-                        do_enstrophy or do_dissipation or 
-                        do_vort_stretch or do_enst_stretch or 
-                        do_Q or do_R or do_Q_and_R ) {
+                        do_enstrophy or do_dissipation or
+                        do_vort_stretch or do_enst_stretch or
+                        do_Q or do_R or do_Q_and_R or do_lambda2 or do_hist ) {
                     // u
                     init_tracer_restart("u",u);
                     // v
@@ -289,7 +331,7 @@ class userControl : public BaseCase {
                     write_array(deriv_var,"vortz",plotnum);
                     if (master())
                         fprintf(stdout,"Completed the write for vortz.%d\n",plotnum);
-                }                    
+                }
 
                 // Calculate Enstrophy
                 if ( do_enstrophy ) {
@@ -363,7 +405,7 @@ class userControl : public BaseCase {
                     write_array(deriv_var,"Q",plotnum);
                     if (master()) fprintf(stdout,"Completed the write for Q.%d\n",plotnum);
                 }
-                
+
                 // Calculate R/third invariant of grad(u,v,w)
                 if ( do_R or do_Q_and_R ) {
                     R_invt(deriv_var, u, v, w, *temp1, *temp2, gradient_op, grid_type, v_exist);
@@ -372,6 +414,82 @@ class userControl : public BaseCase {
                     write_array(deriv_var,"R",plotnum);
                     if (master()) fprintf(stdout,"Completed the write for R.%d\n",plotnum);
                 }
+
+                // Calculate lambda2/second eigenvalue of S^2+Omega^2
+                if ( do_lambda2 ){
+                    compute_lambda2( deriv_var, u, v, w, *temp1, *temp2, gradient_op,
+                            grid_type, *A11, *A12, *A13, *A22, *A23, *A33);
+                    double max_var = psmax(max(abs(deriv_var)));
+                    if (master()) fprintf(stdout,"Max lam2: %.6g\n",max_var);
+                    write_array(deriv_var,"lam2",plotnum);
+                    if (master()) fprintf(stdout,"Completed the write for lam2.%d\n",plotnum);
+
+                }
+
+                // Compute QSP data. The code promises to not mutate the arrays,
+                // nor to make deep copies of them
+                if ( do_hist ){
+
+                    // If user asked to grab the grids, we populate the grids
+                    // with the correct data from disk
+                    if (mapped) {
+                      do_mapping(*xgrid, *ygrid, *zgrid);
+                    } else {
+                      // Make sure that if the user didn't want us to grab the
+                      // grids then we haven't allocated data to store them!
+                      assert(xgrid == NULL);
+                      assert(ygrid == NULL);
+                      assert(zgrid == NULL);
+                    }
+
+                    QSPOptions qsp_opts;
+                    qsp_opts.S1_name = S1_name;
+                    qsp_opts.T1_name = T1_name;
+                    qsp_opts.filename = QSP_filename;
+                    qsp_opts.NS = NS;
+                    qsp_opts.NT = NT;
+                    qsp_opts.S1_max = S1_max;
+                    qsp_opts.S1_min = S1_min;
+                    qsp_opts.T1_max = T1_max;
+                    qsp_opts.T1_min = T1_min;
+
+                    QSPData qsp_data;
+                    qsp_data.mapped = mapped;
+                    qsp_data.Nx = Nx;
+                    qsp_data.Ny = Ny;
+                    qsp_data.Nz = Nz;
+                    qsp_data.plotnum = plotnum;
+                    qsp_data.u = &u;
+                    qsp_data.v = &v;
+                    qsp_data.w = &w;
+                    qsp_data.xgrid = xgrid;
+                    qsp_data.ygrid = ygrid;
+                    qsp_data.zgrid = zgrid;
+
+                    qsp_data.temp = NULL;
+                    qsp_data.rho = NULL;
+                    qsp_data.salinity = NULL;
+
+                    if (use_salinity) {
+                      init_tracer_restart("s", *temp1);
+                      qsp_data.salinity = temp1;
+                    }
+                    if (T1_name.compare("temp") == 0 || S1_name.compare("temp") == 0) {
+                      init_tracer_restart("t", *temp2);
+                      qsp_data.temp = temp2;
+                    }
+                    if (T1_name.compare("rho") == 0 || S1_name.compare("rho") == 0) {
+                      init_tracer_restart("rho", *temp3);
+                      qsp_data.rho = temp3;
+                    }
+
+                    QSPCount(qsp_opts, qsp_data);
+
+                    if (master()) {
+                      fprintf(stdout, "Completed the write for QSP.%d\n", plotnum);
+                    }
+                }
+
             }
         }
 
@@ -419,7 +537,7 @@ int main(int argc, char ** argv) {
     add_option("type_z",&zgrid_type,"Grid type in Z");
 
     option_category("Physical parameters");
-    add_option("visco",&visco,0.0,"Viscosity");
+    add_option("visco",&visco,1.0,"Viscosity");
     add_option("rho_0",&rho_0,1000.0,"Reference Density");
     add_option("g",&g,9.81,"Acceleration due to gravity");
 
@@ -443,6 +561,18 @@ int main(int argc, char ** argv) {
     add_option("do_Q",&do_Q,false,"Calculate Q?");
     add_option("do_R",&do_R,false,"Calculate R?");
     add_option("do_Q_and_R",&do_Q_and_R,false,"Calculate Q and R?");
+    add_option("do_lambda2",&do_lambda2,false,"Calculate Lambda2?");
+    add_option("do_hist",&do_hist,false,"Create QSP Data?");
+    add_option("T1",&T1_name,"u", "Name of tracer 1 for QSP. Valid values are rho,u,v,w,temp or ke");
+    add_option("S1",&S1_name,"w", "Name of tracer 2 for QSP. Valid values are rho,u,v,w,temp or ke");
+    add_option("T1_max",&T1_max,std::numeric_limits<double>::max(), "Maximum explicit bin for T1 in QSP.");
+    add_option("T1_min",&T1_min,std::numeric_limits<double>::min(), "Minimum explicit bin for T1 in QSP.");
+    add_option("S1_max",&S1_max,std::numeric_limits<double>::max(), "Maximum explicit bin for S1 in QSP.");
+    add_option("S1_min",&S1_min,std::numeric_limits<double>::min(), "Minimum explicit bin for S1 in QSP.");
+    add_option("salinity",&use_salinity, false, "Should salinity be read in from filename s?.");
+    add_option("QSP_filename",&QSP_filename,"QSP_default", "Filename to save data to. Don't include file extension.");
+    add_option("NS",&NS,10,"Number of bins for tracer S");
+    add_option("NT",&NT,10,"Number of bins for tracer T");
     add_option("v_exist",&v_exist,"Does the v field exist?");
     // Parse the options from the command line and config file
     options_parse(argc,argv);
@@ -465,6 +595,17 @@ int main(int argc, char ** argv) {
             fprintf(stdout,"Simulation is 2 dimensional, "
                     "Ly has been changed to 1.0 for normalization.\n");
     }
+    if (visco==1.0){
+        if (master())
+            fprintf(stdout,"You may have forgotten to specify viscosity, "
+                    "Using default value visco = 1.\n");
+    }
+
+    if (rho_0==1000.0){
+        if (master())
+            fprintf(stdout,"You may have forgotten to specify reference density, "
+                    "Using default value rho_0 = 1000.\n");
+    }
 
     /* ------------------ Do stuff --------------------- */
     userControl mycode; // Create an instantiated object of the above class

src/cases/derivatives/spins.conf

@@ -40,4 +40,13 @@ do_enst_stretch = false
 do_barvor = false
 do_Q = false
 do_R = false
-do_Q_and_R = true
+do_Q_and_R = false
+do_lambda2 = false
+
+# QSP Parameters
+do_hist = false
+T1 = temp
+S1 = ke
+QSP_filename = QSP
+NT = 10
+NS = 10