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, &