Merge pull request #30 from matt-frey/ape-calculation

APE calculation

configure.ac

@@ -155,4 +155,36 @@ else
     AC_MSG_RESULT([no])
 fi
 
+
+ENABLE_IW_TEST_CASE='no'
+AC_ARG_ENABLE([iw-test-case],
+              [AS_HELP_STRING([--enable-iw-test-case], [enable IW test case (default=no)])],
+              [ENABLE_IW_TEST_CASE=$enableval])
+
+AM_CONDITIONAL([ENABLE_IW_TEST_CASE], [test "$ENABLE_IW_TEST_CASE" = "yes"])
+
+AC_MSG_CHECKING([whether we are enabling the IW test case])
+if test "x$ENABLE_IW_TEST_CASE" = "xyes"; then
+    AC_MSG_RESULT([yes])
+    FCFLAGS="$FCFLAGS -DENABLE_IW_TEST_CASE"
+else
+    AC_MSG_RESULT([no])
+fi
+
+
+ENABLE_RT_TEST_CASE='no'
+AC_ARG_ENABLE([rt-test-case],
+              [AS_HELP_STRING([--enable-rt-test-case], [enable RT test case (default=no)])],
+              [ENABLE_RT_TEST_CASE=$enableval])
+
+AM_CONDITIONAL([ENABLE_RT_TEST_CASE], [test "$ENABLE_RT_TEST_CASE" = "yes"])
+
+AC_MSG_CHECKING([whether we are enabling the RT test case])
+if test "x$ENABLE_RT_TEST_CASE" = "xyes"; then
+    AC_MSG_RESULT([yes])
+    FCFLAGS="$FCFLAGS -DENABLE_RT_TEST_CASE"
+else
+    AC_MSG_RESULT([no])
+fi
+
 AC_OUTPUT

src/advance.f90

@@ -204,7 +204,7 @@ subroutine adapt(t)
 #ifdef ENABLE_BUOYANCY
             double precision             :: yp(0:nz, 0:ny-1, 0:nx-1)        ! derivatives in y physical space
             double precision             :: zp(0:nz, 0:ny-1, 0:nx-1)        ! derivatives in z physical space
-            double precision             :: pe
+            double precision             :: ape
 #endif
             double precision             :: strain(3, 3), eigs(3)
             double precision             :: dudx(0:nz, 0:ny-1, 0:nx-1)      ! du/dx in physical space
@@ -280,8 +280,8 @@ subroutine adapt(t)
             en = get_enstrophy()
 #ifdef ENABLE_BUOYANCY
             call field_combine_physical(sbuoy, buoy)
-            pe = get_potential_energy()
-            write(WRITE_ECOMP, '(1x,f13.6,3(1x,1p,e14.7))') t , ke, pe, en
+            ape = get_available_potential_energy()
+            write(WRITE_ECOMP, '(1x,f13.6,3(1x,1p,e14.7))') t , ke, ape, en
 #else
             write(WRITE_ECOMP, '(1x,f13.6,2(1x,1p,e14.7))') t , ke, en
 #endif

src/fields/field_diagnostics_netcdf.f90

@@ -21,7 +21,7 @@ module field_diagnostics_netcdf
     integer            :: ncid
     integer            :: t_axis_id, t_dim_id, n_writes, ke_id, en_id
 #ifdef ENABLE_BUOYANCY
-    integer            :: pe_id, bmax_id, bmin_id
+    integer            :: ape_id, bmax_id, bmin_id
 #endif
 
     double precision   :: restart_time
