Add computation of energy conversion from internal to BPE

This is the last energy conversion term to get a complete
energy budget.

src/Science.cpp

@@ -454,6 +454,46 @@ void compute_Background_PE(double & BPE_tot, TArrayn::DTArray & rho,
     delete[] local_vols_r;
 }
 
+// Internal energy converted to Background potential energy
+// See Winters et al. 1995 (Available potential energy and mixing in density-stratified fluids)
+// phi_i = - kappa * g * (surface integral of density at top boundary
+//                        minus the surface integral of density at bottom boundary)
+void compute_BPE_from_internal(double & phi_i, TArrayn::DTArray & rho,
+        double kappa_rho, double rho_0, double g, int Nz,
+        bool dimensional_rho, bool mapped, TArrayn::DTArray * Hprime) {
+    // Tensor variables for indexing
+    blitz::firstIndex ii;
+    blitz::secondIndex jj;
+    blitz::Range all = blitz::Range::all();
+
+    if ( !mapped ) {
+        // Unmapped case
+        if ( !dimensional_rho ) {
+            phi_i = -kappa_rho * g * pssum(sum(
+                        rho_0 * (rho(all,all,Nz-1) - rho(all,all,0))
+                        * (*get_quad_x())(ii) * (*get_quad_y())(jj) ));
+        } else {
+            phi_i = -kappa_rho * g * pssum(sum(
+                        (rho(all,all,Nz-1) - rho(all,all,0))
+                        * (*get_quad_x())(ii) * (*get_quad_y())(jj) ));
+        }
+    } else {
+        fprintf(stderr,"BPE from internal energy calculation is not available in mapped cases.\n");
+        assert(0);
+        // Mapped case
+        // Need to fix
+        /*if ( !dimensional_rho ) {
+            phi_i = -kappa_rho * g * pssum(sum(
+                        rho_0 * ( rho(all,all,Nz-1) - rho(all,all,0)*pow(1+pow(*Hprime,2),0.5) )
+                        * (*get_quad_x())(ii) * (*get_quad_y())(jj) ));
+        } else {
+            phi_i = -kappa_rho * g * pssum(sum(
+                        ( rho(all,all,Nz-1) - rho(all,all,0)*pow(1+pow(*Hprime,2),0.5) )
+                        * (*get_quad_x())(ii) * (*get_quad_y())(jj) ));
+        }*/
+    }
+}
+
 // Global arrays to store quadrature weights
 Array<double,1> _quadw_x, _quadw_y, _quadw_z;
 

src/Science.hpp

@@ -39,10 +39,14 @@ void dissipation(TArrayn::DTArray & diss, TArrayn::DTArray & u, TArrayn::DTArray
         TArrayn::DTArray & w, TArrayn::Grad * gradient_op, const string * grid_type,
         const int Nx, const int Ny, const int Nz, const double visco);
 
-// Background Potential Energy
+// Background Potential Energy (BPE)
 void compute_Background_PE(double & BPE_tot, TArrayn::DTArray & rho, int Nx, int Ny, int Nz,
         double Lx, double Ly, double Lz, double g, double rho_0, int iter,
         bool mapped = false, Array<double,1> hill = Array<double,1>(), bool dimensional_rho = false);
+// Internal energy converted to BPE
+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, bool mapped = false, TArrayn::DTArray * Hprime = NULL);
 
 // Quadrature weights
 void compute_quadweights(int szx, int szy, int szz, 

src/cases/gravity_current/gravity_current.cpp

@@ -234,6 +234,9 @@ class userControl : public BaseCase {
                 diss_tot = pssum(sum((*temp1)*
                             (*get_quad_x())(ii)*(*get_quad_y())(jj)*(*get_quad_z())(kk)));
             }
+            // Conversion from internal energy to background potential energy
+            double phi_i;
+            compute_BPE_from_internal(phi_i, *tracers[RHO], kappa_rho, rho_0, g, Nz);
             // Vorticity / Enstrophy
             double max_vort_x, enst_x_tot;
             double max_vort_y, enst_y_tot;
@@ -285,6 +288,7 @@ class userControl : public BaseCase {
                 add_diagnostic("KE_x", ke_x,            header, line);
                 add_diagnostic("KE_z", ke_z,            header, line);
                 add_diagnostic("PE_tot", pe_tot,        header, line);
+                add_diagnostic("BPE_from_int", phi_i,   header, line);
                 if (compute_BPE) {
                     add_diagnostic("BPE_tot", BPE_tot,  header, line);
                 }