From 931ab78e9b7c2fb7c8348299db1e2e5c268dc6ac Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Mon, 15 Dec 2025 19:35:42 -0500 Subject: [PATCH] Adds t17d and t20d diagnostics t17d and t20d are the depth of the 17 and 20 degC isotherms respectively, measured from the surface. Interpolation is done via Recon1d reconstruction and then a new function in the class that solves that inverts the reconstruction. This solver was not in the original class because I had originally only thought it would be used in the new grid generator which interpolates density and not for the reconstructed field. Since some other diagnostics for CMIP will make use of remapping or T, I figured we might as well use the reconstructions to do the interpolation. - adds to new registrations and posts in MOM_diagnostics - adds function `%x(self, k, val)` to the Recon1d class. --- src/ALE/Recon1d_PCM.F90 | 31 ++++++++++++ src/ALE/Recon1d_PLM_CW.F90 | 38 ++++++++++++++ src/ALE/Recon1d_PLM_CWK.F90 | 1 + src/ALE/Recon1d_PLM_WLS.F90 | 11 ++++ src/ALE/Recon1d_PLM_hybgen.F90 | 11 ++++ src/ALE/Recon1d_PPM_CW.F90 | 14 +++++- src/ALE/Recon1d_PPM_CWK.F90 | 51 ++++++++++++++++++- src/ALE/Recon1d_PPM_hybgen.F90 | 2 +- src/ALE/Recon1d_type.F90 | 75 ++++++++++++++++++++++++++- src/diagnostics/MOM_diagnostics.F90 | 78 +++++++++++++++++++++++++++-- 10 files changed, 303 insertions(+), 9 deletions(-) diff --git a/src/ALE/Recon1d_PCM.F90 b/src/ALE/Recon1d_PCM.F90 index 3b64844983..d7b386fc23 100644 --- a/src/ALE/Recon1d_PCM.F90 +++ b/src/ALE/Recon1d_PCM.F90 @@ -17,6 +17,7 @@ module Recon1d_PCM !! average() *locally defined !! f() *locally defined !! dfdx() *locally defined +!! - x() *locally defined !! check_reconstruction() *locally defined !! unit_tests() *locally defined !! destroy() *locally defined @@ -36,6 +37,8 @@ module Recon1d_PCM procedure :: f => f !> Implementation of the derivative of the PCM reconstruction at a point [A] procedure :: dfdx => dfdx + !> Implementation of solver for x: f(x)=t + procedure :: x => x !> Implementation of deallocation for PCM procedure :: destroy => destroy !> Implementation of check reconstruction for the PCM reconstruction @@ -105,6 +108,24 @@ real function dfdx(this, k, x) end function dfdx +!> Solver for x: f(x)=t +real function x(this, k, t) + class(PCM), intent(in) :: this !< This reconstruction + integer, intent(in) :: k !< Cell number + real, intent(in) :: t !< Value to solve for [A] + real :: slp ! Difference across cell [A] + + slp = this%u_mean(min(k+1,this%n)) - this%u_mean(max(k-1,1)) + if ( abs(slp) > 0. ) slp = sign(1., slp) + x = 0.5 ! Fall back if t==u_mean + ! if t>u_mean & slp=1 then x=1 + ! if tu_mean & slp=-1 then x=0 + ! if t 0. ) x = 0.5 + slp * sign(0.5, t - this%u_mean(k)) +end function x + !> Average between xa and xb for cell k of a 1D PCM reconstruction [A] real function average(this, k, xa, xb) class(PCM), intent(in) :: this !< This reconstruction @@ -181,6 +202,16 @@ logical function unit_tests(this, verbose, stdout, stderr) call test%real_arr(3, um, (/0.,0.,0./), 'dfdx in center') call test%real_arr(3, ur, (/0.,0.,0./), 'dfdx on right edge') + call test%real_scalar( this%x(1,0.), 0., 'f-1(1,0)=0') + call test%real_scalar( this%x(1,1.), 0.5, 'f-1(1,1)=0.5') + call test%real_scalar( this%x(1,3.), 1., 'f-1(1,3)=1') + call test%real_scalar( this%x(2,1.), 0., 'f-1(2,1)=0') + call test%real_scalar( this%x(2,3.), 0.5, 'f-1(2,3)=0.5') + call test%real_scalar( this%x(2,5.), 1., 'f-1(2,5)=1') + call test%real_scalar( this%x(3,3.), 0., 'f-1(3,3)=0') + call test%real_scalar( this%x(3,5.), 0.5, 'f-1(3,5)=0.5') + call test%real_scalar( this%x(3,7.), 1., 'f-1(3,7)=1') + do k = 1, 3 um(k) = this%average(k, 0.5, 0.75) enddo diff --git a/src/ALE/Recon1d_PLM_CW.F90 b/src/ALE/Recon1d_PLM_CW.F90 index 0c53246286..95e943266e 100644 --- a/src/ALE/Recon1d_PLM_CW.F90 +++ b/src/ALE/Recon1d_PLM_CW.F90 @@ -23,6 +23,7 @@ module Recon1d_PLM_CW !! - average() *locally defined !! - f() *locally defined !! - dfdx() *locally defined +!! - x() *locally defined !! - check_reconstruction() *locally defined !! - unit_tests() *locally defined !! - destroy() *locally defined @@ -45,6 +46,8 @@ module Recon1d_PLM_CW procedure :: f => f !> Implementation of the derivative of the PLM_CW reconstruction at a point [A] procedure :: dfdx => dfdx + !> Implementation of solver for x: f(x)=t + procedure :: x => x !> Implementation of deallocation for PLM_CW procedure :: destroy => destroy !> Implementation of check reconstruction for the PLM_CW reconstruction @@ -194,6 +197,31 @@ real function dfdx(this, k, x) end function dfdx +!> Solver for x such that f(x)=t +real function x(this, k, t) + class(PLM_CW), intent(in) :: this !< This reconstruction + integer, intent(in) :: k !< Cell number + real, intent(in) :: t !< Value to solve for [A] + real :: slp ! Difference across cell [A] + + slp = this%ur(k) - this%ul(k) + if ( abs(slp) > 0. ) then + x = ( t - this%ul(k) ) / slp + x = max( 0., min( x, 1. ) ) + else + slp = this%ul(min(k+1,this%n)) - this%ur(max(k-1,1)) + if ( abs(slp) > 0. ) slp = sign(1., slp) + x = 0.5 ! Fall back if t==u_mean + ! if t>u_mean & slp=1 then x=1 + ! if tu_mean & slp=-1 then x=0 + ! if t 0. ) x = 0.5 + slp * sign(0.5, t - this%u_mean(k)) + endif +end function x + !> Average between xa and xb for cell k of a 1D PLM reconstruction [A] real function average(this, k, xa, xb) class(PLM_CW), intent(in) :: this !< This reconstruction @@ -334,6 +362,16 @@ logical function unit_tests(this, verbose, stdout, stderr) call test%real_arr(3, um, (/0.,2.,0./), 'dfdx in center') call test%real_arr(3, ur, (/0.,2.,0./), 'dfdx on right edge') + call test%real_scalar( this%x(1,0.), 0., 'f-1(1,0)=0') + call test%real_scalar( this%x(1,1.), 0.5, 'f-1(1,1)=0.5') + call test%real_scalar( this%x(1,3.), 1., 'f-1(1,3)=1') + call test%real_scalar( this%x(2,1.), 0., 'f-1(2,1)=0') + call test%real_scalar( this%x(2,3.), 0.5, 'f-1(2,3)=0.5') + call test%real_scalar( this%x(2,5.), 1., 'f-1(2,5)=1') + call test%real_scalar( this%x(3,3.), 0., 'f-1(3,3)=0') + call test%real_scalar( this%x(3,5.), 0.5, 'f-1(3,5)=0.5') + call test%real_scalar( this%x(3,7.), 1., 'f-1(3,7)=1') + do k = 1, 3 um(k) = this%average(k, 0.5, 0.75) ! Average from x=0.25 to 0.75 in each cell enddo diff --git a/src/ALE/Recon1d_PLM_CWK.F90 b/src/ALE/Recon1d_PLM_CWK.F90 index b30af80aa1..76d1d286c9 100644 --- a/src/ALE/Recon1d_PLM_CWK.F90 +++ b/src/ALE/Recon1d_PLM_CWK.F90 @@ -29,6 +29,7 @@ module Recon1d_PLM_CWK !! - average() -> recon1d_plm_cw.average() !! - f() -> recon1d_plm_cw.f() !! - dfdx() -> recon1d_plm_cw.dfdx() +!! - x() -> recon1d_plm_cw.x() !! - check_reconstruction() -> recon1d_plm_cw.check_reconstruction() !! - unit_tests() -> recon1d_plm_cw.unit_tests() !! - destroy() -> recon1d_plm_cw.destroy() diff --git a/src/ALE/Recon1d_PLM_WLS.F90 b/src/ALE/Recon1d_PLM_WLS.F90 index fa38c782aa..24d6988f24 100644 --- a/src/ALE/Recon1d_PLM_WLS.F90 +++ b/src/ALE/Recon1d_PLM_WLS.F90 @@ -17,6 +17,7 @@ module Recon1d_PLM_WLS !! - average() *locally defined !! - f() *locally defined !! - dfdx() *locally defined +!! - x() -> recon1d_type.x() !! - check_reconstruction() *locally defined !! - unit_tests() *locally defined !! - destroy() *locally defined @@ -371,6 +372,16 @@ logical function unit_tests(this, verbose, stdout, stderr) enddo call test%real_arr(3, um, (/1.25,3.25,5.25/), "Return interval average") + call test%real_scalar( this%x(1,0.), 0., 'f-1(1,0)=0') + call test%real_scalar( this%x(1,1.), 0.5, 'f-1(1,1)=0.5') + call test%real_scalar( this%x(1,3.), 1., 'f-1(1,3)=1') + call test%real_scalar( this%x(2,1.), 0., 'f-1(2,1)=0') + call test%real_scalar( this%x(2,3.), 0.5, 'f-1(2,3)=0.5') + call test%real_scalar( this%x(2,5.), 1., 'f-1(2,5)=1') + call test%real_scalar( this%x(3,3.), 0., 'f-1(3,3)=0') + call test%real_scalar( this%x(3,5.), 0.5, 'f-1(3,5)=0.5') + call test%real_scalar( this%x(3,7.), 1., 'f-1(3,7)=1') + call this%destroy() deallocate( um, ul, ur, ull, urr ) diff --git a/src/ALE/Recon1d_PLM_hybgen.F90 b/src/ALE/Recon1d_PLM_hybgen.F90 index 0cf2e8e001..9cae0e4218 100644 --- a/src/ALE/Recon1d_PLM_hybgen.F90 +++ b/src/ALE/Recon1d_PLM_hybgen.F90 @@ -29,6 +29,7 @@ module Recon1d_PLM_hybgen !! - average() *locally defined !! - f() *locally defined !! - dfdx() *locally defined +!! - x() -> recon1d_plm_cw.x() !! - check_reconstruction() *locally defined !! - unit_tests() *locally defined !! - destroy() *locally defined @@ -358,6 +359,16 @@ logical function unit_tests(this, verbose, stdout, stderr) call test%real_arr(3, um, (/0.,2.,0./), 'dfdx in center') call test%real_arr(3, ur, (/0.,2.,0./), 'dfdx on right edge') + call test%real_scalar( this%x(1,0.), 0., 'f-1(1,0)=0') + call test%real_scalar( this%x(1,1.), 0.5, 'f-1(1,1)=0.5') + call test%real_scalar( this%x(1,3.), 1., 'f-1(1,3)=1') + call test%real_scalar( this%x(2,1.), 0., 'f-1(2,1)=0') + call test%real_scalar( this%x(2,3.), 0.5, 'f-1(2,3)=0.5') + call test%real_scalar( this%x(2,5.), 1., 'f-1(2,5)=1') + call test%real_scalar( this%x(3,3.), 0., 'f-1(3,3)=0') + call test%real_scalar( this%x(3,5.), 0.5, 'f-1(3,5)=0.5') + call test%real_scalar( this%x(3,7.), 1., 'f-1(3,7)=1') + do k = 1, 3 um(k) = this%average(k, 0.5, 0.75) ! Average from x=0.25 to 0.75 in each cell enddo diff --git a/src/ALE/Recon1d_PPM_CW.F90 b/src/ALE/Recon1d_PPM_CW.F90 index 7e25bc3d49..1241ffa861 100644 --- a/src/ALE/Recon1d_PPM_CW.F90 +++ b/src/ALE/Recon1d_PPM_CW.F90 @@ -26,6 +26,7 @@ module Recon1d_PPM_CW !! - average() *locally defined !! - f() *locally defined !! - dfdx() *locally defined +!! - x() *locally defined !! - check_reconstruction() *locally defined !! - unit_tests() *locally defined !! - destroy() *locally defined @@ -49,6 +50,8 @@ module Recon1d_PPM_CW procedure :: f => f !> Implementation of the derivative of the PPM_CW reconstruction at a point [A] procedure :: dfdx => dfdx + !> Implementation of solver for x: f(x)=t +! procedure :: x => x !> Implementation of deallocation for PPM_CW procedure :: destroy => destroy !> Implementation of check reconstruction for the PPM_CW reconstruction @@ -152,7 +155,7 @@ subroutine reconstruct(this, h, u) this%ur(n) = u(n) ! PCM this%ul(n) = u(n) ! PCM - do K = 2, n ! K=2 is interface between cells 1 and 2 + do K = 2, n-1 ! K=2 is interface between cells 1 and 2 u0 = u(k-1) u1 = u(k) u2 = u(k+1) @@ -333,6 +336,7 @@ logical function unit_tests(this, verbose, stdout, stderr) call test%set( stdout=stdout ) ! Sets the stdout channel in test call test%set( stderr=stderr ) ! Sets the stderr channel in test call test%set( verbose=verbose ) ! Sets the verbosity flag in test +call test%set( stop_instantly=.true. ) if (verbose) write(stdout,'(a)') 'PPM_CW:unit_tests testing with linear fn' @@ -366,6 +370,10 @@ logical function unit_tests(this, verbose, stdout, stderr) call test%real_arr(5, um, (/0.,3.,3.,3.,0./), 'dfdx in center') call test%real_arr(5, ur, (/0.,3.,3.,3.,0./), 'dfdx on right edge') + call test%real_scalar( this%x(2,1.), 0., 'f-1(2,1)=0') + call test%real_scalar( this%x(2,4.), 0.5, 'f-1(2,4)=0.5') + call test%real_scalar( this%x(2,5.5), 1., 'f-1(2,5.5)=1') + do k = 1, 5 um(k) = this%average(k, 0.5, 0.75) ! Average from x=0.25 to 0.75 in each cell enddo @@ -392,6 +400,10 @@ logical function unit_tests(this, verbose, stdout, stderr) call test%real_arr(5, um, (/1.,6.75,18.75,36.75,61./), 'Return center') call test%real_arr(5, ur, (/1.,12.,27.,48.,61./), 'Return right edge') + call test%real_scalar( this%x(3,12.), 0., 'f-1(3,12)=0') + call test%real_scalar( this%x(3,18.75), 0.5, 'f-1(3,18.75)=0.5', robits=1) + call test%real_scalar( this%x(3,27.), 1., 'f-1(3,27)=1') + ! x = 3 i i=0 at origin ! f(x) = x^2 / 3 = 3 i^2 ! f[i] = [ ( 3 i )^3 - ( 3 i - 3 )^3 ] i=1,2,3,4,5 diff --git a/src/ALE/Recon1d_PPM_CWK.F90 b/src/ALE/Recon1d_PPM_CWK.F90 index 7e0d613e7a..923756cac0 100644 --- a/src/ALE/Recon1d_PPM_CWK.F90 +++ b/src/ALE/Recon1d_PPM_CWK.F90 @@ -27,6 +27,7 @@ module Recon1d_PPM_CWK !! - average() *locally defined !! - f() *locally defined !! - dfdx() *locally defined +!! - x() *locally defined !! - check_reconstruction() *locally defined !! - unit_tests() *locally defined !! - destroy() *locally defined @@ -50,6 +51,8 @@ module Recon1d_PPM_CWK procedure :: f => f !> Implementation of the derivative of the PPM_CWK reconstruction at a point [A] procedure :: dfdx => dfdx + !> Implementation of solver for x: f(x)=t + procedure :: x => x !> Implementation of deallocation for PPM_CWK procedure :: destroy => destroy !> Implementation of check reconstruction for the PPM_CWK reconstruction @@ -137,7 +140,7 @@ subroutine reconstruct(this, h, u) this%ur(n) = u(n) ! PCM this%ul(n) = u(n) ! PCM - do K = 2, n ! K=2 is interface between cells 1 and 2 + do K = 2, n-1 ! K=2 is interface between cells 1 and 2 u0 = u(k-1) u1 = u(k) u2 = u(k+1) @@ -215,6 +218,48 @@ real function dfdx(this, k, x) end function dfdx +!> Solver for x: f(x)=t +real function x(this, k, t) + class(PPM_CWK), intent(in) :: this !< This reconstruction + integer, intent(in) :: k !< Cell number + real, intent(in) :: t !< Value to solve for [A] + real :: xl, xr, xo ! Left/right bounds and guess [nondim] + real :: fl, fr ! Left right values [A] + real :: slp ! Difference across cell or derivative wrt nondim x [A] + real :: f_at_x ! Value at current x [A] + integer :: iter + + slp = this%ur(k) - this%ul(k) + if ( abs(slp) > 0. ) then + xl = 0. + xr = 1. + xo = ( t - this%ul(k) ) / slp ! First guess by regula falsi + f_at_x = this%f(k, xo) + do iter = 1,5 + slp = this%dfdx(k, xo) + x = xo - ( f_at_x - t ) / slp ! Newton-Raphson step + if ( x < xl ) x = 0.5 * ( xl + xo ) ! Replace with bi-section + if ( x > xr ) x = 0.5 * ( xr + xo ) ! Replace with bi-section + f_at_x = this%f(k, x) + if ( abs(f_at_x - t) <= 0. .or. abs(x - xo) < this%x_tolerance ) return + if ( f_at_x < t ) xl = x ! Replace left bound + if ( f_at_x > t ) xr = x ! Replace right bound + xo = x + enddo + else + x = 0.5 ! Fall back if t==u_mean + slp = this%ul(min(k+1,this%n)) - this%ur(max(k-1,1)) + if ( abs(slp) > 0. ) slp = sign(1., slp) + ! if t>u_mean & slp=1 then x=1 + ! if tu_mean & slp=-1 then x=0 + ! if t Average between xa and xb for cell k of a 1D PPM reconstruction [A] real function average(this, k, xa, xb) class(PPM_CWK), intent(in) :: this !< This reconstruction @@ -373,6 +418,10 @@ logical function unit_tests(this, verbose, stdout, stderr) call test%real_arr(5, um, (/1.,6.75,18.75,36.75,61./), 'Return center') call test%real_arr(5, ur, (/1.,12.,27.,48.,61./), 'Return right edge') + call test%real_scalar( this%x(3,12.), 0., 'f-1(3,12)=0') + call test%real_scalar( this%x(3,18.75), 0.5, 'f-1(3,18.75)=0.5', robits=1) + call test%real_scalar( this%x(3,27.), 1., 'f-1(3,27)=1') + ! x = 3 i i=0 at origin ! f(x) = x^2 / 3 = 3 i^2 ! f[i] = [ ( 3 i )^3 - ( 3 i - 3 )^3 ] i=1,2,3,4,5 diff --git a/src/ALE/Recon1d_PPM_hybgen.F90 b/src/ALE/Recon1d_PPM_hybgen.F90 index 058c0a80dc..e3ee4b8def 100644 --- a/src/ALE/Recon1d_PPM_hybgen.F90 +++ b/src/ALE/Recon1d_PPM_hybgen.F90 @@ -127,7 +127,7 @@ subroutine reconstruct(this, h, u) this%ur(n) = u(n) ! PCM this%ul(n) = u(n) ! PCM - do K = 2, n ! K=2 is interface between cells 1 and 2 + do K = 2, n-1 ! K=2 is interface between cells 1 and 2 u0 = u(k-1) u1 = u(k) u2 = u(k+1) diff --git a/src/ALE/Recon1d_type.F90 b/src/ALE/Recon1d_type.F90 index c11d880cc8..f475480a53 100644 --- a/src/ALE/Recon1d_type.F90 +++ b/src/ALE/Recon1d_type.F90 @@ -16,6 +16,7 @@ module Recon1d_type integer :: n = 0 !< Number of cells in column real, allocatable, dimension(:) :: u_mean !< Cell mean [A] real :: h_neglect = 0. !< A negligibly small width used in cell reconstructions in the same units as h [H] + real :: x_tolerance = 1. * epsilon(1.) !< Solver tolerance for x in element (0,1) [nondim] logical :: check = .false. !< If true, enable some consistency checking logical :: debug = .false. !< If true, dump info as calculations are made (do not enable) @@ -50,6 +51,8 @@ module Recon1d_type ! The following functions/subroutines are shared across all reconstructions and provided by this module ! unless replaced for the purpose of optimization + !> Solves for x such that f(x)=t + procedure :: x => x !> Remaps the column to subgrid h_sub procedure :: remap_to_sub_grid => remap_to_sub_grid !> Set debugging @@ -97,7 +100,7 @@ end function i_average !> Point-wise value of reconstruction [A] !! - !! THe function is only valid for 0 <= x <= 1. x is effectively clipped to this range. + !! The function is only valid for 0 <= x <= 1. x is effectively clipped to this range. real function i_f(this, k, x) import :: Recon1d class(Recon1d), intent(in) :: this !< This reconstruction @@ -107,7 +110,7 @@ end function i_f !> Point-wise value of derivative reconstruction [A] !! - !! THe function is only valid for 0 <= x <= 1. x is effectively clipped to this range. + !! The function is only valid for 0 <= x <= 1. x is effectively clipped to this range. real function i_dfdx(this, k, x) import :: Recon1d class(Recon1d), intent(in) :: this !< This reconstruction @@ -115,6 +118,18 @@ real function i_dfdx(this, k, x) real, intent(in) :: x !< Non-dimensional position within element [nondim] end function i_dfdx + !> Point-wise solver for x: f(x)=t [nondim] + !! + !! The function solves for the non-dimensional position x within the cell where + !! the reconstruction f(x)=t. The solver returns x=0 or x=1 if the target, t, + !! is outside of the cell. + real function i_x(this, k, t) + import :: Recon1d + class(Recon1d), intent(in) :: this !< This reconstruction + integer, intent(in) :: k !< Cell number + real, intent(in) :: t !< Value to solve for [A] + end function i_x + !> Returns true if some inconsistency is detected, false otherwise !! !! The nature of "consistency" is defined by the implementations @@ -164,6 +179,62 @@ end function i_unit_tests contains +!> Solve for x such that f(x)=t +!! +!! This solver uses bounded Newton-Raphson method with a fixed +!! number of iterations +real function x(this, k, t) + class(Recon1d), intent(in) :: this !< This reconstruction + integer, intent(in) :: k !< Cell number + real, intent(in) :: t !< Value to solve for [A] + real :: xl, xr, xo ! Left/right bounds and guess [nondim] + real :: fl, fr ! Left right values [A] + real :: slp ! Difference across cell or derivative wrt nondim x [A] + real :: f_at_x ! Value at current x [A] + integer :: iter + + x = 0.5 ! Fall back for special conditions + fl = this%f(k, 0.) + fr = this%f(k, 1.) + slp = fr - fl + if ( ( fl - t ) * ( t - fr ) > 0. ) then + ! t is inside the range fl..fr + xl = 0. + xr = 1. + xo = ( t - this%f(k, 0.) ) / slp ! First guess by regula falsi + f_at_x = this%f(k, xo) + do iter = 1,10 + slp = this%dfdx(k, xo) + x = xo - ( f_at_x - t ) / slp ! Newton-Raphson step + if ( x < xl ) x = 0.5 * ( xl + xo ) ! Replace with bi-section + if ( x > xr ) x = 0.5 * ( xr + xo ) ! Replace with bi-section + f_at_x = this%f(k, x) + if ( abs(f_at_x - t) <= 0. .or. abs(x - xo) < this%x_tolerance ) return + if ( f_at_x < t ) xl = x ! Replace left bound + if ( f_at_x > t ) xr = x ! Replace right bound + xo = x + enddo + elseif ( abs(slp) > 0. ) then + slp = sign(1., slp) + ! if t>u_mean & slp=1 then x=1 + ! if tu_mean & slp=-1 then x=0 + ! if t 0. ) slp = sign(1., slp) + ! if t>u_mean & slp=1 then x=1 + ! if tu_mean & slp=-1 then x=0 + ! if t 0. ) x = 0.5 + slp * sign(0.5, t - this%u_mean(k)) + endif +end function x + !> Remaps the column to subgrid h_sub !! !! It is assumed that h_sub is a perfect sub-grid of h0, meaning each h0 cell diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 6c220c79cf..36437b8603 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -20,6 +20,7 @@ module MOM_diagnostics use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type use MOM_domains, only : To_North, To_East use MOM_EOS, only : calculate_density, calculate_density_derivs, EOS_domain +use MOM_EOS, only : calculate_spec_vol use MOM_EOS, only : cons_temp_to_pot_temp, abs_saln_to_prac_saln use MOM_error_handler, only : MOM_error, FATAL, WARNING use MOM_file_parser, only : get_param, log_version, param_file_type @@ -33,6 +34,7 @@ module MOM_diagnostics use MOM_variables, only : accel_diag_ptrs, cont_diag_ptrs, surface use MOM_verticalGrid, only : verticalGrid_type, get_thickness_units, get_flux_units use MOM_wave_speed, only : wave_speed, wave_speed_CS, wave_speed_init +use Recon1d_EPPM_CWK, only : EPPM_CWK implicit none ; private @@ -114,6 +116,7 @@ module MOM_diagnostics integer :: id_drho_dT = -1, id_drho_dS = -1 integer :: id_h_pre_sync = -1 integer :: id_tosq = -1, id_sosq = -1 + integer :: id_t20d = -1, id_t17d = -1 !>@} type(wave_speed_CS) :: wave_speed !< Wave speed control struct @@ -890,7 +893,7 @@ subroutine calculate_vertical_integrals(h, tv, p_surf, G, GV, US, CS) !! as setting the surface pressure to 0. type(diagnostics_CS), intent(inout) :: CS !< Control structure returned by a !! previous call to diagnostics_init. - + ! Local variables real, dimension(SZI_(G),SZJ_(G)) :: & z_top, & ! Height of the top of a layer or the ocean [Z ~> m]. z_bot, & ! Height of the bottom of a layer (for id_mass) or the @@ -902,10 +905,17 @@ subroutine calculate_vertical_integrals(h, tv, p_surf, G, GV, US, CS) btm_pres,&! The pressure at the ocean bottom, or CMIP variable 'pbo'. ! This is the column mass multiplied by gravity plus the pressure ! at the ocean surface [R L2 T-2 ~> Pa]. - tr_int ! vertical integral of a tracer times density, + tr_int,& ! vertical integral of a tracer times density, ! (Rho_0 in a Boussinesq model) [Conc R Z ~> Conc kg m-2]. - real :: IG_Earth ! Inverse of gravitational acceleration [T2 Z L-2 ~> s2 m-1]. - + d17,& ! Depth of 17 degC isotherm [Z ~> m] + d20 ! Depth of 20 degC isotherm [Z ~> m] + real :: IG_Earth ! Inverse of gravitational acceleration [T2 Z L-2 ~> s2 m-1]. + real :: Ttop, Tbot ! Temperature at top/bottom of cell [C ~> degC] + type(EPPM_CWK) :: PPM ! Class for reconstruction + real :: d_from_ssh(0:GV%ke) ! eta-z (Distance from surface) [Z ~> m] + real :: dz ! Layer thickness in Z [Z ~> m] + real :: p(0:GV%ke) ! Hydrostatic pressure [R L2 T-2 ~> Pa] + real :: spv(GV%ke) ! Specific volume [R-1 ~> m3 kg-1] integer :: i, j, k, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -950,6 +960,57 @@ subroutine calculate_vertical_integrals(h, tv, p_surf, G, GV, US, CS) endif if (CS%id_col_mass > 0) call post_data(CS%id_col_mass, mass, CS%diag) endif + if (CS%id_t20d > 0 .or. CS%id_t17d > 0) then + call PPM%init(GV%ke, h_neglect=0.) + do j=js,je ; do i=is,ie + ! Pre-calculate the interface depths relative to the surface + if (GV%Boussinesq) then + d_from_ssh(0) = 0. + do k=1,nz + d_from_ssh(k) = d_from_ssh(k-1) + h(i,j,k) * GV%H_to_Z + enddo + else + ! Non-Boussinesq + p(0) = 0. + do k=1,nz + p(k) = p(k-1) + h(i,j,k) * ( GV%H_to_RZ * GV%g_Earth ) ! Pressure at lower interface + p(k-1) = 0.5 * ( p(k-1) + p(k) ) ! EOS fn needs pressure in the layer center + enddo + call calculate_spec_vol(tv%T(i,j,:), tv%S(i,j,:), p(0:nz-1), spv, tv%eqn_of_state) + d_from_ssh(0) = 0. + do k=1,nz + d_from_ssh(k) = d_from_ssh(k-1) + ( h(i,j,k) * GV%H_to_RZ ) * spv(k) + enddo + endif + call PPM%reconstruct(h(i,j,:), tv%T(i,j,:)) + d17(i,j) = d_from_ssh(nz) + d20(i,j) = d_from_ssh(nz) + do k=nz,1,-1 + Ttop = PPM%f(k, 0.) + Tbot = PPM%f(k, 1.) + if ( Tbot>Ttop ) cycle ! The cell is inverted, skip to next + if ( 20.=0 + if ( Tbot<=17. .and. 17.<=Ttop ) then + ! The 17 degC isotherm is within the cell which is non-negatively stratified + d17(i,j) = d_from_ssh(k-1) + dz * PPM%x(k, 17.) + elseif ( Ttop<17. ) then + ! The 17 degC isotherm is above the top of the cell + d17(i,j) = d_from_ssh(k-1) + endif + if ( Tbot<=20. .and. 20.<=Ttop ) then + ! The 20 degC isotherm is within the cell which is non-negatively stratified + d20(i,j) = d_from_ssh(k-1) + dz * PPM%x(k, 20.) + elseif ( Ttop<20. ) then + ! The 20 degC isotherm is above the top of the cell + d20(i,j) = d_from_ssh(k-1) + endif + enddo + enddo ; enddo + call PPM%destroy() + if (CS%id_t17d > 0) call post_data(CS%id_t17d, d17, CS%diag) + if (CS%id_t20d > 0) call post_data(CS%id_t20d, d20, CS%diag) + endif end subroutine calculate_vertical_integrals @@ -1891,6 +1952,15 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag CS%id_abssosga = register_scalar_field('ocean_model', 'ssabss_global', Time, diag, & long_name='Global Area Average Sea Surface Absolute Salinity', & units='psu', conversion=US%S_to_ppt, standard_name='sea_surface_absolute_salinity') + + CS%id_t20d = register_diag_field('ocean_model', 't20d', diag%axesT1, Time, & + 'Depth of 20 degree Celsius Isotherm', & + units='m', conversion=US%Z_to_m, & + standard_name='depth_of_isosurface_of_sea_water_potential_temperature') + CS%id_t17d = register_diag_field('ocean_model', 't17d', diag%axesT1, Time, & + 'Depth of 17 degree Celsius Isotherm', & + units='m', conversion=US%Z_to_m, & + standard_name='depth_of_isosurface_of_sea_water_potential_temperature') endif CS%id_u = register_diag_field('ocean_model', 'u', diag%axesCuL, Time, &