@@ -98,13 +98,13 @@ subroutine create_netcdf_field_stats_file(basename, overwrite)
 #ifdef ENABLE_BUOYANCY
             call define_netcdf_dataset(                                     &
                 ncid=ncid,                                                  &
-                name='pe',                                                  &
-                long_name='domain-averaged potential energy',               &
+                name='ape',                                                 &
+                long_name='domain-averaged available potential energy',     &
                 std_name='',                                                &
                 unit='m^2/s^2',                                             &
                 dtype=NF90_DOUBLE,                                          &
                 dimids=(/t_dim_id/),                                        &
-                varid=pe_id)
+                varid=ape_id)
 
             call define_netcdf_dataset(                                     &
                 ncid=ncid,                                                  &
@@ -142,7 +142,7 @@ subroutine read_netcdf_field_stats_content
             call get_var_id(ncid, 'en', en_id)
 
 #ifdef ENABLE_BUOYANCY
-            call get_var_id(ncid, 'pe', pe_id)
+            call get_var_id(ncid, 'ape', ape_id)
 
             call get_var_id(ncid, 'min_buoyancy', bmin_id)
 
@@ -157,7 +157,7 @@ subroutine write_netcdf_field_stats(t)
             double precision, intent(in) :: t
             double precision             :: ke, en
 #ifdef ENABLE_BUOYANCY
-            double precision             :: bmin, bmax, pe
+            double precision             :: bmin, bmax, ape
 #endif
 
             call start_timer(field_stats_io_timer)
@@ -172,7 +172,7 @@ subroutine write_netcdf_field_stats(t)
 
 #ifdef ENABLE_BUOYANCY
             call field_combine_physical(sbuoy, buoy)
-            pe = get_potential_energy()
+            ape = get_available_potential_energy()
             bmin = minval(buoy)
             bmax = maxval(buoy)
 #endif
@@ -188,7 +188,7 @@ subroutine write_netcdf_field_stats(t)
             call write_netcdf_scalar(ncid, ke_id, ke, n_writes)
             call write_netcdf_scalar(ncid, en_id, en, n_writes)
 #ifdef ENABLE_BUOYANCY
-            call write_netcdf_scalar(ncid, pe_id, pe, n_writes)
+            call write_netcdf_scalar(ncid, ape_id, ape, n_writes)
             call write_netcdf_scalar(ncid, bmin_id, bmin, n_writes)
             call write_netcdf_scalar(ncid, bmax_id, bmax, n_writes)
 #endif

src/fields/fields.f90

@@ -3,12 +3,13 @@
 !     and functions.
 ! =============================================================================
 module fields
-    use parameters, only : nx, ny, nz, vcell, ncelli, dx, lower, extent, ngrid, vdomaini
+    use parameters, only : nx, ny, nz, ncelli, dx, lower, extent
     use constants, only : zero, f12, f14, one
     use merge_sort
     use inversion_utils, only : fftxys2p, diffx, diffy, central_diffz   &
                               , field_combine_semi_spectral             &
                               , field_decompose_semi_spectral
+    use ape_density, only : ape_den
     implicit none
 
     ! x: zonal
@@ -39,8 +40,6 @@ module fields
     ! initial \xi and \eta mean
     double precision :: ini_vor_mean(2)
 
-    double precision :: peref
-
     contains
 
         ! Allocate all fields
@@ -88,62 +87,31 @@ subroutine field_default
             ini_vor_mean = zero
         end subroutine field_default
 
-        subroutine calculate_peref
-#ifdef ENABLE_BUOYANCY
-            integer          :: ii(ngrid), i, j, k, n, m
-            double precision :: b(ngrid)
-            double precision :: gam, zmean
-
-            n = 1
-            m = nz+1
-            do i = 0, nx-1
-                do j = 0, ny-1
-                    b(n:m) = buoy(:, j, i)
-                    n = m + 1
-                    m = n + nz
-                enddo
-            enddo
-
-            call msort(b, ii)
-
-            gam = one / (extent(1) * extent(2))
-            zmean = f12 * gam * vcell
-
-            peref = - b(1) * vcell * zmean
-
-            do k = 2, ngrid
-                zmean = zmean + gam * vcell
-                peref = peref - b(k) * vcell * zmean
-            enddo
-
-            ! divide by domain volume to get domain-averaged peref
-            peref = peref * vdomaini
-#endif
-        end subroutine calculate_peref
-
-        ! domain-averaged potential energy
-        function get_potential_energy() result(pe)
-            double precision :: pe
+        ! domain-averaged available potential energy
+        function get_available_potential_energy() result(ape)
+            double precision :: ape
 #ifdef ENABLE_BUOYANCY
             integer          :: i, j, k
             double precision :: z(0:nz)
 
             do k = 0, nz
-                z(k) = dble(k) * dx(3)
+                z(k) = lower(3) + dble(k) * dx(3)
             enddo
 
-            pe = zero
+            ape = zero
             do i = 0, nx-1
                 do j = 0, ny-1
-                    pe = pe - sum(buoy(:, j, i) * z(:))
+                    ape = ape + sum(ape_den(buoy(1:nz-1, j, i), z(1:nz-1))) &
+                        + f12 *     ape_den(buoy(0,      j, i), z(0))       &
+                        + f12 *     ape_den(buoy(nz,     j, i), z(nz))
                 enddo
             enddo
 
-            pe = pe * ncelli - peref
+            ape = ape * ncelli
 #else
-            pe = zero
+            ape = zero
 #endif
-        end function get_potential_energy
+        end function get_available_potential_energy
 
         ! domain-averaged kinetic energy
         function get_kinetic_energy() result(ke)

src/ps3d.f90

@@ -45,7 +45,7 @@ subroutine pre_run
                               , output              &
                               , read_config_file    &
                               , time
-            double precision :: bbdif, ke, pe, te, en
+            double precision :: bbdif, ke, ape, te, en
 
             call register_timer('ps', ps_timer)
             call register_timer('field I/O', field_io_timer)
@@ -82,17 +82,15 @@ subroutine pre_run
             ! calculate the initial \xi and \eta mean and save it in ini_vor_mean:
             ini_vor_mean = calc_vorticity_mean()
 
-            call calculate_peref
-
             call vor2vel
 #ifdef ENABLE_BUOYANCY
             bbdif = maxval(buoy) - minval(buoy)
 #else
             bbdif = zero
 #endif
             ke = get_kinetic_energy()
-            pe = get_potential_energy()
-            te = ke + pe
+            ape = get_available_potential_energy()
+            te = ke + ape
             en = get_enstrophy()
 
 #ifdef ENABLE_BUOYANCY

src/utils/Makefile.am

@@ -10,5 +10,6 @@ libps_utils_la_SOURCES = 	\
 	config.f90		\
 	physics.f90		\
 	jacobi.f90		\
-	merge_sort.f90
+	merge_sort.f90		\
+	ape_density.f90
 

src/utils/ape_density.f90

@@ -0,0 +1,29 @@
+module ape_density
+    use constants
+    implicit none
+
+    public :: ape_den
+
+    contains
+
+        elemental function ape_den(b, z) result(a)
+            double precision, intent(in) :: b       ! buoyancy value
+            double precision, intent(in) :: z       ! height
+            double precision             :: a       ! APE density
+            double precision             :: br
+
+#ifdef ENABLE_IW_TEST_CASE
+            a = f18 * (b - four * z) ** 2
+#elif ENABLE_RT_TEST_CASE
+            br = max(b, -one)
+            br = min(br, one)
+            a = br * dasin(br) + dsqrt(one - br ** 2) - z * br - dcos(z)
+#else
+            ! dummy line to avoid compiler warning of 'unused variables'
+            a = b + z
+
+            a = zero
+#endif
+        end function ape_den
+
+end module ape_density