From d6a30520fb091a15eac8718216e7098874ec93c4 Mon Sep 17 00:00:00 2001 From: Matthias Frey Date: Wed, 1 Feb 2023 15:26:39 +0000 Subject: [PATCH 1/3] APE calculation using APE density --- configure.ac | 32 ++++++++++++++ src/advance.f90 | 6 +-- src/fields/field_diagnostics_netcdf.f90 | 16 +++---- src/fields/fields.f90 | 58 ++++++------------------- src/ps3d.f90 | 8 ++-- src/utils/Makefile.am | 3 +- src/utils/ape_density.f90 | 26 +++++++++++ 7 files changed, 87 insertions(+), 62 deletions(-) create mode 100644 src/utils/ape_density.f90 diff --git a/configure.ac b/configure.ac index 762058d..53bc6b6 100644 --- a/configure.ac +++ b/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 diff --git a/src/advance.f90 b/src/advance.f90 index 6de864a..1713792 100644 --- a/src/advance.f90 +++ b/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 diff --git a/src/fields/field_diagnostics_netcdf.f90 b/src/fields/field_diagnostics_netcdf.f90 index c4951cd..04f5c9c 100644 --- a/src/fields/field_diagnostics_netcdf.f90 +++ b/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 diff --git a/src/fields/fields.f90 b/src/fields/fields.f90 index d7878d4..cdade64 100644 --- a/src/fields/fields.f90 +++ b/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, two 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 = gam * vcell - - peref = - b(1) * vcell * zmean - - do k = 2, ngrid - zmean = zmean + gam * two * 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) diff --git a/src/ps3d.f90 b/src/ps3d.f90 index 4b1e94a..fe2a16b 100644 --- a/src/ps3d.f90 +++ b/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,8 +82,6 @@ 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) @@ -91,8 +89,8 @@ subroutine pre_run 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 diff --git a/src/utils/Makefile.am b/src/utils/Makefile.am index 5c200e6..a0e7e0d 100644 --- a/src/utils/Makefile.am +++ b/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 diff --git a/src/utils/ape_density.f90 b/src/utils/ape_density.f90 new file mode 100644 index 0000000..7192f46 --- /dev/null +++ b/src/utils/ape_density.f90 @@ -0,0 +1,26 @@ +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 + +#ifdef ENABLE_IW_TEST_CASE + a = f18 * (b - four * z) ** 2 +#elseif ENABLE_RT_TEST_CASE + a = b * dasin(b) + dsqrt(one - b ** 2) - z * b - 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 From e66be9ad87feb9cdef060d659dfbb3791890c2a0 Mon Sep 17 00:00:00 2001 From: Matthias Frey Date: Wed, 1 Feb 2023 15:57:31 +0000 Subject: [PATCH 2/3] fix preprocessort directive --- src/utils/ape_density.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utils/ape_density.f90 b/src/utils/ape_density.f90 index 7192f46..5eed501 100644 --- a/src/utils/ape_density.f90 +++ b/src/utils/ape_density.f90 @@ -13,7 +13,7 @@ elemental function ape_den(b, z) result(a) #ifdef ENABLE_IW_TEST_CASE a = f18 * (b - four * z) ** 2 -#elseif ENABLE_RT_TEST_CASE +#elif ENABLE_RT_TEST_CASE a = b * dasin(b) + dsqrt(one - b ** 2) - z * b - dcos(z) #else ! dummy line to avoid compiler warning of 'unused variables' From b7058c48692db3bce365a5890a7eb48e0198db42 Mon Sep 17 00:00:00 2001 From: Matthias Frey Date: Fri, 3 Mar 2023 08:47:43 +0000 Subject: [PATCH 3/3] restrict buoyancy in RT --- src/utils/ape_density.f90 | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/utils/ape_density.f90 b/src/utils/ape_density.f90 index 5eed501..35a928a 100644 --- a/src/utils/ape_density.f90 +++ b/src/utils/ape_density.f90 @@ -10,11 +10,14 @@ 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 - a = b * dasin(b) + dsqrt(one - b ** 2) - z * b - dcos(z) + 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