diff --git a/src/ALE/MOM_remapping.F90 b/src/ALE/MOM_remapping.F90 index bf68782438..6e404ff2dd 100644 --- a/src/ALE/MOM_remapping.F90 +++ b/src/ALE/MOM_remapping.F90 @@ -37,8 +37,8 @@ module MOM_remapping ! The following routines are visible to the outside world public remapping_core_h, remapping_core_w -public initialize_remapping, end_remapping, remapping_set_param -public remapping_unit_tests +public initialize_remapping, end_remapping, remapping_set_param, extract_member_remapping_CS +public remapping_unit_tests, build_reconstructions_1d, average_value_ppoly public dzFromH1H2 ! The following are private parameter constants @@ -109,6 +109,26 @@ subroutine remapping_set_param(CS, remapping_scheme, boundary_extrapolation, & endif end subroutine remapping_set_param +subroutine extract_member_remapping_CS(CS, remapping_scheme, degree, boundary_extrapolation, check_reconstruction, & + check_remapping, force_bounds_in_subcell) + type(remapping_CS), intent(in) :: CS !< Control structure for remapping module + integer, optional, intent(out) :: remapping_scheme !< Determines which reconstruction scheme to use + integer, optional, intent(out) :: degree !< Degree of polynomial reconstruction + logical, optional, intent(out) :: boundary_extrapolation !< If true, extrapolate boundaries + logical, optional, intent(out) :: check_reconstruction !< If true, reconstructions are checked for consistency. + logical, optional, intent(out) :: check_remapping !< If true, the result of remapping are checked + !! for conservation and bounds. + logical, optional, intent(out) :: force_bounds_in_subcell !< If true, the intermediate values used in + !! remapping are forced to be bounded. + + if (present(remapping_scheme)) remapping_scheme = CS%remapping_scheme + if (present(degree)) degree = CS%degree + if (present(boundary_extrapolation)) boundary_extrapolation = CS%boundary_extrapolation + if (present(check_reconstruction)) check_reconstruction = CS%check_reconstruction + if (present(check_remapping)) check_remapping = CS%check_remapping + if (present(force_bounds_in_subcell)) force_bounds_in_subcell = CS%force_bounds_in_subcell + +end subroutine extract_member_remapping_CS !> Calculate edge coordinate x from cell width h subroutine buildGridFromH(nz, h, x) integer, intent(in) :: nz !< Number of cells @@ -175,8 +195,7 @@ subroutine remapping_core_h(CS, n0, h0, u0, n1, h1, u1) integer :: k real :: eps, h0tot, h0err, h1tot, h1err, u0tot, u0err, u0min, u0max, u1tot, u1err, u1min, u1max, uh_err - call build_reconstructions_1d( n0, h0, u0, CS%remapping_scheme, CS%degree, CS%boundary_extrapolation, & - ppoly_r_coefficients, ppoly_r_E, ppoly_r_S, iMethod ) + call build_reconstructions_1d( CS, n0, h0, u0, ppoly_r_coefficients, ppoly_r_E, ppoly_r_S, iMethod ) if (CS%check_reconstruction) call check_reconstructions_1d(n0, h0, u0, CS%degree, & CS%boundary_extrapolation, ppoly_r_coefficients, ppoly_r_E, ppoly_r_S) @@ -243,8 +262,7 @@ subroutine remapping_core_w( CS, n0, h0, u0, n1, dx, u1 ) real :: u0tot, u0err, u0min, u0max, u1tot, u1err, u1min, u1max, uh_err real, dimension(n1) :: h1 !< Cell widths on target grid - call build_reconstructions_1d( n0, h0, u0, CS%remapping_scheme, CS%degree, CS%boundary_extrapolation, & - ppoly_r_coefficients, ppoly_r_E, ppoly_r_S, iMethod ) + call build_reconstructions_1d( CS, n0, h0, u0, ppoly_r_coefficients, ppoly_r_E, ppoly_r_S, iMethod ) if (CS%check_reconstruction) call check_reconstructions_1d(n0, h0, u0, CS%degree, & CS%boundary_extrapolation, ppoly_r_coefficients, ppoly_r_E, ppoly_r_S) @@ -301,20 +319,19 @@ subroutine remapping_core_w( CS, n0, h0, u0, n1, dx, u1 ) end subroutine remapping_core_w !> Creates polynomial reconstructions of u0 on the source grid h0. -subroutine build_reconstructions_1d( n0, h0, u0, remapping_scheme, deg, boundary_extrapolation, & - ppoly_r_coefficients, ppoly_r_E, ppoly_r_S, iMethod ) - integer, intent(in) :: n0 !< Number of cells on source grid - real, dimension(n0), intent(in) :: h0 !< Cell widths on source grid - real, dimension(n0), intent(in) :: u0 !< Cell averages on source grid - integer, intent(in) :: remapping_scheme !< Remapping scheme - integer, intent(in) :: deg !< Degree of polynomial reconstruction - logical, intent(in) :: boundary_extrapolation !< Extrapolate at boundaries if true - real, dimension(n0,deg+1),intent(out) :: ppoly_r_coefficients !< Coefficients of polynomial - real, dimension(n0,2), intent(out) :: ppoly_r_E !< Edge value of polynomial - real, dimension(n0,2), intent(out) :: ppoly_r_S !< Edge slope of polynomial - integer, intent(out) :: iMethod !< Integration method +subroutine build_reconstructions_1d( CS, n0, h0, u0, ppoly_r_coefficients, ppoly_r_E, ppoly_r_S, iMethod ) + type(remapping_CS), intent(in) :: CS + integer, intent(in) :: n0 !< Number of cells on source grid + real, dimension(n0), intent(in) :: h0 !< Cell widths on source grid + real, dimension(n0), intent(in) :: u0 !< Cell averages on source grid + real, dimension(n0,CS%degree+1), intent(out) :: ppoly_r_coefficients !< Coefficients of polynomial + real, dimension(n0,2), intent(out) :: ppoly_r_E !< Edge value of polynomial + real, dimension(n0,2), intent(out) :: ppoly_r_S !< Edge slope of polynomial + integer, intent(out) :: iMethod !< Integration method ! Local variables integer :: local_remapping_scheme + integer :: remapping_scheme !< Remapping scheme + logical :: boundary_extrapolation !< Extrapolate at boundaries if true ! Reset polynomial ppoly_r_E(:,:) = 0.0 @@ -322,7 +339,7 @@ subroutine build_reconstructions_1d( n0, h0, u0, remapping_scheme, deg, boundary ppoly_r_coefficients(:,:) = 0.0 iMethod = -999 - local_remapping_scheme = remapping_scheme + local_remapping_scheme = CS%remapping_scheme if (n0<=1) then local_remapping_scheme = REMAPPING_PCM elseif (n0<=3) then @@ -336,21 +353,21 @@ subroutine build_reconstructions_1d( n0, h0, u0, remapping_scheme, deg, boundary iMethod = INTEGRATION_PCM case ( REMAPPING_PLM ) call PLM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_coefficients ) - if ( boundary_extrapolation) then + if ( CS%boundary_extrapolation ) then call PLM_boundary_extrapolation( n0, h0, u0, ppoly_r_E, ppoly_r_coefficients) end if iMethod = INTEGRATION_PLM case ( REMAPPING_PPM_H4 ) call edge_values_explicit_h4( n0, h0, u0, ppoly_r_E ) call PPM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_coefficients ) - if ( boundary_extrapolation) then + if ( CS%boundary_extrapolation ) then call PPM_boundary_extrapolation( n0, h0, u0, ppoly_r_E, ppoly_r_coefficients ) end if iMethod = INTEGRATION_PPM case ( REMAPPING_PPM_IH4 ) call edge_values_implicit_h4( n0, h0, u0, ppoly_r_E ) call PPM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_coefficients ) - if ( boundary_extrapolation) then + if ( CS%boundary_extrapolation ) then call PPM_boundary_extrapolation( n0, h0, u0, ppoly_r_E, ppoly_r_coefficients ) end if iMethod = INTEGRATION_PPM @@ -358,7 +375,7 @@ subroutine build_reconstructions_1d( n0, h0, u0, remapping_scheme, deg, boundary call edge_values_implicit_h4( n0, h0, u0, ppoly_r_E ) call edge_slopes_implicit_h3( n0, h0, u0, ppoly_r_S ) call PQM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_S, ppoly_r_coefficients ) - if ( boundary_extrapolation) then + if ( CS%boundary_extrapolation ) then call PQM_boundary_extrapolation_v1( n0, h0, u0, ppoly_r_E, ppoly_r_S, ppoly_r_coefficients ) end if iMethod = INTEGRATION_PQM @@ -366,7 +383,7 @@ subroutine build_reconstructions_1d( n0, h0, u0, remapping_scheme, deg, boundary call edge_values_implicit_h6( n0, h0, u0, ppoly_r_E ) call edge_slopes_implicit_h5( n0, h0, u0, ppoly_r_S ) call PQM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_S, ppoly_r_coefficients ) - if ( boundary_extrapolation) then + if ( CS%boundary_extrapolation ) then call PQM_boundary_extrapolation_v1( n0, h0, u0, ppoly_r_E, ppoly_r_S, ppoly_r_coefficients ) end if iMethod = INTEGRATION_PQM diff --git a/src/ALE/polynomial_functions.F90 b/src/ALE/polynomial_functions.F90 index ee90ec7e60..5fdc63e1ca 100644 --- a/src/ALE/polynomial_functions.F90 +++ b/src/ALE/polynomial_functions.F90 @@ -12,7 +12,7 @@ module polynomial_functions implicit none ; private -public :: evaluation_polynomial, integration_polynomial +public :: evaluation_polynomial, integration_polynomial, first_derivative_polynomial ! ----------------------------------------------------------------------------- ! This module contains the following routines @@ -51,6 +51,37 @@ real function evaluation_polynomial( coefficients, nb_coefficients, x ) end function evaluation_polynomial +!> Calculates the first derivative of a polynomial with coefficients as above +!! evaluated at a point x +real function first_derivative_polynomial( coefficients, nb_coefficients, x ) +! ----------------------------------------------------------------------------- +! The polynomial is defined by the coefficients contained in the +! array of the same name, as follows: C(1) + C(2)x + C(3)x^2 + C(4)x^3 + ... +! where C refers to the array 'coefficients'. +! The number of coefficients is given by nb_coefficients and x +! is the coordinate where the polynomial's derivative is to be evaluated. +! +! The function returns the value of the polynomial at x. +! ----------------------------------------------------------------------------- + + ! Arguments + real, dimension(:), intent(in) :: coefficients + integer, intent(in) :: nb_coefficients + real, intent(in) :: x + + ! Local variables + integer :: k + real :: f ! value of polynomial at x + + f = 0.0 + do k = 2,nb_coefficients + f = f + REAL(k-1)*coefficients(k) * ( x**(k-2) ) + end do + + first_derivative_polynomial = f + +end function first_derivative_polynomial + ! ----------------------------------------------------------------------------- ! Exact integration of polynomial of degree n ! ----------------------------------------------------------------------------- diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90 index 149fe01bb1..d3d6d7dc33 100644 --- a/src/equation_of_state/MOM_EOS.F90 +++ b/src/equation_of_state/MOM_EOS.F90 @@ -3,20 +3,26 @@ module MOM_EOS ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_EOS_linear, only : calculate_density_scalar_linear, calculate_density_array_linear -use MOM_EOS_linear, only : calculate_density_derivs_linear, calculate_specvol_derivs_linear, int_density_dz_linear +use MOM_EOS_linear, only : calculate_density_linear +use MOM_EOS_linear, only : calculate_density_derivs_linear +use MOM_EOS_linear, only : calculate_specvol_derivs_linear, int_density_dz_linear +use MOM_EOS_linear, only : calculate_density_second_derivs_linear use MOM_EOS_linear, only : calculate_compress_linear, int_spec_vol_dp_linear -use MOM_EOS_Wright, only : calculate_density_scalar_wright, calculate_density_array_wright -use MOM_EOS_Wright, only : calculate_density_derivs_wright, calculate_specvol_derivs_wright, int_density_dz_wright +use MOM_EOS_Wright, only : calculate_density_wright, calculate_density_wright +use MOM_EOS_Wright, only : calculate_density_derivs_wright +use MOM_EOS_Wright, only : calculate_specvol_derivs_wright, int_density_dz_wright use MOM_EOS_Wright, only : calculate_compress_wright, int_spec_vol_dp_wright -use MOM_EOS_UNESCO, only : calculate_density_scalar_unesco, calculate_density_array_unesco +use MOM_EOS_Wright, only : calculate_density_second_derivs_wright +use MOM_EOS_UNESCO, only : calculate_density_unesco use MOM_EOS_UNESCO, only : calculate_density_derivs_unesco, calculate_density_unesco use MOM_EOS_UNESCO, only : calculate_compress_unesco -use MOM_EOS_NEMO, only : calculate_density_scalar_nemo, calculate_density_array_nemo +use MOM_EOS_NEMO, only : calculate_density_nemo use MOM_EOS_NEMO, only : calculate_density_derivs_nemo, calculate_density_nemo use MOM_EOS_NEMO, only : calculate_compress_nemo -use MOM_EOS_TEOS10, only : calculate_density_scalar_teos10, calculate_density_array_teos10 -use MOM_EOS_TEOS10, only : calculate_density_derivs_teos10, calculate_specvol_derivs_teos10 +use MOM_EOS_TEOS10, only : calculate_density_teos10 +use MOM_EOS_TEOS10, only : calculate_density_derivs_teos10 +use MOM_EOS_TEOS10, only : calculate_specvol_derivs_teos10 +use MOM_EOS_TEOS10, only : calculate_density_second_derivs_teos10 use MOM_EOS_TEOS10, only : calculate_compress_teos10 use MOM_EOS_TEOS10, only : gsw_sp_from_sr, gsw_pt_from_ct use MOM_TFreeze, only : calculate_TFreeze_linear, calculate_TFreeze_Millero, calculate_TFreeze_teos10 @@ -30,8 +36,8 @@ module MOM_EOS #include public calculate_compress, calculate_density, query_compressible -public calculate_density_derivs, calculate_specific_vol_derivs -public EOS_init, EOS_end, EOS_allocate +public calculate_density_derivs, calculate_specific_vol_derivs, calculate_density_second_derivs +public EOS_init, EOS_manual_init, EOS_end, EOS_allocate public EOS_use_linear public int_density_dz, int_specific_vol_dp public int_density_dz_generic_plm, int_density_dz_generic_ppm @@ -40,12 +46,21 @@ module MOM_EOS public calculate_TFreeze public convert_temp_salt_for_TEOS10 public gsw_sp_from_sr, gsw_pt_from_ct +public extract_member_EOS !> Calculates density of sea water from T, S and P interface calculate_density module procedure calculate_density_scalar, calculate_density_array end interface calculate_density +interface calculate_density_derivs + module procedure calculate_density_derivs_scalar, calculate_density_derivs_array +end interface calculate_density_derivs + +interface calculate_density_second_derivs + module procedure calculate_density_second_derivs_scalar, calculate_density_second_derivs_array +end interface calculate_density_second_derivs + !> Calculates the freezing point of sea water from T, S and P interface calculate_TFreeze module procedure calculate_TFreeze_scalar, calculate_TFreeze_array @@ -72,11 +87,11 @@ module MOM_EOS end type EOS_type ! The named integers that might be stored in eqn_of_state_type%form_of_EOS. -integer, parameter :: EOS_LINEAR = 1 -integer, parameter :: EOS_UNESCO = 2 -integer, parameter :: EOS_WRIGHT = 3 -integer, parameter :: EOS_TEOS10 = 4 -integer, parameter :: EOS_NEMO = 5 +integer, parameter, public :: EOS_LINEAR = 1 +integer, parameter, public :: EOS_UNESCO = 2 +integer, parameter, public :: EOS_WRIGHT = 3 +integer, parameter, public :: EOS_TEOS10 = 4 +integer, parameter, public :: EOS_NEMO = 5 character*(10), parameter :: EOS_LINEAR_STRING = "LINEAR" character*(10), parameter :: EOS_UNESCO_STRING = "UNESCO" @@ -108,16 +123,16 @@ subroutine calculate_density_scalar(T, S, pressure, rho, EOS) select case (EOS%form_of_EOS) case (EOS_LINEAR) - call calculate_density_scalar_linear(T, S, pressure, rho, & + call calculate_density_linear(T, S, pressure, rho, & EOS%Rho_T0_S0, EOS%dRho_dT, EOS%dRho_dS) case (EOS_UNESCO) - call calculate_density_scalar_unesco(T, S, pressure, rho) + call calculate_density_unesco(T, S, pressure, rho) case (EOS_WRIGHT) - call calculate_density_scalar_wright(T, S, pressure, rho) + call calculate_density_wright(T, S, pressure, rho) case (EOS_TEOS10) - call calculate_density_scalar_teos10(T, S, pressure, rho) + call calculate_density_teos10(T, S, pressure, rho) case (EOS_NEMO) - call calculate_density_scalar_nemo(T, S, pressure, rho) + call calculate_density_nemo(T, S, pressure, rho) case default call MOM_error(FATAL, & "calculate_density_scalar: EOS is not valid.") @@ -140,16 +155,16 @@ subroutine calculate_density_array(T, S, pressure, rho, start, npts, EOS) select case (EOS%form_of_EOS) case (EOS_LINEAR) - call calculate_density_array_linear(T, S, pressure, rho, start, npts, & + call calculate_density_linear(T, S, pressure, rho, start, npts, & EOS%Rho_T0_S0, EOS%dRho_dT, EOS%dRho_dS) case (EOS_UNESCO) - call calculate_density_array_unesco(T, S, pressure, rho, start, npts) + call calculate_density_unesco(T, S, pressure, rho, start, npts) case (EOS_WRIGHT) - call calculate_density_array_wright(T, S, pressure, rho, start, npts) + call calculate_density_wright(T, S, pressure, rho, start, npts) case (EOS_TEOS10) - call calculate_density_array_teos10(T, S, pressure, rho, start, npts) + call calculate_density_teos10(T, S, pressure, rho, start, npts) case (EOS_NEMO) - call calculate_density_array_nemo (T, S, pressure, rho, start, npts) + call calculate_density_nemo (T, S, pressure, rho, start, npts) case default call MOM_error(FATAL, & "calculate_density_array: EOS%form_of_EOS is not valid.") @@ -210,7 +225,7 @@ subroutine calculate_TFreeze_array(S, pressure, T_fr, start, npts, EOS) end subroutine calculate_TFreeze_array !> Calls the appropriate subroutine to calculate density derivatives for 1-D array inputs. -subroutine calculate_density_derivs(T, S, pressure, drho_dT, drho_dS, start, npts, EOS) +subroutine calculate_density_derivs_array(T, S, pressure, drho_dT, drho_dS, start, npts, EOS) real, dimension(:), intent(in) :: T !< Potential temperature referenced to the surface (degC) real, dimension(:), intent(in) :: S !< Salinity (PSU) real, dimension(:), intent(in) :: pressure !< Pressure (Pa) @@ -225,8 +240,8 @@ subroutine calculate_density_derivs(T, S, pressure, drho_dT, drho_dS, start, npt select case (EOS%form_of_EOS) case (EOS_LINEAR) - call calculate_density_derivs_linear(T, S, pressure, drho_dT, drho_dS, start, & - npts, EOS%Rho_T0_S0, EOS%dRho_dT, EOS%dRho_dS) + call calculate_density_derivs_linear(T, S, pressure, drho_dT, drho_dS, EOS%Rho_T0_S0, EOS%dRho_dT, EOS%dRho_dS, & + start, npts) case (EOS_UNESCO) call calculate_density_derivs_unesco(T, S, pressure, drho_dT, drho_dS, start, npts) case (EOS_WRIGHT) @@ -235,12 +250,106 @@ subroutine calculate_density_derivs(T, S, pressure, drho_dT, drho_dS, start, npt call calculate_density_derivs_teos10(T, S, pressure, drho_dT, drho_dS, start, npts) case (EOS_NEMO) call calculate_density_derivs_nemo(T, S, pressure, drho_dT, drho_dS, start, npts) + case default + call MOM_error(FATAL, & + "calculate_density_derivs_array: EOS%form_of_EOS is not valid.") + end select + +end subroutine calculate_density_derivs_array + +!> Calls the appropriate subroutines to calculate density derivatives by promoting a scalar to a one-element array +subroutine calculate_density_derivs_scalar(T, S, pressure, drho_dT, drho_dS, EOS) + real, intent(in) :: T !< Potential temperature referenced to the surface (degC) + real, intent(in) :: S !< Salinity (PSU) + real, intent(in) :: pressure !< Pressure (Pa) + real, intent(out) :: drho_dT !< The partial derivative of density with potential tempetature, in kg m-3 K-1. + real, intent(out) :: drho_dS !< The partial derivative of density with salinity, in kg m-3 psu-1. + type(EOS_type), pointer :: EOS !< Equation of state structure + if (.not.associated(EOS)) call MOM_error(FATAL, & + "calculate_density_derivs called with an unassociated EOS_type EOS.") + + select case (EOS%form_of_EOS) + case (EOS_LINEAR) + call calculate_density_derivs_linear(T, S, pressure, drho_dT, drho_dS, & + EOS%Rho_T0_S0, EOS%dRho_dT, EOS%dRho_dS) + case (EOS_WRIGHT) + call calculate_density_derivs_wright(T, S, pressure, drho_dT, drho_dS) + case (EOS_TEOS10) + call calculate_density_derivs_teos10(T, S, pressure, drho_dT, drho_dS) + case default + call MOM_error(FATAL, & + "calculate_density_derivs_scalar: EOS%form_of_EOS is not valid.") + end select + +end subroutine calculate_density_derivs_scalar + +!> Calls the appropriate subroutine to calculate density second derivatives for 1-D array inputs. +subroutine calculate_density_second_derivs_array(T, S, pressure, drho_dS_dS, drho_dS_dT, drho_dT_dT, drho_dS_dP, drho_dT_dP, & + start, npts, EOS) + real, dimension(:), intent(in) :: T !< Potential temperature referenced to the surface (degC) + real, dimension(:), intent(in) :: S !< Salinity (PSU) + real, dimension(:), intent(in) :: pressure !< Pressure (Pa) + real, dimension(:), intent(out) :: drho_dS_dS !< Partial derivative of beta with respect to S + real, dimension(:), intent(out) :: drho_dS_dT !< Partial derivative of beta with resepct to T + real, dimension(:), intent(out) :: drho_dT_dT !< Partial derivative of alpha with respect to T + real, dimension(:), intent(out) :: drho_dS_dP !< Partial derivative of beta with respect to pressure + real, dimension(:), intent(out) :: drho_dT_dP !< Partial derivative of alpha with respect to pressure + integer, intent(in) :: start !< Starting index within the array + integer, intent(in) :: npts !< The number of values to calculate + type(EOS_type), pointer :: EOS !< Equation of state structure + !! + if (.not.associated(EOS)) call MOM_error(FATAL, & + "calculate_density_derivs called with an unassociated EOS_type EOS.") + + select case (EOS%form_of_EOS) + case (EOS_LINEAR) + call calculate_density_second_derivs_linear(T, S, pressure, drho_dS_dS, drho_dS_dT, drho_dT_dT, drho_dS_dP, & + drho_dT_dP, start, npts) + case (EOS_WRIGHT) + call calculate_density_second_derivs_wright(T, S, pressure, drho_dS_dS, drho_dS_dT, drho_dT_dT, drho_dS_dP, & + drho_dT_dP, start, npts) + case (EOS_TEOS10) + call calculate_density_second_derivs_teos10(T, S, pressure, drho_dS_dS, drho_dS_dT, drho_dT_dT, drho_dS_dP, & + drho_dT_dP, start, npts) case default call MOM_error(FATAL, & "calculate_density_derivs: EOS%form_of_EOS is not valid.") end select -end subroutine calculate_density_derivs +end subroutine calculate_density_second_derivs_array + +!> Calls the appropriate subroutine to calculate density second derivatives for scalar nputs. +subroutine calculate_density_second_derivs_scalar(T, S, pressure, drho_dS_dS, drho_dS_dT, drho_dT_dT, drho_dS_dP, drho_dT_dP, & + EOS) + real, intent(in) :: T !< Potential temperature referenced to the surface (degC) + real, intent(in) :: S !< Salinity (PSU) + real, intent(in) :: pressure !< Pressure (Pa) + real, intent(out) :: drho_dS_dS !< Partial derivative of beta with respect to S + real, intent(out) :: drho_dS_dT !< Partial derivative of beta with resepct to T + real, intent(out) :: drho_dT_dT !< Partial derivative of alpha with respect to T + real, intent(out) :: drho_dS_dP !< Partial derivative of beta with respect to pressure + real, intent(out) :: drho_dT_dP !< Partial derivative of alpha with respect to pressure + type(EOS_type), pointer :: EOS !< Equation of state structure + !! + if (.not.associated(EOS)) call MOM_error(FATAL, & + "calculate_density_derivs called with an unassociated EOS_type EOS.") + + select case (EOS%form_of_EOS) + case (EOS_LINEAR) + call calculate_density_second_derivs_linear(T, S, pressure, drho_dS_dS, drho_dS_dT, drho_dT_dT, drho_dS_dP, & + drho_dT_dP) + case (EOS_WRIGHT) + call calculate_density_second_derivs_wright(T, S, pressure, drho_dS_dS, drho_dS_dT, drho_dT_dT, drho_dS_dP, & + drho_dT_dP) + case (EOS_TEOS10) + call calculate_density_second_derivs_teos10(T, S, pressure, drho_dS_dS, drho_dS_dT, drho_dT_dT, drho_dS_dP, & + drho_dT_dP) + case default + call MOM_error(FATAL, & + "calculate_density_derivs: EOS%form_of_EOS is not valid.") + end select + +end subroutine calculate_density_second_derivs_scalar !> Calls the appropriate subroutine to calculate specific volume derivatives for an array. subroutine calculate_specific_vol_derivs(T, S, pressure, dSV_dT, dSV_dS, start, npts, EOS) @@ -545,7 +654,7 @@ subroutine EOS_init(param_file, EOS) units="deg C Pa-1", default=0.0) endif - if (EOS%form_of_EOS == EOS_TEOS10 .OR. EOS%form_of_EOS == EOS_NEMO .AND. EOS%form_of_TFreeze /= TFREEZE_TEOS10) then + if ((EOS%form_of_EOS == EOS_TEOS10 .OR. EOS%form_of_EOS == EOS_NEMO) .AND. EOS%form_of_TFreeze /= TFREEZE_TEOS10) then call MOM_error(FATAL, "interpret_eos_selection: EOS_TEOS10 or EOS_NEMO \n" //& "should only be used along with TFREEZE_FORM = TFREEZE_TEOS10 .") endif @@ -553,6 +662,34 @@ subroutine EOS_init(param_file, EOS) end subroutine EOS_init +!> Manually initialized an EOS type (intended for unit testing of routines which need a specific EOS) +subroutine EOS_manual_init(EOS, form_of_EOS, form_of_TFreeze, EOS_quadrature, Compressible, Rho_T0_S0, drho_dT, & + dRho_dS, TFr_S0_P0, dTFr_dS, dTFr_dp) + type(EOS_type), pointer :: EOS + integer, optional, intent(in ) :: form_of_EOS + integer, optional, intent(in ) :: form_of_TFreeze + logical, optional, intent(in ) :: EOS_quadrature + logical, optional, intent(in ) :: Compressible + real , optional, intent(in ) :: Rho_T0_S0 + real , optional, intent(in ) :: drho_dT + real , optional, intent(in ) :: dRho_dS + real , optional, intent(in ) :: TFr_S0_P0 + real , optional, intent(in ) :: dTFr_dS + real , optional, intent(in ) :: dTFr_dp + + if (present(form_of_EOS )) EOS%form_of_EOS = form_of_EOS + if (present(form_of_TFreeze)) EOS%form_of_TFreeze = form_of_TFreeze + if (present(EOS_quadrature )) EOS%EOS_quadrature = EOS_quadrature + if (present(Compressible )) EOS%Compressible = Compressible + if (present(Rho_T0_S0 )) EOS%Rho_T0_S0 = Rho_T0_S0 + if (present(drho_dT )) EOS%drho_dT = drho_dT + if (present(dRho_dS )) EOS%dRho_dS = dRho_dS + if (present(TFr_S0_P0 )) EOS%TFr_S0_P0 = TFr_S0_P0 + if (present(dTFr_dS )) EOS%dTFr_dS = dTFr_dS + if (present(dTFr_dp )) EOS%dTFr_dp = dTFr_dp + +end subroutine EOS_manual_init + !> Allocates EOS_type subroutine EOS_allocate(EOS) type(EOS_type), pointer :: EOS !< Equation of state structure @@ -2161,7 +2298,33 @@ subroutine convert_temp_salt_for_TEOS10(T, S, press, G, kd, mask_z, EOS) enddo; enddo; enddo end subroutine convert_temp_salt_for_TEOS10 - +! Extractor routine for the EOS type if the members need to be accessed outside this module +subroutine extract_member_EOS(EOS, form_of_EOS, form_of_TFreeze, EOS_quadrature, Compressible, Rho_T0_S0, drho_dT, & + dRho_dS, TFr_S0_P0, dTFr_dS, dTFr_dp) + type(EOS_type), pointer :: EOS + integer, optional, intent(out) :: form_of_EOS + integer, optional, intent(out) :: form_of_TFreeze + logical, optional, intent(out) :: EOS_quadrature + logical, optional, intent(out) :: Compressible + real , optional, intent(out) :: Rho_T0_S0 + real , optional, intent(out) :: drho_dT + real , optional, intent(out) :: dRho_dS + real , optional, intent(out) :: TFr_S0_P0 + real , optional, intent(out) :: dTFr_dS + real , optional, intent(out) :: dTFr_dp + + if (present(form_of_EOS )) form_of_EOS = EOS%form_of_EOS + if (present(form_of_TFreeze)) form_of_TFreeze = EOS%form_of_TFreeze + if (present(EOS_quadrature )) EOS_quadrature = EOS%EOS_quadrature + if (present(Compressible )) Compressible = EOS%Compressible + if (present(Rho_T0_S0 )) Rho_T0_S0 = EOS%Rho_T0_S0 + if (present(drho_dT )) drho_dT = EOS%drho_dT + if (present(dRho_dS )) dRho_dS = EOS%dRho_dS + if (present(TFr_S0_P0 )) TFr_S0_P0 = EOS%TFr_S0_P0 + if (present(dTFr_dS )) dTFr_dS = EOS%dTFr_dS + if (present(dTFr_dp )) dTFr_dp = EOS%dTFr_dp + +end subroutine extract_member_EOS end module MOM_EOS diff --git a/src/equation_of_state/MOM_EOS_NEMO.F90 b/src/equation_of_state/MOM_EOS_NEMO.F90 index 3c6bd3973a..4755b4029c 100644 --- a/src/equation_of_state/MOM_EOS_NEMO.F90 +++ b/src/equation_of_state/MOM_EOS_NEMO.F90 @@ -42,6 +42,10 @@ module MOM_EOS_NEMO module procedure calculate_density_scalar_nemo, calculate_density_array_nemo end interface calculate_density_nemo +interface calculate_density_derivs_nemo + module procedure calculate_density_derivs_scalar_nemo, calculate_density_derivs_array_nemo +end interface calculate_density_derivs_nemo + real, parameter :: Pa2db = 1.e-4 real, parameter :: rdeltaS = 32. real, parameter :: r1_S0 = 0.875/35.16504 @@ -274,7 +278,7 @@ subroutine calculate_density_array_nemo(T, S, pressure, rho, start, npts) enddo end subroutine calculate_density_array_nemo -subroutine calculate_density_derivs_nemo(T, S, pressure, drho_dT, drho_dS, start, npts) +subroutine calculate_density_derivs_array_nemo(T, S, pressure, drho_dT, drho_dS, start, npts) real, intent(in), dimension(:) :: T !< Conservative temperature in C. real, intent(in), dimension(:) :: S !< Absolute salinity in g/kg. real, intent(in), dimension(:) :: pressure !< Pressure in Pa. @@ -352,7 +356,27 @@ subroutine calculate_density_derivs_nemo(T, S, pressure, drho_dT, drho_dS, start drho_dS(j) = zn / zs enddo -end subroutine calculate_density_derivs_nemo +end subroutine calculate_density_derivs_array_nemo + +!> Wrapper to calculate_density_derivs_array for scalar inputs +subroutine calculate_density_derivs_scalar_nemo(T, S, pressure, drho_dt, drho_ds) + real, intent(in) :: T, S, pressure + real, intent(out) :: drho_dt + real, intent(out) :: drho_ds + ! Local variables + real :: al0, p0, lambda + integer :: j + real, dimension(1) :: T0, S0, pressure0 + real, dimension(1) :: drdt0, drds0 + + T0(1) = T + S0(1) = S + pressure0(1) = pressure + + call calculate_density_derivs_array_nemo(T0, S0, pressure0, drdt0, drds0, 1, 1) + drho_dt = drdt0(1) + drho_ds = drds0(1) +end subroutine calculate_density_derivs_scalar_nemo subroutine calculate_compress_nemo(T, S, pressure, rho, drho_dp, start, npts) real, intent(in), dimension(:) :: T !< Conservative temperature in C. diff --git a/src/equation_of_state/MOM_EOS_TEOS10.F90 b/src/equation_of_state/MOM_EOS_TEOS10.F90 index 6dc7af1706..f573b6e3ee 100644 --- a/src/equation_of_state/MOM_EOS_TEOS10.F90 +++ b/src/equation_of_state/MOM_EOS_TEOS10.F90 @@ -27,19 +27,29 @@ module MOM_EOS_TEOS10 use gsw_mod_toolbox, only : gsw_sp_from_sr, gsw_pt_from_ct use gsw_mod_toolbox, only : gsw_rho, gsw_rho_first_derivatives, gsw_specvol_first_derivatives +use gsw_mod_toolbox, only : gsw_rho_second_derivatives !use gsw_mod_toolbox, only : gsw_sr_from_sp, gsw_ct_from_pt implicit none ; private public calculate_compress_teos10, calculate_density_teos10 -public calculate_density_derivs_teos10, calculate_specvol_derivs_teos10 -public calculate_density_scalar_teos10, calculate_density_array_teos10 +public calculate_density_derivs_teos10 +public calculate_specvol_derivs_teos10 +public calculate_density_second_derivs_teos10 public gsw_sp_from_sr, gsw_pt_from_ct interface calculate_density_teos10 module procedure calculate_density_scalar_teos10, calculate_density_array_teos10 end interface calculate_density_teos10 +interface calculate_density_derivs_teos10 + module procedure calculate_density_derivs_scalar_teos10, calculate_density_derivs_array_teos10 +end interface calculate_density_derivs_teos10 + +interface calculate_density_second_derivs_teos10 + module procedure calculate_density_second_derivs_scalar_teos10, calculate_density_second_derivs_array_teos10 +end interface calculate_density_second_derivs_teos10 + real, parameter :: Pa2db = 1.e-4 ! The conversion factor from Pa to dbar. contains @@ -110,7 +120,7 @@ subroutine calculate_density_array_teos10(T, S, pressure, rho, start, npts) enddo end subroutine calculate_density_array_teos10 -subroutine calculate_density_derivs_teos10(T, S, pressure, drho_dT, drho_dS, start, npts) +subroutine calculate_density_derivs_array_teos10(T, S, pressure, drho_dT, drho_dS, start, npts) real, intent(in), dimension(:) :: T !< Conservative temperature in C. real, intent(in), dimension(:) :: S !< Absolute salinity in g/kg. real, intent(in), dimension(:) :: pressure !< Pressure in Pa. @@ -141,7 +151,20 @@ subroutine calculate_density_derivs_teos10(T, S, pressure, drho_dT, drho_dS, sta call gsw_rho_first_derivatives(zs, zt, zp, drho_dsa=drho_dS(j), drho_dct=drho_dT(j)) enddo -end subroutine calculate_density_derivs_teos10 +end subroutine calculate_density_derivs_array_teos10 + +subroutine calculate_density_derivs_scalar_teos10(T, S, pressure, drho_dT, drho_dS) + real, intent(in) :: T, S, pressure + real, intent(out) :: drho_dT, drho_dS + ! Local variables + real :: zs,zt,zp + !Conversions + zs = S !gsw_sr_from_sp(S) !Convert practical salinity to absolute salinity + zt = T !gsw_ct_from_pt(S,T) !Convert potantial temp to conservative temp + zp = pressure* Pa2db !Convert pressure from Pascal to decibar + if(S.lt.-1.0e-10) return !Can we assume safely that this is a missing value? + call gsw_rho_first_derivatives(zs, zt, zp, drho_dsa=drho_dS, drho_dct=drho_dT) +end subroutine calculate_density_derivs_scalar_teos10 subroutine calculate_specvol_derivs_teos10(T, S, pressure, dSV_dT, dSV_dS, start, npts) real, intent(in), dimension(:) :: T !< Conservative temperature in C. @@ -176,6 +199,66 @@ subroutine calculate_specvol_derivs_teos10(T, S, pressure, dSV_dT, dSV_dS, start end subroutine calculate_specvol_derivs_teos10 +!> Calculate the 5 second derivatives of the equation of state for scalar inputs +subroutine calculate_density_second_derivs_scalar_teos10(T, S, pressure, drho_dS_dS, drho_dS_dT, drho_dT_dT, drho_dS_dP, & + drho_dT_dP) + real, intent(in) :: T, S, pressure + real, intent(out) :: drho_dS_dS !< Partial derivative of beta with respect to S + real, intent(out) :: drho_dS_dT !< Partial derivative of beta with resepct to T + real, intent(out) :: drho_dT_dT !< Partial derivative of alpha with respect to T + real, intent(out) :: drho_dS_dP !< Partial derivative of beta with respect to pressure + real, intent(out) :: drho_dT_dP !< Partial derivative of alpha with respect to pressure +! * Arguments: T - conservative temperature in C. * +! * (in) S - absolute salinity in g/kg. * +! * (in) pressure - pressure in Pa. * +! * (out) drho_dT - the partial derivative of density with * +! * potential temperature, in kg m-3 K-1. * +! * (out) drho_dS - the partial derivative of density with * +! * salinity, in kg m-3 psu-1. * + real :: zs,zt,zp + + !Conversions + zs = S !gsw_sr_from_sp(S) !Convert practical salinity to absolute salinity + zt = T !gsw_ct_from_pt(S,T) !Convert potantial temp to conservative temp + zp = pressure* Pa2db !Convert pressure from Pascal to decibar + if(S.lt.-1.0e-10) return !Can we assume safely that this is a missing value? + call gsw_rho_second_derivatives(zs, zt, zp, rho_sa_sa=drho_dS_dS, rho_sa_ct=drho_dS_dT, & + rho_ct_ct=drho_dT_dT, rho_sa_p=drho_dS_dP, rho_ct_p=drho_dT_dP) + +end subroutine calculate_density_second_derivs_scalar_teos10 + +!> Calculate the 5 second derivatives of the equation of state for scalar inputs +subroutine calculate_density_second_derivs_array_teos10(T, S, pressure, drho_dS_dS, drho_dS_dT, drho_dT_dT, drho_dS_dP, & + drho_dT_dP, start, npts) + real, dimension(:), intent(in) :: T, S, pressure + real, dimension(:), intent(out) :: drho_dS_dS !< Partial derivative of beta with respect to S + real, dimension(:), intent(out) :: drho_dS_dT !< Partial derivative of beta with resepct to T + real, dimension(:), intent(out) :: drho_dT_dT !< Partial derivative of alpha with respect to T + real, dimension(:), intent(out) :: drho_dS_dP !< Partial derivative of beta with respect to pressure + real, dimension(:), intent(out) :: drho_dT_dP !< Partial derivative of alpha with respect to pressure + integer, intent(in) :: start !< The starting point in the arrays. + integer, intent(in) :: npts !< The number of values to calculate. +! * Arguments: T - conservative temperature in C. * +! * (in) S - absolute salinity in g/kg. * +! * (in) pressure - pressure in Pa. * +! * (out) drho_dT - the partial derivative of density with * +! * potential temperature, in kg m-3 K-1. * +! * (out) drho_dS - the partial derivative of density with * +! * salinity, in kg m-3 psu-1. * + real :: zs,zt,zp + integer :: j + do j=start,start+npts-1 + !Conversions + zs = S(j) !gsw_sr_from_sp(S) !Convert practical salinity to absolute salinity + zt = T(j) !gsw_ct_from_pt(S,T) !Convert potantial temp to conservative temp + zp = pressure(j)* Pa2db !Convert pressure from Pascal to decibar + if(zs .lt. -1.0e-10) return !Can we assume safely that this is a missing value? + call gsw_rho_second_derivatives(zs, zt, zp, rho_sa_sa=drho_dS_dS(j), rho_sa_ct=drho_dS_dT(j), & + rho_ct_ct=drho_dT_dT(j), rho_sa_p=drho_dS_dP(j), rho_ct_p=drho_dT_dP(j)) + enddo + +end subroutine calculate_density_second_derivs_array_teos10 + !> This subroutine computes the in situ density of sea water (rho in * !! units of kg/m^3) and the compressibility (drho/dp = C_sound^-2) * !! (drho_dp in units of s2 m-2) from salinity (sal in psu), potential* diff --git a/src/equation_of_state/MOM_EOS_Wright.F90 b/src/equation_of_state/MOM_EOS_Wright.F90 index b0354b0815..e868b5f81f 100644 --- a/src/equation_of_state/MOM_EOS_Wright.F90 +++ b/src/equation_of_state/MOM_EOS_Wright.F90 @@ -33,13 +33,21 @@ module MOM_EOS_Wright public calculate_compress_wright, calculate_density_wright public calculate_density_derivs_wright, calculate_specvol_derivs_wright -public calculate_density_scalar_wright, calculate_density_array_wright +public calculate_density_second_derivs_wright public int_density_dz_wright, int_spec_vol_dp_wright interface calculate_density_wright module procedure calculate_density_scalar_wright, calculate_density_array_wright end interface calculate_density_wright +interface calculate_density_derivs_wright + module procedure calculate_density_derivs_scalar_wright, calculate_density_derivs_array_wright +end interface + +interface calculate_density_second_derivs_wright + module procedure calculate_density_second_derivs_scalar_wright, calculate_density_second_derivs_array_wright +end interface + !real :: a0, a1, a2, b0, b1, b2, b3, b4, b5, c0, c1, c2, c3, c4, c5 ! One of the two following blocks of values should be commented out. ! Following are the values for the full range formula. @@ -140,8 +148,8 @@ subroutine calculate_density_array_wright(T, S, pressure, rho, start, npts) enddo end subroutine calculate_density_array_wright -! #@# This subroutine needs a doxygen description. -subroutine calculate_density_derivs_wright(T, S, pressure, drho_dT, drho_dS, start, npts) +!> For a given thermodynamic state, return the thermal/haline expansion coefficients +subroutine calculate_density_derivs_array_wright(T, S, pressure, drho_dT, drho_dS, start, npts) real, intent(in), dimension(:) :: T !< Potential temperature relative to the surface !! in C. real, intent(in), dimension(:) :: S !< Salinity in PSU. @@ -180,7 +188,105 @@ subroutine calculate_density_derivs_wright(T, S, pressure, drho_dT, drho_dS, sta (pressure(j)+p0) * ( (pressure(j)+p0)*a2 + (c4 + c5*T(j)) )) enddo -end subroutine calculate_density_derivs_wright +end subroutine calculate_density_derivs_array_wright + +!> The scalar version of calculate_density_derivs which promotes scalar inputs to a 1-element array and then +!! demotes the output back to a scalar +subroutine calculate_density_derivs_scalar_wright(T, S, pressure, drho_dT, drho_dS) + real, intent(in) :: T !< Potential temperature relative to the surface + !! in C. + real, intent(in) :: S !< Salinity in PSU. + real, intent(in) :: pressure !< Pressure in Pa. + real, intent(out) :: drho_dT !< The partial derivative of density with potential + !! temperature, in kg m-3 K-1. + real, intent(out) :: drho_dS !< The partial derivative of density with salinity, + !! in kg m-3 psu-1. + + ! Local variables needed to promote the input/output scalars to 1-element arrays + real, dimension(1) :: T0, S0, P0 + real, dimension(1) :: drdt0, drds0 + + T0(1) = T + S0(1) = S + P0(1) = pressure + call calculate_density_derivs_array_wright(T0, S0, P0, drdt0, drds0, 1, 1) + drho_dT = drdt0(1) + drho_dS = drds0(1) + +end subroutine calculate_density_derivs_scalar_wright + +!> Second derivatives of density with respect to temperature, salinity, and pressure +subroutine calculate_density_second_derivs_array_wright(T, S, P, drho_ds_ds, drho_ds_dt, drho_dt_dt, & + drho_ds_dp, drho_dt_dp, start, npts) + real, dimension(:), intent(in ) :: T !< Potential temperature referenced to 0 dbar + real, dimension(:), intent(in ) :: S !< Salinity in PSU + real, dimension(:), intent(in ) :: P !< Pressure in Pa + real, dimension(:), intent( out) :: drho_ds_ds !< Partial derivative of beta with respect to S + real, dimension(:), intent( out) :: drho_ds_dt !< Partial derivative of beta with resepct to T + real, dimension(:), intent( out) :: drho_dt_dt !< Partial derivative of alpha with respect to T + real, dimension(:), intent( out) :: drho_ds_dp !< Partial derivative of beta with respect to pressure + real, dimension(:), intent( out) :: drho_dt_dp !< Partial derivative of alpha with respect to pressure + integer, intent(in ) :: start !< Starting index in T,S,P + integer, intent(in ) :: npts !< Number of points to loop over + + integer :: j + ! Based on the above expression with common terms factored, there probably exists a more numerically stable + ! and/or efficient expression + real :: z0, z1, z2, z3, z4, z5, z6 ,z7, z8, z9, z10, z11, z2_2, z2_3 + + do j = start,start+npts-1 + z0 = T(j)*(b1 + b5*S(j) + T(j)*(b2 + b3*T(j))) + z1 = (b0 + P(j) + b4*S(j) + z0) + z3 = (b1 + b5*S(j) + T(j)*(2.*b2 + 2.*b3*T(j))) + z4 = (c0 + c4*S(j) + T(j)*(c1 + c5*S(j) + T(j)*(c2 + c3*T(j)))) + z5 = (b1 + b5*S(j) + T(j)*(b2 + b3*T(j)) + T(j)*(b2 + 2.*b3*T(j))) + z6 = c1 + c5*S(j) + T(j)*(c2 + c3*T(j)) + T(j)*(c2 + 2.*c3*T(j)) + z7 = (c4 + c5*T(j) + a2*z1) + z8 = (c1 + c5*S(j) + T(j)*(2.*c2 + 3.*c3*T(j)) + a1*z1) + z9 = (a0 + a2*S(j) + a1*T(j)) + z10 = (b4 + b5*T(j)) + z11 = (z10*z4 - z1*z7) + z2 = (c0 + c4*S(j) + T(j)*(c1 + c5*S(j) + T(j)*(c2 + c3*T(j))) + z9*z1) + z2_2 = z2*z2 + z2_3 = z2_2*z2 + + drho_ds_ds(j) = (z10*(c4 + c5*T(j)) - a2*z10*z1 - z10*z7)/z2_2 - (2.*(c4 + c5*T(j) + z9*z10 + a2*z1)*z11)/z2_3 + drho_ds_dt(j) = (z10*z6 - z1*(c5 + a2*z5) + b5*z4 - z5*z7)/z2_2 - (2.*(z6 + z9*z5 + a1*z1)*z11)/z2_3 + drho_dt_dt(j) = (z3*z6 - z1*(2.*c2 + 6.*c3*T(j) + a1*z5) + (2.*b2 + 4.*b3*T(j))*z4 - z5*z8)/z2_2 - & + (2.*(z6 + z9*z5 + a1*z1)*(z3*z4 - z1*z8))/z2_3 + drho_ds_dp(j) = (-c4 - c5*T(j) - 2.*a2*z1)/z2_2 - (2.*z9*z11)/z2_3 + drho_dt_dp(j) = (-c1 - c5*S(j) - T(j)*(2.*c2 + 3.*c3*T(j)) - 2.*a1*z1)/z2_2 - (2.*z9*(z3*z4 - z1*z8))/z2_3 + enddo + +end subroutine calculate_density_second_derivs_array_wright + +!> Second derivatives of density with respect to temperature, salinity, and pressure for scalar inputs. Inputs +!! promoted to 1-element array and output demoted to scalar +subroutine calculate_density_second_derivs_scalar_wright(T, S, P, drho_ds_ds, drho_ds_dt, drho_dt_dt, & + drho_ds_dp, drho_dt_dp) + real, intent(in ) :: T !< Potential temperature referenced to 0 dbar + real, intent(in ) :: S !< Salinity in PSU + real, intent(in ) :: P !< Pressure in Pa + real, intent( out) :: drho_ds_ds !< Partial derivative of beta with respect to S + real, intent( out) :: drho_ds_dt !< Partial derivative of beta with resepct to T + real, intent( out) :: drho_dt_dt !< Partial derivative of alpha with respect to T + real, intent( out) :: drho_ds_dp !< Partial derivative of beta with respect to pressure + real, intent( out) :: drho_dt_dp !< Partial derivative of alpha with respect to pressure + ! Local variables + real, dimension(1) :: T0, S0, P0 + real, dimension(1) :: drdsds, drdsdt, drdtdt, drdsdp, drdtdp + + T0(1) = T + S0(1) = S + P0(1) = P + call calculate_density_second_derivs_array_wright(T0, S0, P0, drdsds, drdsdt, drdtdt, drdsdp, drdtdp, 1, 1) + drho_ds_ds = drdsds(1) + drho_ds_dt = drdsdt(1) + drho_dt_dt = drdtdt(1) + drho_ds_dp = drdsdp(1) + drho_dt_dp = drdtdp(1) + +end subroutine calculate_density_second_derivs_scalar_wright subroutine calculate_specvol_derivs_wright(T, S, pressure, dSV_dT, dSV_dS, start, npts) real, intent(in), dimension(:) :: T !< Potential temperature relative to the surface diff --git a/src/equation_of_state/MOM_EOS_linear.F90 b/src/equation_of_state/MOM_EOS_linear.F90 index bb3acbf56a..0e51e336b3 100644 --- a/src/equation_of_state/MOM_EOS_linear.F90 +++ b/src/equation_of_state/MOM_EOS_linear.F90 @@ -32,14 +32,24 @@ module MOM_EOS_linear #include public calculate_compress_linear, calculate_density_linear -public calculate_density_derivs_linear, calculate_specvol_derivs_linear +public calculate_density_derivs_linear, calculate_density_derivs_scalar_linear +public calculate_specvol_derivs_linear public calculate_density_scalar_linear, calculate_density_array_linear +public calculate_density_second_derivs_linear public int_density_dz_linear, int_spec_vol_dp_linear interface calculate_density_linear module procedure calculate_density_scalar_linear, calculate_density_array_linear end interface calculate_density_linear +interface calculate_density_derivs_linear + module procedure calculate_density_derivs_scalar_linear, calculate_density_derivs_array_linear +end interface calculate_density_derivs_linear + +interface calculate_density_second_derivs_linear + module procedure calculate_density_second_derivs_scalar_linear, calculate_density_second_derivs_array_linear +end interface calculate_density_second_derivs_linear + contains !> This subroutine computes the density of sea water with a trivial @@ -115,8 +125,8 @@ end subroutine calculate_density_array_linear !> This subroutine calculates the partial derivatives of density * !! with potential temperature and salinity. -subroutine calculate_density_derivs_linear(T, S, pressure, drho_dT_out, & - drho_dS_out, start, npts, Rho_T0_S0, dRho_dT, dRho_dS) +subroutine calculate_density_derivs_array_linear(T, S, pressure, drho_dT_out, & + drho_dS_out, Rho_T0_S0, dRho_dT, dRho_dS, start, npts) real, intent(in), dimension(:) :: T !< Potential temperature relative to the surface !! in C. real, intent(in), dimension(:) :: S !< Salinity in PSU. @@ -125,12 +135,12 @@ subroutine calculate_density_derivs_linear(T, S, pressure, drho_dT_out, & !! potential temperature, in kg m-3 K-1. real, intent(out), dimension(:) :: drho_dS_out !< The partial derivative of density with !! salinity, in kg m-3 psu-1. - integer, intent(in) :: start !< The starting point in the arrays. - integer, intent(in) :: npts !< The number of values to calculate. real, intent(in) :: Rho_T0_S0 !< The density at T=0, S=0, in kg m-3. real, intent(in) :: dRho_dT, dRho_dS !< The derivatives of density with !! temperature and salinity, in kg m-3 C-1 !! and kg m-3 psu-1. + integer, intent(in) :: start !< The starting point in the arrays. + integer, intent(in) :: npts !< The number of values to calculate. ! * This subroutine calculates the partial derivatives of density * ! * with potential temperature and salinity. * @@ -154,7 +164,76 @@ subroutine calculate_density_derivs_linear(T, S, pressure, drho_dT_out, & drho_dS_out(j) = dRho_dS enddo -end subroutine calculate_density_derivs_linear +end subroutine calculate_density_derivs_array_linear + +!> This subroutine calculates the partial derivatives of density * +!! with potential temperature and salinity for a single point. +subroutine calculate_density_derivs_scalar_linear(T, S, pressure, drho_dT_out, & + drho_dS_out, Rho_T0_S0, dRho_dT, dRho_dS) + real, intent(in) :: T !< Potential temperature relative to the surface + !! in C. + real, intent(in) :: S !< Salinity in PSU. + real, intent(in) :: pressure !< Pressure in Pa. + real, intent(out) :: drho_dT_out !< The partial derivative of density with + !! potential temperature, in kg m-3 K-1. + real, intent(out) :: drho_dS_out !< The partial derivative of density with + !! salinity, in kg m-3 psu-1. + real, intent(in) :: Rho_T0_S0 !< The density at T=0, S=0, in kg m-3. + real, intent(in) :: dRho_dT, dRho_dS !< The derivatives of density with + !! temperature and salinity, in kg m-3 C-1 + !! and kg m-3 psu-1. + drho_dT_out = dRho_dT + drho_dS_out = dRho_dS + +end subroutine calculate_density_derivs_scalar_linear + +!> This subroutine calculates the five, partial second derivatives of density w.r.t. +!! potential temperature and salinity and pressure which for a linear equation of state should all be 0. +subroutine calculate_density_second_derivs_scalar_linear(T, S,pressure, drho_dS_dS, drho_dS_dT, drho_dT_dT,& + drho_dS_dP, drho_dT_dP) + real, intent(in) :: T !< Potential temperature relative to the surface + !! in C. + real, intent(in) :: S !< Salinity in PSU. + real, intent(in) :: pressure !< Pressure in Pa. + real, intent(out) :: drho_dS_dS !< The partial derivative of density with + real, intent(out) :: drho_dS_dT !< The partial derivative of density with + real, intent(out) :: drho_dT_dT !< The partial derivative of density with + real, intent(out) :: drho_dS_dP !< The partial derivative of density with + real, intent(out) :: drho_dT_dP !< The partial derivative of density with + + drho_dS_dS = 0. + drho_dS_dT = 0. + drho_dT_dT = 0. + drho_dS_dP = 0. + drho_dT_dP = 0. + +end subroutine calculate_density_second_derivs_scalar_linear + +!> This subroutine calculates the five, partial second derivatives of density w.r.t. +!! potential temperature and salinity and pressure which for a linear equation of state should all be 0. +subroutine calculate_density_second_derivs_array_linear(T, S,pressure, drho_dS_dS, drho_dS_dT, drho_dT_dT,& + drho_dS_dP, drho_dT_dP, start, npts) + real, dimension(:), intent(in) :: T !< Potential temperature relative to the surface + !! in C. + real, dimension(:), intent(in) :: S !< Salinity in PSU. + real, dimension(:), intent(in) :: pressure !< Pressure in Pa. + real, dimension(:), intent(out) :: drho_dS_dS !< The partial derivative of density with + real, dimension(:), intent(out) :: drho_dS_dT !< The partial derivative of density with + real, dimension(:), intent(out) :: drho_dT_dT !< The partial derivative of density with + real, dimension(:), intent(out) :: drho_dS_dP !< The partial derivative of density with + real, dimension(:), intent(out) :: drho_dT_dP !< The partial derivative of density with + integer, intent(in) :: start !< The starting point in the arrays. + integer, intent(in) :: npts !< The number of values to calculate. + integer :: j + do j=start,start+npts-1 + drho_dS_dS(j) = 0. + drho_dS_dT(j) = 0. + drho_dT_dT(j) = 0. + drho_dS_dP(j) = 0. + drho_dT_dP(j) = 0. + enddo + +end subroutine calculate_density_second_derivs_array_linear ! #@# This subroutine needs a doxygen description. subroutine calculate_specvol_derivs_linear(T, S, pressure, dSV_dT, dSV_dS, & diff --git a/src/equation_of_state/TEOS10/gsw_rho_second_derivatives.f90 b/src/equation_of_state/TEOS10/gsw_rho_second_derivatives.f90 new file mode 100644 index 0000000000..fdf75e7a0a --- /dev/null +++ b/src/equation_of_state/TEOS10/gsw_rho_second_derivatives.f90 @@ -0,0 +1,78 @@ +!========================================================================== +elemental subroutine gsw_rho_second_derivatives (sa, ct, p, rho_sa_sa, & + rho_sa_ct, rho_ct_ct, rho_sa_p, rho_ct_p) +!========================================================================== +! +! Calculates five second-order derivatives of rho. Note that this function +! uses the using the computationally-efficient expression for specific +! volume (Roquet et al., 2014). +! +! SA = Absolute Salinity [ g/kg ] +! CT = Conservative Temperature (ITS-90) [ deg C ] +! p = sea pressure [ dbar ] +! ( i.e. absolute pressure - 10.1325 dbar ) +! +! rho_SA_SA = The second-order derivative of rho with respect to +! Absolute Salinity at constant CT & p. [ J/(kg (g/kg)^2) ] +! rho_SA_CT = The second-order derivative of rho with respect to +! SA and CT at constant p. [ J/(kg K(g/kg)) ] +! rho_CT_CT = The second-order derivative of rho with respect to CT at +! constant SA & p +! rho_SA_P = The second-order derivative with respect to SA & P at +! constant CT. +! rho_CT_P = The second-order derivative with respect to CT & P at +! constant SA. +!-------------------------------------------------------------------------- + +use gsw_mod_toolbox, only : gsw_specvol, gsw_specvol_first_derivatives +use gsw_mod_toolbox, only : gsw_specvol_second_derivatives + +use gsw_mod_kinds + +implicit none + +real (r8), intent(in) :: sa, ct, p +real (r8), intent(out), optional :: rho_sa_sa, rho_sa_ct, rho_ct_ct +real (r8), intent(out), optional :: rho_sa_p, rho_ct_p + +integer :: iflag1, iflag2 +real (r8) :: rec_v, rec_v2, rec_v3, v_ct, v_ct_ct, v_ct_p, v_p, v_sa, v_sa_ct +real (r8) :: v_sa_p, v_sa_sa + +iflag1 = 0 +if (present(rho_sa_sa) .or. present(rho_sa_ct) & + .or. present(rho_sa_p)) iflag1 = ibset(iflag1,1) +if (present(rho_sa_ct) .or. present(rho_ct_ct) & + .or. present(rho_ct_p)) iflag1 = ibset(iflag1,2) +if (present(rho_sa_p) .or. present(rho_ct_p)) iflag1 = ibset(iflag1,3) + +call gsw_specvol_first_derivatives(sa,ct,p,v_sa,v_ct,v_p,iflag=iflag1) + +iflag2 = 0 +if (present(rho_sa_sa)) iflag2 = ibset(iflag2,1) +if (present(rho_sa_ct)) iflag2 = ibset(iflag2,2) +if (present(rho_ct_ct)) iflag2 = ibset(iflag2,3) +if (present(rho_sa_p)) iflag2 = ibset(iflag2,4) +if (present(rho_ct_p)) iflag2 = ibset(iflag2,5) + +call gsw_specvol_second_derivatives(sa,ct,p,v_sa_sa,v_sa_ct,v_ct_ct, & + v_sa_p,v_ct_p,iflag=iflag2) + +rec_v = 1.0_r8/gsw_specvol(sa,ct,p) +rec_v2 = rec_v**2 +rec_v3 = rec_v2*rec_v + +if (present(rho_sa_sa)) rho_sa_sa = -v_sa_sa*rec_v2 + 2.0_r8*v_sa*v_sa*rec_v3 + +if (present(rho_sa_ct)) rho_sa_ct = -v_sa_ct*rec_v2 + 2.0_r8*v_sa*v_ct*rec_v3 + +if (present(rho_ct_ct)) rho_ct_ct = -v_ct_ct*rec_v2 + 2.0_r8*v_ct*v_ct*rec_v3 + +if (present(rho_sa_p)) rho_sa_p = -v_sa_p*rec_v2 + 2.0_r8*v_sa*v_p*rec_v3 + +if (present(rho_ct_p)) rho_ct_p = -v_ct_p*rec_v2 + 2.0_r8*v_ct*v_p*rec_v3 + +return +end subroutine + +!-------------------------------------------------------------------------- diff --git a/src/equation_of_state/TEOS10/gsw_specvol_second_derivatives.f90 b/src/equation_of_state/TEOS10/gsw_specvol_second_derivatives.f90 new file mode 100644 index 0000000000..39096109e9 --- /dev/null +++ b/src/equation_of_state/TEOS10/gsw_specvol_second_derivatives.f90 @@ -0,0 +1,131 @@ +!========================================================================== +elemental subroutine gsw_specvol_second_derivatives (sa, ct, p, v_sa_sa, & + v_sa_ct, v_ct_ct, v_sa_p, v_ct_p, iflag) +! ========================================================================= +! +! Calculates five second-order derivatives of specific volume (v). +! Note that this function uses the computationally-efficient +! expression for specific volume (Roquet et al., 2014). +! +! SA = Absolute Salinity [ g/kg ] +! CT = Conservative Temperature (ITS-90) [ deg C ] +! p = sea pressure [ dbar ] +! ( i.e. absolute pressure - 10.1325 dbar ) +! +! v_SA_SA = The second derivative of specific volume with respect to +! Absolute Salinity at constant CT & p. [ J/(kg (g/kg)^2) ] +! v_SA_CT = The second derivative of specific volume with respect to +! SA and CT at constant p. [ J/(kg K(g/kg)) ] +! v_CT_CT = The second derivative of specific volume with respect to +! CT at constant SA and p. [ J/(kg K^2) ] +! v_SA_P = The second derivative of specific volume with respect to +! SA and P at constant CT. [ J/(kg K(g/kg)) ] +! v_CT_P = The second derivative of specific volume with respect to +! CT and P at constant SA. [ J/(kg K(g/kg)) ] +!-------------------------------------------------------------------------- + +use gsw_mod_teos10_constants, only : gsw_sfac, offset + +use gsw_mod_specvol_coefficients + +use gsw_mod_kinds + +implicit none + +real (r8), intent(in) :: sa, ct, p +integer, intent(in), optional :: iflag +real (r8), intent(out), optional :: v_sa_sa, v_sa_ct, v_ct_ct, v_sa_p, v_ct_p + +integer :: i +logical :: flags(5) +real (r8) :: v_ct_ct_part, v_ct_p_part, v_sa_ct_part, v_sa_p_part +real (r8) :: v_sa_sa_part, xs, xs2, ys, z + +xs2 = gsw_sfac*sa + offset +xs = sqrt(xs2) +ys = ct*0.025_r8 +z = p*1e-4_r8 + +if (present(iflag)) then + do i = 1, 5 + flags(i) = btest(iflag,i) + end do +else + flags = .true. +end if + +if (present(v_sa_sa) .and. flags(1)) then + + v_sa_sa_part = (-b000 + xs2*(b200 + xs*(2.0_r8*b300 + xs*(3.0_r8*b400 & + + 4.0_r8*b500*xs))) + ys*(-b010 + xs2*(b210 + xs*(2.0_r8*b310 & + + 3.0_r8*b410*xs)) + ys*(-b020 + xs2*(b220 + 2.0_r8*b320*xs) & + + ys*(-b030 + b230*xs2 + ys*(-b040 - b050*ys)))) + z*(-b001 & + + xs2*(b201 + xs*(2.0_r8*b301 + 3.0_r8*b401*xs)) + ys*(-b011 & + + xs2*(b211 + 2.0_r8*b311*xs) + ys*(-b021 + b221*xs2 & + + ys*(-b031 - b041*ys))) + z*(-b002 + xs2*(b202 + 2.0_r8*b302*xs) & + + ys*(-b012 + b212*xs2 + ys*(-b022 - b032*ys)) + z*(-b003 & + - b013*ys - b004*z))))/xs2 + + v_sa_sa = 0.25_r8*gsw_sfac*gsw_sfac*v_sa_sa_part/xs + +end if + +if (present(v_sa_ct) .and. flags(2)) then + + v_sa_ct_part = (b010 + xs*(b110 + xs*(b210 + xs*(b310 + b410*xs))) & + + ys*(2.0_r8*(b020 + xs*(b120 + xs*(b220 + b320*xs))) & + + ys*(3.0_r8*(b030 + xs*(b130 + b230*xs)) + ys*(4.0_r8*(b040 + b140*xs) & + + 5.0_r8*b050*ys))) + z*(b011 + xs*(b111 + xs*(b211 + b311*xs)) & + + ys*(2.0_r8*(b021 + xs*(b121 + b221*xs)) + ys*(3.0_r8*(b031 + b131*xs) & + + 4.0_r8*b041*ys)) + z*(b012 + xs*(b112 + b212*xs) + ys*(2.0_r8*(b022 & + + b122*xs) + 3.0_r8*b032*ys) + b013*z)))/xs + + v_sa_ct = 0.025_r8*0.5_r8*gsw_sfac*v_sa_ct_part + +end if + +if (present(v_ct_ct) .and. flags(3)) then + + v_ct_ct_part = a010 + xs*(a110 + xs*(a210 + xs*(a310 + a410*xs))) & + + ys*(2.0_r8*(a020 + xs*(a120 + xs*(a220 + a320*xs))) & + + ys*(3.0_r8*(a030 + xs*(a130 + a230*xs)) + ys*(4.0_r8*(a040 & + + a140*xs) + 5.0_r8*a050*ys))) + z*( a011 + xs*(a111 + xs*(a211 & + + a311*xs)) + ys*(2.0_r8*(a021 + xs*(a121 + a221*xs)) & + + ys*(3.0_r8*(a031 + a131*xs) + 4.0_r8*a041*ys)) + z*(a012 & + + xs*(a112 + a212*xs) + ys*(2.0_r8*(a022 + a122*xs) & + + 3.0_r8*a032*ys) + a013*z)) + + v_ct_ct = 0.025_r8*0.025_r8*v_ct_ct_part + +end if + +if (present(v_sa_p) .and. flags(4)) then + + v_sa_p_part = b001 + xs*(b101 + xs*(b201 + xs*(b301 + b401*xs))) + ys*(b011 & + + xs*(b111 + xs*(b211 + b311*xs)) + ys*(b021 + xs*(b121 + b221*xs) & + + ys*(b031 + b131*xs + b041*ys))) + z*(2.0_r8*(b002 + xs*(b102 & + + xs*(b202 + b302*xs)) + ys*(b012 + xs*(b112 + b212*xs) + ys*(b022 & + + b122*xs + b032*ys))) + z*(3.0_r8*(b003 + b103*xs + b013*ys) & + + 4.0_r8*b004*z)) + + v_sa_p = 1e-8_r8*0.5_r8*gsw_sfac*v_sa_p_part + +end if + +if (present(v_ct_p) .and. flags(5)) then + + v_ct_p_part = a001 + xs*(a101 + xs*(a201 + xs*(a301 + a401*xs))) + ys*(a011 & + + xs*(a111 + xs*(a211 + a311*xs)) + ys*(a021 + xs*(a121 + a221*xs) & + + ys*(a031 + a131*xs + a041*ys))) + z*(2.0_r8*(a002 + xs*(a102 & + + xs*(a202 + a302*xs)) + ys*(a012 + xs*(a112 + a212*xs) + ys*(a022 & + + a122*xs + a032*ys))) + z*(3.0_r8*(a003 + a103*xs + a013*ys) & + + 4.0_r8*a004*z)) + + v_ct_p = 1e-8_r8*0.025_r8*v_ct_p_part + +end if + +return +end subroutine + +!-------------------------------------------------------------------------- diff --git a/src/tracer/MOM_neutral_diffusion.F90 b/src/tracer/MOM_neutral_diffusion.F90 index 0093ebb953..f7501b3b95 100644 --- a/src/tracer/MOM_neutral_diffusion.F90 +++ b/src/tracer/MOM_neutral_diffusion.F90 @@ -3,17 +3,25 @@ module MOM_neutral_diffusion ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end -use MOM_cpu_clock, only : CLOCK_MODULE, CLOCK_ROUTINE -use MOM_diag_mediator, only : diag_ctrl, time_type -use MOM_diag_mediator, only : post_data, register_diag_field -use MOM_EOS, only : EOS_type, calculate_compress, calculate_density_derivs -use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg, is_root_pe -use MOM_file_parser, only : get_param, log_version, param_file_type -use MOM_file_parser, only : openParameterBlock, closeParameterBlock -use MOM_grid, only : ocean_grid_type -use MOM_tracer_registry,only : tracer_registry_type -use MOM_verticalGrid, only : verticalGrid_type +use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end +use MOM_cpu_clock, only : CLOCK_MODULE, CLOCK_ROUTINE +use MOM_diag_mediator, only : diag_ctrl, time_type +use MOM_diag_mediator, only : post_data, register_diag_field +use MOM_EOS, only : EOS_type, EOS_manual_init, calculate_compress, calculate_density_derivs +use MOM_EOS, only : calculate_density_second_derivs +use MOM_EOS, only : extract_member_EOS, EOS_LINEAR, EOS_TEOS10, EOS_WRIGHT +use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg, is_root_pe +use MOM_file_parser, only : get_param, log_version, param_file_type +use MOM_file_parser, only : openParameterBlock, closeParameterBlock +use MOM_grid, only : ocean_grid_type +use MOM_remapping, only : remapping_CS, initialize_remapping +use MOM_remapping, only : extract_member_remapping_CS, build_reconstructions_1d +use MOM_remapping, only : average_value_ppoly, remappingSchemesDoc, remappingDefaultScheme +use MOM_tracer_registry, only : tracer_registry_type +use MOM_verticalGrid, only : verticalGrid_type +use polynomial_functions, only : evaluation_polynomial, first_derivative_polynomial +use PPM_functions, only : PPM_reconstruction, PPM_boundary_extrapolation +use regrid_edge_values, only : edge_values_implicit_h4 implicit none ; private @@ -28,8 +36,15 @@ module MOM_neutral_diffusion type, public :: neutral_diffusion_CS ; private integer :: nkp1 ! Number of interfaces for a column = nk + 1 - integer :: nkp1X2 ! Number of intersecting interfaces between columns = 2 * nkp1 - + integer :: nsurf ! Number of neutral surfaces + integer :: ppoly_deg = 2 ! Degree of polynomial used for reconstructions + logical :: continuous_reconstruction = .true. ! True if using continuous PPM reconstruction at interfaces + logical :: boundary_extrap = .true. + logical :: refine_position = .false. + integer :: max_iter ! Maximum number of iterations if refine_position is defined + real :: tolerance ! Convergence criterion representing difference from true neutrality + + ! Positions of neutral surfaces in both the u, v directions real, allocatable, dimension(:,:,:) :: uPoL ! Non-dimensional position with left layer uKoL-1, u-point real, allocatable, dimension(:,:,:) :: uPoR ! Non-dimensional position with right layer uKoR-1, u-point integer, allocatable, dimension(:,:,:) :: uKoL ! Index of left interface corresponding to neutral surface, u-point @@ -40,6 +55,20 @@ module MOM_neutral_diffusion integer, allocatable, dimension(:,:,:) :: vKoL ! Index of left interface corresponding to neutral surface, v-point integer, allocatable, dimension(:,:,:) :: vKoR ! Index of right interface corresponding to neutral surface, v-point real, allocatable, dimension(:,:,:) :: vHeff ! Effective thickness at v-point (H units) + ! Coefficients of polynomial reconstructions for temperature and salinity + real, allocatable, dimension(:,:,:,:) :: ppoly_coeffs_T !< Polynomial coefficients for temperature + real, allocatable, dimension(:,:,:,:) :: ppoly_coeffs_S !< Polynomial coefficients for temperature + ! Variables needed for continuous reconstructions + real, allocatable, dimension(:,:,:) :: dRdT ! dRho/dT (kg/m3/degC) at interfaces + real, allocatable, dimension(:,:,:) :: dRdS ! dRho/dS (kg/m3/ppt) at interfaces + real, allocatable, dimension(:,:,:) :: Tint ! Interface T (degC) + real, allocatable, dimension(:,:,:) :: Sint ! Interface S (ppt) + real, allocatable, dimension(:,:,:) :: Pint ! Interface pressure (Pa) + ! Variables needed for discontinuous reconstructions + real, allocatable, dimension(:,:,:,:) :: T_i ! Top edge reconstruction of temperature (degC) + real, allocatable, dimension(:,:,:,:) :: S_i ! Top edge reconstruction of salinity (ppt) + real, allocatable, dimension(:,:,:,:) :: dRdT_i ! dRho/dT (kg/m3/degC) at top edge + real, allocatable, dimension(:,:,:,:) :: dRdS_i ! dRho/dS (kg/m3/ppt) at top edge type(diag_ctrl), pointer :: diag ! structure to regulate output integer, allocatable, dimension(:) :: id_neutral_diff_tracer_conc_tend ! tracer concentration tendency @@ -50,13 +79,14 @@ module MOM_neutral_diffusion real :: C_p ! heat capacity of seawater (J kg-1 K-1) + type(remapping_CS) :: remap_CS end type neutral_diffusion_CS ! This include declares and sets the variable "version". #include "version_variable.h" character(len=40) :: mdl = "MOM_neutral_diffusion" ! module name -logical, parameter :: debug_this_module = .false. +logical :: debug_this_module = .false. ! If true, verbose output of find neutral position contains @@ -70,6 +100,7 @@ logical function neutral_diffusion_init(Time, G, param_file, diag, CS) ! Local variables character(len=256) :: mesg ! Message for error messages. + character(len=80) :: string ! Temporary strings if (associated(CS)) then call MOM_error(FATAL, "neutral_diffusion_init called with associated control structure.") @@ -89,31 +120,88 @@ logical function neutral_diffusion_init(Time, G, param_file, diag, CS) allocate(CS) CS%diag => diag - + ! call openParameterBlock(param_file,'NEUTRAL_DIFF') ! Read all relevant parameters and write them to the model log. -! call openParameterBlock(param_file,'NEUTRAL_DIFF') + call get_param(param_file, mdl, "NDIFF_CONTINUOUS", CS%continuous_reconstruction, & + "If true, uses a continuous reconstruction of T and S when \n"// & + "finding neutral surfaces along which diffusion will happen.\n"// & + "If false, a PPM discontinuous reconstruction of T and S \n"// & + "is done which results in a higher order routine but exacts \n"// & + "a higher computational cost.", default=.true.) + ! Initialize and configure remapping + if (CS%continuous_reconstruction .eqv. .false.) then + call get_param(param_file, mdl, "NDIFF_BOUNDARY_EXTRAP", CS%boundary_extrap, & + "Uses a rootfinding approach to find the position of a\n"// & + "neutral surface within a layer taking into account the\n"// & + "nonlinearity of the equation of state and the\n"// & + "polynomial reconstructions of T/S.", & + default=.false.) + call get_param(param_file, mdl, "NDIFF_REMAPPING_SCHEME", string, & + "This sets the reconstruction scheme used\n"//& + "for vertical remapping for all variables.\n"//& + "It can be one of the following schemes:\n"//& + trim(remappingSchemesDoc), default=remappingDefaultScheme) + call initialize_remapping( CS%remap_CS, string, boundary_extrapolation = CS%boundary_extrap ) + call extract_member_remapping_CS(CS%remap_CS, degree=CS%ppoly_deg) + call get_param(param_file, mdl, "NDIFF_REFINE_POSITION", CS%refine_position, & + "Uses a rootfinding approach to find the position of a\n"// & + "neutral surface within a layer taking into account the\n"// & + "nonlinearity of the equation of state and the\n"// & + "polynomial reconstructions of T/S.", & + default=.false.) + if (CS%refine_position) then + call get_param(param_file, mdl, "NDIFF_TOLERANCE", CS%tolerance, & + "Sets the convergence criterion for finding the neutral\n"// & + "position within a layer in kg m-3.", & + default=1.e-10) + call get_param(param_file, mdl, "NDIFF_MAX_ITER", CS%max_iter, & + "The maximum number of iterations to be done before \n"// & + "exiting the iterative loop to find the neutral surface", & + default=10) + endif + call get_param(param_file, mdl, "NDIFF_DEBUG", debug_this_module, & + "Turns on verbose output for discontinuous neutral \n"// & + "diffusion routines.", & + default = .false.) + endif + ! call get_param(param_file, mdl, "KHTR", CS%KhTr, & ! "The background along-isopycnal tracer diffusivity.", & ! units="m2 s-1", default=0.0) -! call closeParameterBlock(param_file) - +! call closeParameterBlock(param_file) + if (CS%continuous_reconstruction) then + CS%nsurf = 2*G%ke+2 ! Continuous reconstruction means that every interface has two connections + allocate(CS%dRdT(SZI_(G),SZJ_(G),SZK_(G)+1)) ; CS%dRdT(:,:,:) = 0. + allocate(CS%dRdS(SZI_(G),SZJ_(G),SZK_(G)+1)) ; CS%dRdS(:,:,:) = 0. + else + CS%nsurf = 4*G%ke ! Discontinuous means that every interface has four connections + allocate(CS%T_i(SZI_(G),SZJ_(G),SZK_(G),2)) ; CS%T_i(:,:,:,:) = 0. + allocate(CS%S_i(SZI_(G),SZJ_(G),SZK_(G),2)) ; CS%S_i(:,:,:,:) = 0. + allocate(CS%dRdT_i(SZI_(G),SZJ_(G),SZK_(G),2)) ; CS%dRdT_i(:,:,:,:) = 0. + allocate(CS%dRdS_i(SZI_(G),SZJ_(G),SZK_(G),2)) ; CS%dRdS_i(:,:,:,:) = 0. + allocate(CS%ppoly_coeffs_T(SZI_(G),SZJ_(G),SZK_(G),CS%ppoly_deg+1)) ; CS%ppoly_coeffs_T(:,:,:,:) = 0. + allocate(CS%ppoly_coeffs_S(SZI_(G),SZJ_(G),SZK_(G),CS%ppoly_deg+1)) ; CS%ppoly_coeffs_S(:,:,:,:) = 0. + endif + ! T-points + allocate(CS%Tint(SZI_(G),SZJ_(G),SZK_(G)+1)) ; CS%Tint(:,:,:) = 0. + allocate(CS%Sint(SZI_(G),SZJ_(G),SZK_(G)+1)) ; CS%Sint(:,:,:) = 0. + allocate(CS%Pint(SZI_(G),SZJ_(G),SZK_(G)+1)) ; CS%Pint(:,:,:) = 0. ! U-points - allocate(CS%uPoL(G%isd:G%ied,G%jsd:G%jed,2*G%ke+2)); CS%uPoL(G%isc-1:G%iec,G%jsc:G%jec,:) = 0. - allocate(CS%uPoR(G%isd:G%ied,G%jsd:G%jed,2*G%ke+2)); CS%uPoR(G%isc-1:G%iec,G%jsc:G%jec,:) = 0. - allocate(CS%uKoL(G%isd:G%ied,G%jsd:G%jed,2*G%ke+2)); CS%uKoL(G%isc-1:G%iec,G%jsc:G%jec,:) = 0 - allocate(CS%uKoR(G%isd:G%ied,G%jsd:G%jed,2*G%ke+2)); CS%uKoR(G%isc-1:G%iec,G%jsc:G%jec,:) = 0 - allocate(CS%uHeff(G%isd:G%ied,G%jsd:G%jed,2*G%ke+1)); CS%uHeff(G%isc-1:G%iec,G%jsc:G%jec,:) = 0 + allocate(CS%uPoL(G%isd:G%ied,G%jsd:G%jed, CS%nsurf)); CS%uPoL(G%isc-1:G%iec,G%jsc:G%jec,:) = 0. + allocate(CS%uPoR(G%isd:G%ied,G%jsd:G%jed, CS%nsurf)); CS%uPoR(G%isc-1:G%iec,G%jsc:G%jec,:) = 0. + allocate(CS%uKoL(G%isd:G%ied,G%jsd:G%jed, CS%nsurf)); CS%uKoL(G%isc-1:G%iec,G%jsc:G%jec,:) = 0 + allocate(CS%uKoR(G%isd:G%ied,G%jsd:G%jed, CS%nsurf)); CS%uKoR(G%isc-1:G%iec,G%jsc:G%jec,:) = 0 + allocate(CS%uHeff(G%isd:G%ied,G%jsd:G%jed,CS%nsurf-1)); CS%uHeff(G%isc-1:G%iec,G%jsc:G%jec,:) = 0 ! V-points - allocate(CS%vPoL(G%isd:G%ied,G%jsd:G%jed,2*G%ke+2)); CS%vPoL(G%isc:G%iec,G%jsc-1:G%jec,:) = 0. - allocate(CS%vPoR(G%isd:G%ied,G%jsd:G%jed,2*G%ke+2)); CS%vPoR(G%isc:G%iec,G%jsc-1:G%jec,:) = 0. - allocate(CS%vKoL(G%isd:G%ied,G%jsd:G%jed,2*G%ke+2)); CS%vKoL(G%isc:G%iec,G%jsc-1:G%jec,:) = 0 - allocate(CS%vKoR(G%isd:G%ied,G%jsd:G%jed,2*G%ke+2)); CS%vKoR(G%isc:G%iec,G%jsc-1:G%jec,:) = 0 - allocate(CS%vHeff(G%isd:G%ied,G%jsd:G%jed,2*G%ke+1)); CS%vHeff(G%isc:G%iec,G%jsc-1:G%jec,:) = 0 + allocate(CS%vPoL(G%isd:G%ied,G%jsd:G%jed, CS%nsurf)); CS%vPoL(G%isc:G%iec,G%jsc-1:G%jec,:) = 0. + allocate(CS%vPoR(G%isd:G%ied,G%jsd:G%jed, CS%nsurf)); CS%vPoR(G%isc:G%iec,G%jsc-1:G%jec,:) = 0. + allocate(CS%vKoL(G%isd:G%ied,G%jsd:G%jed, CS%nsurf)); CS%vKoL(G%isc:G%iec,G%jsc-1:G%jec,:) = 0 + allocate(CS%vKoR(G%isd:G%ied,G%jsd:G%jed, CS%nsurf)); CS%vKoR(G%isc:G%iec,G%jsc-1:G%jec,:) = 0 + allocate(CS%vHeff(G%isd:G%ied,G%jsd:G%jed,CS%nsurf-1)); CS%vHeff(G%isc:G%iec,G%jsc-1:G%jec,:) = 0 end function neutral_diffusion_init - !> Diagnostic handles for neutral diffusion tendencies. subroutine neutral_diffusion_diag_init(Time, G, diag, C_p, Reg, CS) type(time_type),target, intent(in) :: Time !< Time structure @@ -263,54 +351,112 @@ subroutine neutral_diffusion_calc_coeffs(G, GV, h, T, S, EOS, CS) ! Local variables integer :: i, j, k - real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: Tint ! Interface T (degC) - real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: Sint ! Interface S (ppt) - real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: Pint ! Interface pressure (Pa) - real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: dRdT ! Interface thermal expansion coefficient (kg/m3/degC) - real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: dRdS ! Interface haline expansion coefficient (kg/m3/ppt) + ! Variables used for reconstructions + real, dimension(SZK_(G),2) :: ppoly_r_S ! Reconstruction slopes + integer :: iMethod + + if (CS%continuous_reconstruction) then + CS%dRdT(:,:,:) = 0. + CS%dRdS(:,:,:) = 0. + else + CS%T_i(:,:,:,:) = 0. + CS%S_i(:,:,:,:) = 0. + CS%dRdT_i(:,:,:,:) = 0. + CS%dRdS_i(:,:,:,:) = 0. + endif do j = G%jsc-1, G%jec+1 ! Interpolate state to interface do i = G%isc-1, G%iec+1 - call interface_scalar(G%ke, h(i,j,:), T(i,j,:), Tint(i,j,:), 2) - call interface_scalar(G%ke, h(i,j,:), S(i,j,:), Sint(i,j,:), 2) + if (CS%continuous_reconstruction) then + call interface_scalar(G%ke, h(i,j,:), T(i,j,:), CS%Tint(i,j,:), 2) + call interface_scalar(G%ke, h(i,j,:), S(i,j,:), CS%Sint(i,j,:), 2) + else + call build_reconstructions_1d( CS%remap_CS, G%ke, h(i,j,:), T(i,j,:), CS%ppoly_coeffs_T(i,j,:,:), & + CS%T_i(i,j,:,:), ppoly_r_S, iMethod ) + call build_reconstructions_1d( CS%remap_CS, G%ke, h(i,j,:), S(i,j,:), CS%ppoly_coeffs_S(i,j,:,:), & + CS%S_i(i,j,:,:), ppoly_r_S, iMethod ) + endif enddo ! Calculate interface properties - Pint(:,j,1) = 0. ! Assume P=0 (Pa) at surface - needs correcting for atmospheric and ice loading - AJA - do k = 1, G%ke+1 - call calculate_density_derivs(Tint(:,j,k), Sint(:,j,k), Pint(:,j,k), & - dRdT(:,j,k), dRdS(:,j,k), G%isc-1, G%iec-G%isc+3, EOS) - if (k<=G%ke) Pint(:,j,k+1) = Pint(:,j,k) + h(:,j,k) * GV%H_to_Pa ! Pressure at next interface, k+1 (Pa) - enddo + CS%Pint(:,j,1) = 0. ! Assume P=0 (Pa) at surface - needs correcting for atmospheric and ice loading - AJA + ! Continuous reconstruction + if (CS%continuous_reconstruction) then + do k = 1, G%ke+1 + call calculate_density_derivs(CS%Tint(:,j,k), CS%Sint(:,j,k), CS%Pint(:,j,k), & + CS%dRdT(:,j,k), CS%dRdS(:,j,k), G%isc-1, G%iec-G%isc+3, EOS) + if (k<=G%ke) then + CS%Pint(:,j,k+1) = CS%Pint(:,j,k) + h(:,j,k) * GV%H_to_Pa ! Pressure at next interface, k+1 (Pa) + endif + enddo + else ! Discontinuous reconstruction + do k = 1, G%ke + call calculate_density_derivs(CS%T_i(:,j,k,1), CS%S_i(:,j,k,1), CS%Pint(:,j,k), & + CS%dRdT_i(:,j,k,1), CS%dRdS_i(:,j,k,1), G%isc-1, G%iec-G%isc+3, EOS) + if (k<=G%ke) then + CS%Pint(:,j,k+1) = CS%Pint(:,j,k) + h(:,j,k) * GV%H_to_Pa ! Pressure at next interface, k+1 (Pa) + call calculate_density_derivs(CS%T_i(:,j,k,2), CS%S_i(:,j,k,2), CS%Pint(:,j,k+1), & + CS%dRdT_i(:,j,k,2), CS%dRdS_i(:,j,k,2), G%isc-1, G%iec-G%isc+3, EOS) + endif + enddo + endif enddo + CS%uhEff(:,:,:) = 0. + CS%vhEff(:,:,:) = 0. + CS%uPoL(:,:,:) = 0. + CS%vPoL(:,:,:) = 0. + CS%uPoR(:,:,:) = 0. + CS%vPoR(:,:,:) = 0. + CS%uKoL(:,:,:) = 1 + CS%vKoL(:,:,:) = 1 + CS%uKoR(:,:,:) = 1 + CS%vKoR(:,:,:) = 1 + ! Neutral surface factors at U points - do j = G%jsc, G%jec - do I = G%isc-1, G%iec - call find_neutral_surface_positions(G%ke, & - Pint(i,j,:), Tint(i,j,:), Sint(i,j,:), dRdT(i,j,:), dRdS(i,j,:), & - Pint(i+1,j,:), Tint(i+1,j,:), Sint(i+1,j,:), dRdT(i+1,j,:), dRdS(i+1,j,:),& - CS%uPoL(I,j,:), CS%uPoR(I,j,:), CS%uKoL(I,j,:), CS%uKoR(I,j,:), CS%uhEff(I,j,:) ) - enddo - enddo + do j = G%jsc, G%jec ; do I = G%isc-1, G%iec + if (G%mask2dCu(I,j) > 0.) then + if (CS%continuous_reconstruction) then + call find_neutral_surface_positions_continuous(G%ke, & + CS%Pint(i,j,:), CS%Tint(i,j,:), CS%Sint(i,j,:), CS%dRdT(i,j,:), CS%dRdS(i,j,:), & + CS%Pint(i+1,j,:), CS%Tint(i+1,j,:), CS%Sint(i+1,j,:), CS%dRdT(i+1,j,:), CS%dRdS(i+1,j,:), & + CS%uPoL(I,j,:), CS%uPoR(I,j,:), CS%uKoL(I,j,:), CS%uKoR(I,j,:), CS%uhEff(I,j,:) ) + else + call find_neutral_surface_positions_discontinuous(G%ke, CS%ppoly_deg, & + CS%Pint(i,j,:), CS%T_i(i,j,:,:), CS%S_i(i,j,:,:), CS%dRdT_i(i,j,:,:), CS%dRdS_i(i,j,:,:), & + CS%Pint(i+1,j,:), CS%T_i(i+1,j,:,:), CS%S_i(i+1,j,:,:), CS%dRdT_i(i+1,j,:,:), CS%dRdS_i(i+1,j,:,:), & + CS%uPoL(I,j,:), CS%uPoR(I,j,:), CS%uKoL(I,j,:), CS%uKoR(I,j,:), CS%uhEff(I,j,:), & + CS%refine_position, CS%ppoly_coeffs_T(i,j,:,:), CS%ppoly_coeffs_S(i,j,:,:), & + CS%ppoly_coeffs_T(i+1,j,:,:), CS%ppoly_coeffs_S(i+1,j,:,:), EOS, CS%max_iter, CS%tolerance) + endif + endif + enddo ; enddo ! Neutral surface factors at V points - do J = G%jsc-1, G%jec - do i = G%isc, G%iec - call find_neutral_surface_positions(G%ke, & - Pint(i,j,:), Tint(i,j,:), Sint(i,j,:), dRdT(i,j,:), dRdS(i,j,:), & - Pint(i,j+1,:), Tint(i,j+1,:), Sint(i,j+1,:), dRdT(i,j+1,:), dRdS(i,j+1,:),& - CS%vPoL(i,J,:), CS%vPoR(i,J,:), CS%vKoL(i,J,:), CS%vKoR(i,J,:), CS%vhEff(i,J,:) ) - enddo - enddo + do J = G%jsc-1, G%jec ; do i = G%isc, G%iec + if (G%mask2dCv(i,J) > 0.) then + if (CS%continuous_reconstruction) then + call find_neutral_surface_positions_continuous(G%ke, & + CS%Pint(i,j,:), CS%Tint(i,j,:), CS%Sint(i,j,:), CS%dRdT(i,j,:), CS%dRdS(i,j,:), & + CS%Pint(i,j+1,:), CS%Tint(i,j+1,:), CS%Sint(i,j+1,:), CS%dRdT(i,j+1,:), CS%dRdS(i,j+1,:), & + CS%vPoL(i,J,:), CS%vPoR(i,J,:), CS%vKoL(i,J,:), CS%vKoR(i,J,:), CS%vhEff(i,J,:) ) + else + call find_neutral_surface_positions_discontinuous(G%ke, CS%ppoly_deg, & + CS%Pint(i,j,:), CS%T_i(i,j,:,:), CS%S_i(i,j,:,:), CS%dRdT_i(i,j,:,:), CS%dRdS_i(i,j,:,:), & + CS%Pint(i,j+1,:), CS%T_i(i,j+1,:,:), CS%S_i(i,j+1,:,:), CS%dRdT_i(i,j+1,:,:), CS%dRdS_i(i,j+1,:,:), & + CS%vPoL(I,j,:), CS%vPoR(I,j,:), CS%vKoL(I,j,:), CS%vKoR(I,j,:), CS%vhEff(I,j,:), & + CS%refine_position, CS%ppoly_coeffs_T(i,j,:,:), CS%ppoly_coeffs_S(i,j,:,:), & + CS%ppoly_coeffs_T(i,j+1,:,:), CS%ppoly_coeffs_S(i,j+1,:,:), EOS, CS%max_iter, CS%tolerance) + endif + endif + enddo ; enddo CS%uhEff(:,:,:) = CS%uhEff(:,:,:) / GV%H_to_pa CS%vhEff(:,:,:) = CS%vhEff(:,:,:) / GV%H_to_pa end subroutine neutral_diffusion_calc_coeffs - !> Update tracer concentration due to neutral diffusion; layer thickness unchanged by this update. subroutine neutral_diffusion(G, GV, h, Coef_x, Coef_y, Tracer, m, dt, name, CS) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure @@ -325,13 +471,13 @@ subroutine neutral_diffusion(G, GV, h, Coef_x, Coef_y, Tracer, m, dt, name, CS) type(neutral_diffusion_CS), pointer :: CS !< Neutral diffusion control structure ! Local variables - real, dimension(SZIB_(G),SZJ_(G),2*G%ke+1) :: uFlx ! Zonal flux of tracer (concentration * H) - real, dimension(SZI_(G),SZJB_(G),2*G%ke+1) :: vFlx ! Meridional flux of tracer (concentration * H) - real, dimension(SZI_(G),SZJ_(G),G%ke) :: tendency ! tendency array for diagn - real, dimension(SZI_(G),SZJ_(G)) :: tendency_2d ! depth integrated content tendency for diagn - real, dimension(SZIB_(G),SZJ_(G)) :: trans_x_2d ! depth integrated diffusive tracer x-transport diagn - real, dimension(SZI_(G),SZJB_(G)) :: trans_y_2d ! depth integrated diffusive tracer y-transport diagn - real, dimension(G%ke) :: dTracer ! change in tracer concentration due to ndiffusion + real, dimension(SZIB_(G),SZJ_(G),CS%nsurf-1) :: uFlx ! Zonal flux of tracer (concentration * H) + real, dimension(SZI_(G),SZJB_(G),CS%nsurf-1) :: vFlx ! Meridional flux of tracer (concentration * H) + real, dimension(SZI_(G),SZJ_(G),G%ke) :: tendency ! tendency array for diagn + real, dimension(SZI_(G),SZJ_(G)) :: tendency_2d ! depth integrated content tendency for diagn + real, dimension(SZIB_(G),SZJ_(G)) :: trans_x_2d ! depth integrated diffusive tracer x-transport diagn + real, dimension(SZI_(G),SZJB_(G)) :: trans_y_2d ! depth integrated diffusive tracer y-transport diagn + real, dimension(G%ke) :: dTracer ! change in tracer concentration due to ndiffusion integer :: i, j, k, ks, nk real :: ppt2mks, Idt, convert @@ -354,30 +500,30 @@ subroutine neutral_diffusion(G, GV, h, Coef_x, Coef_y, Tracer, m, dt, name, CS) if(trim(name) == 'S') convert = ppt2mks * GV%H_to_kg_m2 endif + uFlx(:,:,:) = 0. + vFlx(:,:,:) = 0. ! x-flux do j = G%jsc,G%jec ; do I = G%isc-1,G%iec if (G%mask2dCu(I,j)>0.) then - call neutral_surface_flux(nk, h(i,j,:), h(i+1,j,:), & + call neutral_surface_flux(nk, CS%nsurf, CS%ppoly_deg, h(i,j,:), h(i+1,j,:), & Tracer(i,j,:), Tracer(i+1,j,:), & CS%uPoL(I,j,:), CS%uPoR(I,j,:), & CS%uKoL(I,j,:), CS%uKoR(I,j,:), & - CS%uhEff(I,j,:), uFlx(I,j,:)) - else - uFlx(I,j,:) = 0. + CS%uhEff(I,j,:), uFlx(I,j,:), & + CS%continuous_reconstruction, CS%remap_CS) endif enddo ; enddo ! y-flux do J = G%jsc-1,G%jec ; do i = G%isc,G%iec if (G%mask2dCv(i,J)>0.) then - call neutral_surface_flux(nk, h(i,j,:), h(i,j+1,:), & + call neutral_surface_flux(nk, CS%nsurf, CS%ppoly_deg, h(i,j,:), h(i,j+1,:), & Tracer(i,j,:), Tracer(i,j+1,:), & CS%vPoL(i,J,:), CS%vPoR(i,J,:), & CS%vKoL(i,J,:), CS%vKoR(i,J,:), & - CS%vhEff(i,J,:), vFlx(i,J,:)) - else - vFlx(I,j,:) = 0. + CS%vhEff(i,J,:), vFlx(i,J,:), & + CS%continuous_reconstruction, CS%remap_CS) endif enddo ; enddo @@ -386,7 +532,7 @@ subroutine neutral_diffusion(G, GV, h, Coef_x, Coef_y, Tracer, m, dt, name, CS) if (G%mask2dT(i,j)>0.) then dTracer(:) = 0. - do ks = 1,2*nk+1 ; + do ks = 1,CS%nsurf-1 ; k = CS%uKoL(I,j,ks) dTracer(k) = dTracer(k) + Coef_x(I,j) * uFlx(I,j,ks) k = CS%uKoR(I-1,j,ks) @@ -419,7 +565,7 @@ subroutine neutral_diffusion(G, GV, h, Coef_x, Coef_y, Tracer, m, dt, name, CS) do j = G%jsc,G%jec ; do I = G%isc-1,G%iec trans_x_2d(I,j) = 0. if (G%mask2dCu(I,j)>0.) then - do ks = 1,2*nk+1 ; + do ks = 1,CS%nsurf-1 ; trans_x_2d(I,j) = trans_x_2d(I,j) - Coef_x(I,j) * uFlx(I,j,ks) enddo trans_x_2d(I,j) = trans_x_2d(I,j) * Idt * convert @@ -434,7 +580,7 @@ subroutine neutral_diffusion(G, GV, h, Coef_x, Coef_y, Tracer, m, dt, name, CS) do J = G%jsc-1,G%jec ; do i = G%isc,G%iec trans_y_2d(i,J) = 0. if (G%mask2dCv(i,J)>0.) then - do ks = 1,2*nk+1 ; + do ks = 1,CS%nsurf-1 ; trans_y_2d(i,J) = trans_y_2d(i,J) - Coef_y(i,J) * vFlx(i,J,ks) enddo trans_y_2d(i,J) = trans_y_2d(i,J) * Idt * convert @@ -471,7 +617,6 @@ subroutine neutral_diffusion(G, GV, h, Coef_x, Coef_y, Tracer, m, dt, name, CS) end subroutine neutral_diffusion - !> Returns interface scalar, Si, for a column of layer values, S. subroutine interface_scalar(nk, h, S, Si, i_method) integer, intent(in) :: nk !< Number of levels @@ -574,7 +719,6 @@ real function ppm_ave(xL, xR, aL, aR, aMean) else ppm_ave = ( aL + xave * ( ( aR - aL ) + a6 ) ) - a6o3 * ( xR**2 + xR * xL + xL**2 ) endif - end function ppm_ave !> A true signum function that returns either -abs(a), when x<0; or abs(a) when x>0; or 0 when x=0. @@ -711,8 +855,9 @@ real function fvlsq_slope(hkm1, hk, hkp1, Skm1, Sk, Skp1) end function fvlsq_slope -!> Returns positions within left/right columns of combined interfaces -subroutine find_neutral_surface_positions(nk, Pl, Tl, Sl, dRdTl, dRdSl, Pr, Tr, Sr, dRdTr, dRdSr, PoL, PoR, KoL, KoR, hEff) +!> Returns positions within left/right columns of combined interfaces using continuous reconstructions of T/S +subroutine find_neutral_surface_positions_continuous(nk, Pl, Tl, Sl, dRdTl, dRdSl, Pr, Tr, Sr, dRdTr, dRdSr, PoL, & + PoR, KoL, KoR, hEff) integer, intent(in) :: nk !< Number of levels real, dimension(nk+1), intent(in) :: Pl !< Left-column interface pressure (Pa) real, dimension(nk+1), intent(in) :: Tl !< Left-column interface potential temperature (degC) @@ -731,6 +876,7 @@ subroutine find_neutral_surface_positions(nk, Pl, Tl, Sl, dRdTl, dRdSl, Pr, Tr, real, dimension(2*nk+1), intent(inout) :: hEff !< Effective thickness between two neutral surfaces (Pa) ! Local variables + integer :: ns ! Number of neutral surfaces integer :: k_surface ! Index of neutral surface integer :: kl ! Index of left interface integer :: kr ! Index of right interface @@ -743,13 +889,14 @@ subroutine find_neutral_surface_positions(nk, Pl, Tl, Sl, dRdTl, dRdSl, Pr, Tr, integer :: lastK_left, lastK_right real :: lastP_left, lastP_right + ns = 2*nk+2 ! Initialize variables for the search kr = 1 ; lastK_right = 1 ; lastP_right = 0. kl = 1 ; lastK_left = 1 ; lastP_left = 0. reached_bottom = .false. ! Loop over each neutral surface, working from top to bottom - neutral_surfaces: do k_surface = 1, 2*nk+2 + neutral_surfaces: do k_surface = 1, ns klm1 = max(kl-1, 1) if (klm1>nk) stop 'find_neutral_surface_positions(): klm1 went out of bounds!' krm1 = max(kr-1, 1) @@ -784,7 +931,7 @@ subroutine find_neutral_surface_positions(nk, Pl, Tl, Sl, dRdTl, dRdSl, Pr, Tr, + ( dRdSl(klm1) + dRdSr(kr) ) * ( Sl(klm1) - Sr(kr) ) ) ! Potential density difference, rho(kl) - rho(kr) (will be positive) dRhoBot = 0.5 * ( ( dRdTl(klm1+1) + dRdTr(kr) ) * ( Tl(klm1+1) - Tr(kr) ) & - + ( dRdSl(klm1+1) + dRdSr(kr) ) * ( Sl(klm1+1) - Sr(kr) ) ) + + ( dRdSl(klm1+1) + dRdSr(kr) ) * ( Sl(klm1+1) - Sr(kr) ) ) ! Because we are looking left, the right surface, kr, is lighter than klm1+1 and should be denser than klm1 ! unless we are still at the top of the left column (kl=1) @@ -874,10 +1021,10 @@ subroutine find_neutral_surface_positions(nk, Pl, Tl, Sl, dRdTl, dRdSl, Pr, Tr, ! NOTE: This would be better expressed in terms of the layers thicknesses rather ! than as differences of position - AJA if (k_surface>1) then - hL = absolute_position(nk,Pl,KoL,PoL,k_surface) - absolute_position(nk,Pl,KoL,PoL,k_surface-1) - hR = absolute_position(nk,Pr,KoR,PoR,k_surface) - absolute_position(nk,Pr,KoR,PoR,k_surface-1) + hL = absolute_position(nk,ns,Pl,KoL,PoL,k_surface) - absolute_position(nk,ns,Pl,KoL,PoL,k_surface-1) + hR = absolute_position(nk,ns,Pr,KoR,PoR,k_surface) - absolute_position(nk,ns,Pr,KoR,PoR,k_surface-1) if ( hL + hR > 0.) then - hEff(k_surface-1) = 2. * hL * hR / ( hL + hR ) + hEff(k_surface-1) = 2. * hL * hR / ( hL + hR ) ! Harmonic mean of layer thicknesses else hEff(k_surface-1) = 0. endif @@ -885,14 +1032,333 @@ subroutine find_neutral_surface_positions(nk, Pl, Tl, Sl, dRdTl, dRdSl, Pr, Tr, enddo neutral_surfaces -end subroutine find_neutral_surface_positions +end subroutine find_neutral_surface_positions_continuous + +!> Higher order version of find_neutral_surface_positions. Returns positions within left/right columns +!! of combined interfaces using intracell reconstructions of T/S +subroutine find_neutral_surface_positions_discontinuous(nk, deg, & + Pres_l, Tl, Sl, dRdT_l, dRdS_l, Pres_r, Tr, Sr, dRdT_r, dRdS_r, PoL, PoR, KoL, KoR, hEff, & + refine_pos_in, ppoly_T_l, ppoly_S_l, ppoly_T_r, ppoly_S_r, EOS, max_iter, tolerance) + integer, intent(in) :: nk !< Number of levels + integer, intent(in) :: deg !< Degree of polynomial used for reconstructions + real, dimension(nk+1), intent(in) :: Pres_l !< Left-column interface pressure (Pa) + real, dimension(nk,2), intent(in) :: Tl !< Left-column top interface potential temperature (degC) + real, dimension(nk,2), intent(in) :: Sl !< Left-column top interface salinity (ppt) + real, dimension(nk,2), intent(in) :: dRdT_l !< Left-column, top interface dRho/dT (kg/m3/degC) + real, dimension(nk,2), intent(in) :: dRdS_l !< Left-column, top interface dRho/dS (kg/m3/ppt) + real, dimension(nk+1), intent(in) :: Pres_r !< Right-column interface pressure (Pa) + real, dimension(nk,2), intent(in) :: Tr !< Right-column top interface potential temperature (degC) + real, dimension(nk,2), intent(in) :: Sr !< Right-column top interface salinity (ppt) + real, dimension(nk,2), intent(in) :: dRdT_r !< Right-column, top interface dRho/dT (kg/m3/degC) + real, dimension(nk,2), intent(in) :: dRdS_r !< Right-column, top interface dRho/dS (kg/m3/ppt) + real, dimension(4*nk), intent(inout) :: PoL !< Fractional position of neutral surface within + !! layer KoL of left column + real, dimension(4*nk), intent(inout) :: PoR !< Fractional position of neutral surface within + !! layer KoR of right column + integer, dimension(4*nk), intent(inout) :: KoL !< Index of first left interface above neutral surface + integer, dimension(4*nk), intent(inout) :: KoR !< Index of first right interface above neutral surface + real, dimension(4*nk-1), intent(inout) :: hEff !< Effective thickness between two neutral surfaces (Pa) + logical, optional, intent(in) :: refine_pos_in !< True if rootfinding is used for position + real, dimension(nk,deg+1), optional, intent(in) :: ppoly_T_l !< Left-column coefficients of T reconstruction + real, dimension(nk,deg+1), optional, intent(in) :: ppoly_S_l !< Left-column coefficients of S reconstruction + real, dimension(nk,deg+1), optional, intent(in) :: ppoly_T_r !< Right-column coefficients of T reconstruction + real, dimension(nk,deg+1), optional, intent(in) :: ppoly_S_r !< Right-column coefficients of S reconstruction + type(EOS_type), optional, pointer :: EOS !< Equation of state structure + integer, optional, intent(in) :: max_iter !< Maximum number of iterations in refine_position + real, optional, intent(in) :: tolerance !< Convergence criterion for refine_position + + ! Local variables + integer :: ns ! Number of neutral surfaces + integer :: k_surface ! Index of neutral surface + integer :: kl_left, kl_right ! Index of layers on the left/right + integer :: ki_left, ki_right ! Index of interfaces on the left/right + logical :: searching_left_column ! True if searching for the position of a right interface in the left column + logical :: searching_right_column ! True if searching for the position of a left interface in the right column + logical :: reached_bottom ! True if one of the bottom-most interfaces has been used as the target + logical :: refine_pos ! Use rootfinding to find the true neutral surface position + real :: dRho, dRhoTop, dRhoBot, dRhoTopm1, hL, hR + integer :: lastK_left, lastK_right, maxK_left, maxK_right + real :: lastP_left, lastP_right + real :: min_bound + logical, dimension(nk) :: top_connected_l, top_connected_r + logical, dimension(nk) :: bot_connected_l, bot_connected_r + + ns = 4*nk + top_connected_l(:) = .false. ; top_connected_r(:) = .false. + bot_connected_l(:) = .false. ; bot_connected_r(:) = .false. + maxK_left = -1 ; maxK_right = -1 + ! Vectors with all the values of the discontinuous reconstruction. + ! Dimensions are [number of layers x number of interfaces]. Second dimension = 1 for top interface, = 2 for bottom +! real, dimension(nk,2) :: Sl, Sr, Tl, Tr, dRdT_l, dRdS_l, dRdT_r, dRdS_r + +! Check to make sure that polynomial reconstructions were passed if refine_pos defined) + refine_pos = .false. + if (present(refine_pos_in)) then + refine_pos = refine_pos_in + if (refine_pos .and. (.not. ( present(ppoly_T_l) .and. present(ppoly_S_l) .and. & + present(ppoly_T_r) .and. present(ppoly_S_r) .and. & + present(tolerance) .and. present(max_iter)) ) ) & + call MOM_error(FATAL, "fine_neutral_surface_positions_discontinuous: refine_pos is requested, but polynomial"// & + "coefficients not available for T and S") + endif + + ! Initialize variables for the search + kl_right = 1 ; ki_right = 1 ; lastK_right = 1 ; lastP_right = -1. + kl_left = 1 ; ki_left = 1 ; lastK_left = 1 ; lastP_left = -1. + + reached_bottom = .false. + searching_left_column = .false. + searching_right_column = .false. + + ! Loop over each neutral surface, working from top to bottom + neutral_surfaces: do k_surface = 1, 4*nk + ! Potential density difference, rho(kr) - rho(kl) + dRho = 0.5 * & + ( ( dRdT_r(kl_right,ki_right) + dRdT_l(kl_left,ki_left) ) * ( Tr(kl_right,ki_right) - Tl(kl_left,ki_left) ) & + + ( dRdS_r(kl_right,ki_right) + dRdS_l(kl_left,ki_left) ) * ( Sr(kl_right,ki_right) - Sl(kl_left,ki_left) ) ) + if (debug_this_module) write(*,'(A,I2,A,E12.4,A,I2,A,I2,A,I2,A,I2)') "k_surface=",k_surface," dRho=",dRho," kl_left=",kl_left, & + " ki_left=",ki_left," kl_right=",kl_right, " ki_right=",ki_right + ! Which column has the lighter surface for the current indexes, kr and kl + if (.not. reached_bottom) then + if (dRho < 0.) then + searching_left_column = .true. + searching_right_column = .false. + elseif (dRho > 0.) then + searching_right_column = .true. + searching_left_column = .false. + else ! dRho == 0. + if ((kl_left + kl_left == 2) .and. (ki_left + ki_right == 2)) then ! Still at surface + searching_left_column = .true. + searching_right_column = .false. + else ! Not the surface so we simply change direction + searching_left_column = .not. searching_left_column + searching_right_column = .not. searching_right_column + endif + endif + endif + + if (searching_left_column) then + ! Determine differences between right column interface and potentially three different parts of the left + ! Potential density difference, rho(kl-1) - rho(kr) (should be negative) + dRhoTop = 0.5 * & + ( ( dRdT_l(kl_left,1) + dRdT_r(kl_right,ki_right) ) * ( Tl(kl_left,1) - Tr(kl_right,ki_right) ) & + + ( dRdS_l(kl_left,1) + dRdS_r(kl_right,ki_right) ) * ( Sl(kl_left,1) - Sr(kl_right,ki_right) ) ) + ! Potential density difference, rho(kl) - rho(kl_right,ki_right) (will be positive) + dRhoBot = 0.5 * & + ( ( dRdT_l(kl_left,2) + dRdT_r(kl_right,ki_right) ) * ( Tl(kl_left,2) - Tr(kl_right,ki_right) ) & + + ( dRdS_l(kl_left,2) + dRdS_r(kl_right,ki_right) ) * ( Sl(kl_left,2) - Sr(kl_right,ki_right) ) ) + if (kl_left>1) then ! Calculate the density difference at top of discontinuity + dRhoTopm1 = 0.5 * & + ( ( dRdT_l(kl_left-1,2) + dRdT_r(kl_right,ki_right) ) * ( Tl(kl_left-1,2) - Tr(kl_right,ki_right) ) & + + ( dRdS_l(kl_left-1,2) + dRdS_r(kl_right,ki_right) ) * ( Sl(kl_left-1,2) - Sr(kl_right,ki_right) ) ) + else + dRhoTopm1 = dRhoTop + endif + if (debug_this_module) then + write(*,'(A,I2,A,E12.4,A,E12.4,A,E12.4)') "Searching left layer ", kl_left, ": dRhoTopm1=", dRhoTopm1, & + " dRhoTop=", dRhoTop, " dRhoBot=", dRhoBot + write(*,'(A,I2,X,I2)') "Searching from right: ", kl_right, ki_right + write(*,*) "Temp/Salt Reference: ", Tr(kl_right,ki_right), Sr(kl_right,ki_right) + write(*,*) "Temp/Salt Top L: ", Tl(kl_left,1), Sl(kl_left,1) + write(*,*) "Temp/Salt Bot L: ", Tl(kl_left,2), Sl(kl_left,2) + endif + + ! Set the position within the starting column + PoR(k_surface) = REAL(ki_right-1) + KoR(k_surface) = REAL(kl_right) + + ! Set position within the searched column + call search_other_column_discontinuous(dRhoTopm1, dRhoTop, dRhoBot, Pres_l(kl_left), Pres_l(kl_left+1), & + lastP_left, lastK_left, kl_left, ki_left, top_connected_l, bot_connected_l, PoL(k_surface), KoL(k_surface)) + if ( refine_pos .and. (PoL(k_surface) > 0.) .and. (PoL(k_surface) < 1.) ) then + min_bound = 0. + if ( (k_surface > 1) .and. ( KoL(k_surface) == KoL(k_surface-1) ) ) min_bound = PoL(k_surface-1) + PoL(k_surface) = refine_nondim_position(max_iter, tolerance, Tr(kl_right,ki_right), Sr(kl_right,ki_right), & + dRdT_r(kl_right,ki_right), dRdS_r(kl_right,ki_right), Pres_l(kl_left), Pres_l(kl_left+1), & + deg, ppoly_T_l(kl_left,:), ppoly_S_l(kl_left,:), EOS, PoL(k_surface), dRhoTop, dRhoBot, min_bound) + endif + if (PoL(k_surface) == 0.) top_connected_l(KoL(k_surface)) = .true. + if (PoL(k_surface) == 1.) bot_connected_l(KoL(k_surface)) = .true. + call increment_interface(nk, kl_right, ki_right, reached_bottom, searching_right_column, searching_left_column) + lastK_left = KoL(k_surface) ; lastP_left = PoL(k_surface) + + elseif (searching_right_column) then + ! Interpolate for the neutral surface position within the right column, layer krm1 + ! Potential density difference, rho(kr-1) - rho(kl) (should be negative) + dRhoTop = 0.5 * & + ( ( dRdT_r(kl_right,1) + dRdT_l(kl_left,ki_left) ) * ( Tr(kl_right,1) - Tl(kl_left,ki_left) ) & + + ( dRdS_r(kl_right,1) + dRdS_l(kl_left,ki_left) ) * ( Sr(kl_right,1) - Sl(kl_left,ki_left) ) ) + dRhoBot = 0.5 * & + ( ( dRdT_r(kl_right,2) + dRdT_l(kl_left,ki_left) ) * ( Tr(kl_right,2) - Tl(kl_left,ki_left) ) & + + ( dRdS_r(kl_right,2) + dRdS_l(kl_left,ki_left) ) * ( Sr(kl_right,2) - Sl(kl_left,ki_left) ) ) + if (kl_right>1) then + dRhoTopm1 = 0.5 * & + ( ( dRdT_r(kl_right-1,2) + dRdT_l(kl_left,ki_left) ) * ( Tr(kl_right-1,2) - Tl(kl_left,ki_left) ) & + + ( dRdS_r(kl_right-1,2) + dRdS_l(kl_left,ki_left) ) * ( Sr(kl_right-1,2) - Sl(kl_left,ki_left) ) ) + else + dRhoTopm1 = dRhoTop + endif + if (debug_this_module) then + write(*,'(A,I2,A,E12.4,A,E12.4,A,E12.4)') "Searching right layer ", kl_right, ": dRhoTopm1=", dRhoTopm1, & + " dRhoTop=", dRhoTop, " dRhoBot=", dRhoBot + write(*,'(A,I2,X,I2)') "Searching from left: ", kl_left, ki_left + write(*,*) "Temp/Salt Reference: ", Tl(kl_left,ki_left), Sl(kl_left,ki_left) + write(*,*) "Temp/Salt Top R: ", Tr(kl_right,1), Sr(kl_right,1) + write(*,*) "Temp/Salt Bot R: ", Tr(kl_right,2), Sr(kl_right,2) + endif + ! Set the position within the starting column + PoL(k_surface) = REAL(ki_left-1) + KoL(k_surface) = REAL(kl_left) + + ! Set position within the searched column + call search_other_column_discontinuous(dRhoTopm1, dRhoTop, dRhoBot, Pres_r(kl_right), Pres_r(kl_right+1), & + lastP_right, lastK_right, kl_right, ki_right, top_connected_r, bot_connected_r, PoR(k_surface), KoR(k_surface)) + if ( refine_pos .and. (PoR(k_surface) > 0. .and. PoR(k_surface) < 1.) ) then + min_bound = 0. + if ( (k_surface > 1) .and. ( KoR(k_surface) == KoR(k_surface-1) ) ) min_bound = PoR(k_surface-1) + PoR(k_surface) = refine_nondim_position(max_iter, tolerance, Tl(kl_left,ki_left), Sl(kl_left,ki_left), & + dRdT_l(kl_left,ki_left), dRdS_l(kl_left,ki_left), Pres_r(kl_right), Pres_r(kl_right+1), & + deg, ppoly_T_r(kl_right,:), ppoly_S_r(kl_right,:), EOS, PoR(k_surface), dRhoTop, dRhoBot, min_bound) + endif + if (PoR(k_surface) == 0.) top_connected_r(KoR(k_surface)) = .true. + if (PoR(k_surface) == 1.) bot_connected_r(KoR(k_surface)) = .true. + call increment_interface(nk, kl_left, ki_left, reached_bottom, searching_left_column, searching_right_column) + lastK_right = KoR(k_surface) ; lastP_right = PoR(k_surface) + + else + stop 'Else what?' + endif + if (debug_this_module) write(*,'(A,I2,A,F6.2,A,I2,A,F6.2)') "KoL:", KoL(k_surface), " PoL:", PoL(k_surface), " KoR:", & + KoR(k_surface), " PoR:", PoR(k_surface) + maxK_left= MAX(KoL(k_surface), maxK_left) + maxK_right= MAX(KoR(k_surface), maxK_right) + ! Effective thickness + ! NOTE: This would be better expressed in terms of the layers thicknesses rather + ! than as differences of position - AJA + if (k_surface>1) then + hL = absolute_position(nk,ns,Pres_l,KoL,PoL,k_surface) - absolute_position(nk,ns,Pres_l,KoL,PoL,k_surface-1) + hR = absolute_position(nk,ns,Pres_r,KoR,PoR,k_surface) - absolute_position(nk,ns,Pres_r,KoR,PoR,k_surface-1) + ! In the case of a layer being unstably stratified, may get a negative thickness. Set the previous position + ! to the current location + if (hL < 0.) then + if ( (KoL(k_surface) 0.) then + hEff(k_surface-1) = 2. * hL * hR / ( hL + hR ) ! Harmonic mean + else + hEff(k_surface-1) = 0. + endif + endif + + enddo neutral_surfaces + +end subroutine find_neutral_surface_positions_discontinuous + +!> Increments the interface which was just connected and also set flags if the bottom is reached +subroutine increment_interface(nk, kl, ki, reached_bottom, searching_this_column, searching_other_column) + integer, intent(in ) :: nk !< Number of vertical levels + integer, intent(inout) :: kl !< Current layer (potentially updated) + integer, intent(inout) :: ki !< Current interface + logical, intent(inout) :: reached_bottom !< Updated when kl == nk and ki == 2 + logical, intent(inout) :: searching_this_column !< Updated when kl == nk and ki == 2 + logical, intent(inout) :: searching_other_column !< Updated when kl == nk and ki == 2 + + if (ki == 1) then + ki = 2 + elseif ((ki == 2) .and. (kl < nk)) then + ki = 1 + kl = kl + 1 + elseif ((kl == nk) .and. (ki==2)) then + reached_bottom = .true. + searching_this_column = .true. + searching_other_column = .false. + else + call MOM_error(FATAL,"Unanticipated eventuality in increment_interface") + endif + + +end subroutine increment_interface + +!> Searches the "other" (searched) column for the position of the neutral surface +subroutine search_other_column_discontinuous(dRhoTopm1, dRhoTop, dRhoBot, Ptop, Pbot, lastP, lastK, kl, ki, & + top_connected, bot_connected, out_P, out_K) + real, intent(in ) :: dRhoTopm1 !< Density difference across previous interface + real, intent(in ) :: dRhoTop !< Density difference across top interface + real, intent(in ) :: dRhoBot !< Density difference across top interface + real, intent(in ) :: Ptop !< Pressure at top interface + real, intent(in ) :: Pbot !< Pressure at bottom interface + real, intent(in ) :: lastP !< Last position connected in the searched column + integer, intent(in ) :: lastK !< Last layer connected in the searched column + integer, intent(in ) :: kl !< Layer in the searched column + integer, intent(in ) :: ki !< Interface of the searched column + logical, dimension(:), intent(inout) :: top_connected !< True if the top interface was pointed to + logical, dimension(:), intent(inout) :: bot_connected !< True if the top interface was pointed to + real, intent( out) :: out_P !< Position within searched column + integer, intent( out) :: out_K !< Layer within searched column + ! Local variables + logical :: search_layer + + search_layer = .true. + ! Bad values to make sure that the particular setup has been processed + out_P = -1. ; out_K = -1 + ! Check if everything in this layer is denser than neutral surface or if at the top of the water column + if ((kl==1 .and. ki==1)) then + if (debug_this_module) write(*,*) "At surface" + out_P = 0. ; out_K = kl + search_layer = .false. + ! Deal with the case where reconstruction is continuous + elseif ( kl>1 ) then + if (dRhoTopm1==dRhoTop .and. dRhoTopm1 == 0. .and. (.not.bot_connected(kl-1)) ) then + out_P = 1 ; out_K = kl-1; + search_layer = .false. + elseif ( (dRhoTopm1 0.) ) then + out_P = 1. ; out_K = kl-1 + search_layer = .false. + endif + endif + + if (search_layer) then + if (dRhoTop > 0.) then + out_P = 0. ; out_K = kl + elseif ( dRhoTop == 0. .and. (.not. top_connected(kl)) ) then + out_P = 0. ; out_K = kl + elseif (dRhoTop >= dRhoBot) then + out_P = 1. ; out_K = kl + else + if (debug_this_module) write(*,*) "Zero crossing point within layer" + out_P = interpolate_for_nondim_position(dRhoTop, Ptop, dRhoBot, Pbot) + out_K = kl + endif + endif + if ( (out_P < 0.) .and. (out_K < 0) ) then + call MOM_error(WARNING, "Unanticipated case in search_other_column_discontinuous") + endif + ! Check to make sure that the layer index is always increasing + if ( (out_K < lastK) .and. lastP==0. .and. out_P == 1. ) then + out_K = lastK ; out_P = 0. + endif + +end subroutine search_other_column_discontinuous !> Converts non-dimensional position within a layer to absolute position (for debugging) -real function absolute_position(n,Pint,Karr,NParr,k_surface) +real function absolute_position(n,ns,Pint,Karr,NParr,k_surface) integer, intent(in) :: n !< Number of levels + integer, intent(in) :: ns !< Number of neutral surfaces real, intent(in) :: Pint(n+1) !< Position of interfaces (Pa) - integer, intent(in) :: Karr(2*n+2) !< Index of interface above position - real, intent(in) :: NParr(2*n+2) !< Non-dimensional position within layer Karr(:) + integer, intent(in) :: Karr(ns) !< Index of interface above position + real, intent(in) :: NParr(ns) !< Non-dimensional position within layer Karr(:) ! Local variables integer :: k_surface, k @@ -904,19 +1370,20 @@ real function absolute_position(n,Pint,Karr,NParr,k_surface) end function absolute_position !> Converts non-dimensional positions within layers to absolute positions (for debugging) -function absolute_positions(n,Pint,Karr,NParr) +function absolute_positions(n,ns,Pint,Karr,NParr) integer, intent(in) :: n !< Number of levels + integer, intent(in) :: ns !< Number of neutral surfaces real, intent(in) :: Pint(n+1) !< Position of interface (Pa) - integer, intent(in) :: Karr(2*n+2) !< Indexes of interfaces about positions - real, intent(in) :: NParr(2*n+2) !< Non-dimensional positions within layers Karr(:) + integer, intent(in) :: Karr(ns) !< Indexes of interfaces about positions + real, intent(in) :: NParr(ns) !< Non-dimensional positions within layers Karr(:) - real, dimension(2*n+2) :: absolute_positions ! Absolute positions (Pa) + real, dimension(ns) :: absolute_positions ! Absolute positions (Pa) ! Local variables integer :: k_surface, k - do k_surface = 1, 2*n+2 - absolute_positions(k_surface) = absolute_position(n,Pint,Karr,NParr,k_surface) + do k_surface = 1, ns + absolute_positions(k_surface) = absolute_position(n,ns,Pint,Karr,NParr,k_surface) enddo end function absolute_positions @@ -952,21 +1419,351 @@ real function interpolate_for_nondim_position(dRhoNeg, Pneg, dRhoPos, Ppos) if ( interpolate_for_nondim_position > 1. ) stop 'interpolate_for_nondim_position: Houston, we have a problem! Pint > Ppos' end function interpolate_for_nondim_position +!> Use root-finding methods to find where dRho = 0, based on the equation of state and the polynomial +!! reconstructions of temperature, salinity. Initial guess is based on the zero crossing of based on linear +!! profiles of dRho, T, and S, between the top and bottom interface. If second derivatives of the EOS are available, +!! it starts with a Newton's method. However, Newton's method is not guaranteed to be bracketed, a check is performed +!! to see if it it diverges outside the interval. In that case (or in the case that second derivatives are not +!! available), Brent's method is used following the implementation found at +!! https://people.sc.fsu.edu/~jburkardt/f_src/brent/brent.f90 +real function refine_nondim_position(max_iter, tolerance, T_ref, S_ref, alpha_ref, beta_ref, P_top, P_bot, deg, & + ppoly_T, ppoly_S, EOS, x0, drho_top, drho_bot, min_bound, force_brent) + integer, intent(in) :: max_iter !< Number of maximum iterations to use + real, intent(in) :: tolerance !< Convergence criterion for delta_rho + real, intent(in) :: T_ref !< Temperature of the neutral surface at the searched from interface + real, intent(in) :: S_ref !< Salinity of the neutral surface at the searched from interface + real, intent(in) :: alpha_ref !< dRho/dT of the neutral surface at the searched from interface + real, intent(in) :: beta_ref !< dRho/dS of the neutral surface at the searched from interface + real, intent(in) :: P_top !< Pressure at the top interface in the layer to be searched + real, intent(in) :: P_bot !< Pressure at the bottom interface in the layer to be searched + integer, intent(in) :: deg !< Order of the polynomimal used for reconstructions + real, dimension(:), intent(in) :: ppoly_T !< Coefficients of the order N polynomial reconstruction of T within + !! the layer to be searched. + real, dimension(:), intent(in) :: ppoly_S !< Coefficients of the order N polynomial reconstruction of T within + !! the layer to be searched. + real, intent(in) :: x0 !< Nondimensional position within the layer where the neutral + !! surface connects. If interpolate_for_nondim_position was + !! previously called, this would be based on linear profile of dRho + real, intent(in) :: drho_top, drho_bot, min_bound + type(EOS_type), pointer :: EOS !< Equation of state structure + logical, optional, intent(in) :: force_brent !< Forces the use of Brent's method instead of Newton's method to find + !! position of neutral surface + ! Local variables + integer :: form_of_EOS + integer :: iter + logical :: do_newton, do_brent + + real :: delta_rho, d_delta_rho_dP ! Terms for the Newton iteration + real :: P_int, P_min ! Interpolated pressure + real :: delta_rho_init, delta_rho_final, x_init + real :: T, S, alpha, beta, alpha_avg, beta_avg + ! Newton's Method variables + real :: dT_dP, dS_dP, delta_T, delta_S, delta_P + real :: dbeta_dS, dbeta_dT, dalpha_dT, dalpha_dS, dbeta_dP, dalpha_dP + ! Brent's Method variables + real :: a, b, c, d, e, f, fa, fb, fc, m, p, q, r, s0, sa, sb, tol, machep + + real :: P_last + logical :: debug = .false. + + delta_P = P_bot-P_top + refine_nondim_position = min_bound + x_init = refine_nondim_position + + call extract_member_EOS(EOS, form_of_EOS = form_of_EOS) + do_newton = (form_of_EOS == EOS_LINEAR) .or. (form_of_EOS == EOS_TEOS10) .or. (form_of_EOS == EOS_WRIGHT) + do_brent = .not. do_newton + if (present(force_brent)) then + do_newton = .not. force_brent + do_brent = force_brent + endif -!> Returns a single column of neutral diffusion fluxes of a tracer. -subroutine neutral_surface_flux(nk, hl, hr, Tl, Tr, PiL, PiR, KoL, KoR, hEff, Flx) - integer, intent(in) :: nk !< Number of levels - real, dimension(nk), intent(in) :: hl !< Left-column layer thickness (Pa) - real, dimension(nk), intent(in) :: hr !< Right-column layer thickness (Pa) - real, dimension(nk), intent(in) :: Tl !< Left-column layer tracer (conc, e.g. degC) - real, dimension(nk), intent(in) :: Tr !< Right-column layer tracer (conc, e.g. degC) - real, dimension(2*nk+2), intent(in) :: PiL !< Fractional position of neutral surface within layer KoL of left column - real, dimension(2*nk+2), intent(in) :: PiR !< Fractional position of neutral surface within layer KoR of right column - integer, dimension(2*nk+2), intent(in) :: KoL !< Index of first left interface above neutral surface - integer, dimension(2*nk+2), intent(in) :: KoR !< Index of first right interface above neutral surface - real, dimension(2*nk+1), intent(in) :: hEff !< Effective thickness between two neutral surfaces (Pa) - real, dimension(2*nk+1), intent(inout) :: Flx !< Flux of tracer between pairs of neutral layers (conc H) + ! Check to make sure that a root exists between the minimum bound and the bottom of the layer + call calc_delta_rho(deg, T_ref, S_ref, alpha_ref, beta_ref, P_top, P_bot, ppoly_T, ppoly_S, refine_nondim_position, & + EOS, delta_rho) + delta_rho_init = delta_rho + if ( SIGN(1.,delta_rho) == SIGN(1.,drho_bot) ) then + ! Return the position of min_bound if closer to 0 than drho_bot + if (ABS(delta_rho) < ABS(drho_bot)) then + refine_nondim_position = min_bound + else + refine_nondim_position = 1. + endif + do_newton = .false. ; do_brent = .false. + endif + if (debug) then + write (*,*) "------" + write (*,*) "Starting delta_rho: ", delta_rho + endif + + ! For now only linear, Wright, and TEOS-10 equations of state have functions providing second derivatives and + ! thus can use Newton's method for the equation of state + if (do_newton) then + ! Set lower bound of pressure + P_min = P_top*(1.-min_bound) + P_bot*(min_bound) + ! Iterate over Newton's method for the function: x0 = x0 - delta_rho/d_delta_rho_dP + do iter = 1, max_iter + ! Evaluate delta_rho(x0) + call calc_delta_rho(deg, T_ref, S_ref, alpha_ref, beta_ref, P_top, P_bot, ppoly_T, ppoly_S, refine_nondim_position, EOS, & + delta_rho, P_int, T, S, alpha_avg, beta_avg, delta_T, delta_S) + ! Check for convergence + if (ABS(delta_rho) <= tolerance) then + do_brent = .false. + exit + endif + ! Evaluate total derivative of delta_rho + call calculate_density_second_derivs( T, S, P_int, dbeta_dS, dbeta_dT, dalpha_dT, dbeta_dP, dalpha_dP, EOS ) + dalpha_dS = dbeta_dT ! Cross derivatives are identicial + ! By chain rule dT_dP= (dT_dz)*(dz/dP) = dT_dz / (Pbot-Ptop) + dT_dP = first_derivative_polynomial( ppoly_T, deg+1, refine_nondim_position ) / delta_P + dS_dP = first_derivative_polynomial( ppoly_S, deg+1, refine_nondim_position ) / delta_P + ! Total derivative of d_delta_rho wrt P + d_delta_rho_dP = 0.5*( delta_S*(dS_dP*dbeta_dS + dT_dP*dbeta_dT + dbeta_dP) + & + ( delta_T*(dS_dP*dalpha_dS + dT_dP*dalpha_dT + dalpha_dP))) + & + dS_dP*beta_avg + dT_dP*alpha_avg + if (d_delta_rho_dP == 0.) then + do_brent = .true. + exit + endif + ! Newton step update + P_last = P_int + P_int = P_int - (delta_rho / d_delta_rho_dP) + if (P_int < P_min .or. P_int > P_bot) then + if (debug) then + write (*,*) "Iteration: ", iter + write (*,*) "delta_rho, d_delta_rho_dP: ", delta_rho, d_delta_rho_dP + write (*,*) "T, T Poly Coeffs: ", T, ppoly_T + write (*,*) "S, S Poly Coeffs: ", S, ppoly_S + write (*,*) "T_ref, alpha_ref: ", T_ref, alpha_ref + write (*,*) "S_ref, beta_ref : ", S_ref, beta_ref + write (*,*) "P, dT_dP, dS_dP:", P_int, dT_dP, dS_dP + write (*,*) "dRhoTop, dRhoBot:", drho_top, drho_bot + write (*,*) "x0: ", x0 + write (*,*) "refine_nondim_position: ", refine_nondim_position + write (*,*) + endif +! call MOM_error(WARNING, "Step went out of bounds") + ! Switch to Brent's method by setting the converged flag to false + do_brent = .true. + ! Reset to first guess if already diverged + if (ABS(delta_rho_init)1.) then + if (debug) then + write (*,*) "T, T Poly Coeffs: ", T, ppoly_T + write (*,*) "S, S Poly Coeffs: ", S, ppoly_S + write (*,*) "T_ref, alpha_ref: ", T_ref, alpha_ref + write (*,*) "S_ref, beta_ref : ", S_ref, beta_ref + write (*,*) "P, dT_dP, dS_dP:", P_int, dT_dP, dS_dP + write (*,*) "x0: ", x0 + write (*,*) "refine_nondim_position: ", refine_nondim_position + endif + call MOM_error(WARNING, "refine_nondim_position>1.") + refine_nondim_position = x0 + endif + + if (refine_nondim_position<0.) then + if (debug) then + write (*,*) "T, T Poly Coeffs: ", T, ppoly_T + write (*,*) "S, S Poly Coeffs: ", S, ppoly_S + write (*,*) "T_ref, alpha_ref: ", T_ref, alpha_ref + write (*,*) "S_ref, beta_ref : ", S_ref, beta_ref + write (*,*) "dT_dP, dS_dP:", dT_dP, dS_dP + write (*,*) "x0: ", x0 + write (*,*) "refine_nondim_position: ", refine_nondim_position + endif + call MOM_error(WARNING, "refine_nondim_position<0.") + refine_nondim_position = x0 + endif + + call calc_delta_rho(deg, T_ref, S_ref, alpha_ref, beta_ref, P_top, P_bot, ppoly_T, ppoly_S, & + refine_nondim_position, EOS, delta_rho) + if (ABS(delta_rho)>=ABS(delta_rho_init)) then + refine_nondim_position = x_init + endif + + + if (debug) then + call calc_delta_rho(deg, T_ref, S_ref, alpha_ref, beta_ref, P_top, P_bot, ppoly_T, ppoly_S, & + refine_nondim_position, EOS, delta_rho) + write (*,*) "End delta_rho: ", delta_rho + write (*,*) "x0, delta_x: ", x0, refine_nondim_position-x0 + write (*,*) "Iterations: ", iter + write (*,*) "******" + endif + +end function refine_nondim_position + +!> Calculate the difference in neutral density between a reference T, S, alpha, and beta +!! and a point on the polynomial reconstructions of T, S +subroutine calc_delta_rho(deg, T_ref, S_ref, alpha_ref, beta_ref, P_top, P_bot, ppoly_T, ppoly_S, x0, EOS, & + delta_rho, P_out, T_out, S_out, alpha_avg_out, beta_avg_out, delta_T_out, delta_S_out) + integer, intent(in) :: deg !< Degree of polynomial reconstruction + real, intent(in) :: T_ref !< Temperature at reference surface + real, intent(in) :: S_ref !< Salinity at reference surface + real, intent(in) :: alpha_ref !< dRho/dT at reference surface + real, intent(in) :: beta_ref !< dRho/dS at reference surface + real, intent(in) :: P_top !< Pressure (Pa) at top interface of layer to be searched + real, intent(in) :: P_bot !< Pressure (Pa) at bottom interface + real, dimension(deg+1), intent(in) :: ppoly_T !< Coefficients of T reconstruction + real, dimension(deg+1), intent(in) :: ppoly_S !< Coefficients of S reconstruciton + real, intent(in) :: x0 !< Nondimensional position to evaluate + type(EOS_type), pointer :: EOS !< Equation of state structure + real, intent(out) :: delta_rho + real, optional, intent(out) :: P_out !< Pressure at point x0 + real, optional, intent(out) :: T_out !< Temperature at point x0 + real, optional, intent(out) :: S_out !< Salinity at point x0 + real, optional, intent(out) :: alpha_avg_out !< Average of alpha between reference and x0 + real, optional, intent(out) :: beta_avg_out !< Average of beta between reference and x0 + real, optional, intent(out) :: delta_T_out !< Difference in temperature between reference and x0 + real, optional, intent(out) :: delta_S_out !< Difference in salinity between reference and x0 + + real :: alpha, beta, alpha_avg, beta_avg, P_int, T, S, delta_T, delta_S + P_int = (1. - x0)*P_top + x0*P_bot + T = evaluation_polynomial( ppoly_T, deg+1, x0 ) + S = evaluation_polynomial( ppoly_S, deg+1, x0 ) + call calculate_density_derivs( T, S, P_int, alpha, beta, EOS ) + + ! Calculate the f(P) term for Newton's method + alpha_avg = 0.5*( alpha + alpha_ref ) + beta_avg = 0.5*( beta + beta_ref ) + delta_T = T - T_ref + delta_S = S - S_ref + delta_rho = alpha_avg*delta_T + beta_avg*delta_S + + ! If doing a Newton step, these quantities are needed, otherwise they can just be optional + if (present(P_out)) P_out = P_int + if (present(T_out)) T_out = T + if (present(S_out)) S_out = S + if (present(alpha_avg_out)) alpha_avg_out = alpha_avg + if (present(beta_avg_out)) beta_avg_out = beta_avg + if (present(delta_T_out)) delta_T_out = delta_T + if (present(delta_S_out)) delta_S_out = delta_S + +end subroutine calc_delta_rho + +!> Returns a single column of neutral diffusion fluxes of a tracer. +subroutine neutral_surface_flux(nk, nsurf, deg, hl, hr, Tl, Tr, PiL, PiR, KoL, KoR, hEff, Flx, continuous, remap_CS) + integer, intent(in) :: nk !< Number of levels + integer, intent(in) :: nsurf !< Number of neutral surfaces + integer, intent(in) :: deg !< Degree of polynomial reconstructions + real, dimension(nk), intent(in) :: hl !< Left-column layer thickness (Pa) + real, dimension(nk), intent(in) :: hr !< Right-column layer thickness (Pa) + real, dimension(nk), intent(in) :: Tl !< Left-column layer tracer (conc, e.g. degC) + real, dimension(nk), intent(in) :: Tr !< Right-column layer tracer (conc, e.g. degC) + real, dimension(nsurf), intent(in) :: PiL !< Fractional position of neutral surface + !! within layer KoL of left column + real, dimension(nsurf), intent(in) :: PiR !< Fractional position of neutral surface + !! within layer KoR of right column + integer, dimension(nsurf), intent(in) :: KoL !< Index of first left interface above neutral surface + integer, dimension(nsurf), intent(in) :: KoR !< Index of first right interface above neutral surface + real, dimension(nsurf-1), intent(in) :: hEff !< Effective thickness between two neutral surfaces (Pa) + real, dimension(nsurf-1), intent(inout) :: Flx !< Flux of tracer between pairs of neutral layers (conc H) + logical, intent(in) :: continuous !< True if using continuous reconstruction + type(remapping_CS), optional, intent(in) :: remap_CS ! Local variables integer :: k_sublayer, klb, klt, krb, krt, k real :: T_right_top, T_right_bottom, T_right_layer @@ -978,57 +1775,73 @@ subroutine neutral_surface_flux(nk, hl, hr, Tl, Tr, PiL, PiR, KoL, KoR, hEff, Fl real, dimension(nk) :: aR_l !< Left-column right edge value of tracer (conc, e.g. degC) real, dimension(nk) :: aL_r !< Right-column left edge value of tracer (conc, e.g. degC) real, dimension(nk) :: aR_r !< Right-column right edge value of tracer (conc, e.g. degC) - - call interface_scalar(nk, hl, Tl, Til, 2) - call interface_scalar(nk, hr, Tr, Tir, 2) + ! Discontinuous reconstruction + integer :: iMethod + real, dimension(nk,2) :: Tid_l !< Left-column interface tracer (conc, e.g. degC) + real, dimension(nk,2) :: Tid_r !< Right-column interface tracer (conc, e.g. degC) + real, dimension(nk,deg+1) :: ppoly_r_coeffs_l + real, dimension(nk,deg+1) :: ppoly_r_coeffs_r + real, dimension(nk,deg+1) :: ppoly_r_S_l + real, dimension(nk,deg+1) :: ppoly_r_S_r ! Setup reconstruction edge values - do k = 1, nk - aL_l(k) = Til(k) - aR_l(k) = Til(k+1) - if ( signum(1., aR_l(k) - Tl(k))*signum(1., Tl(k) - aL_l(k)) <= 0.0 ) then - aL_l(k) = Tl(k) - aR_l(k) = Tl(k) - elseif ( sign(3., aR_l(k) - aL_l(k)) * ( (Tl(k) - aL_l(k)) + (Tl(k) - aR_l(k))) > abs(aR_l(k) - aL_l(k)) ) then - aL_l(k) = Tl(k) + 2.0 * ( Tl(k) - aR_l(k) ) - elseif ( sign(3., aR_l(k) - aL_l(k)) * ( (Tl(k) - aL_l(k)) + (Tl(k) - aR_l(k))) < -abs(aR_l(k) - aL_l(k)) ) then - aR_l(k) = Tl(k) + 2.0 * ( Tl(k) - aL_l(k) ) - endif - aL_r(k) = Tir(k) - aR_r(k) = Tir(k+1) - if ( signum(1., aR_r(k) - Tr(k))*signum(1., Tr(k) - aL_r(k)) <= 0.0 ) then - aL_r(k) = Tr(k) - aR_r(k) = Tr(k) - elseif ( sign(3., aR_r(k) - aL_r(k)) * ( (Tr(k) - aL_r(k)) + (Tr(k) - aR_r(k))) > abs(aR_r(k) - aL_r(k)) ) then - aL_r(k) = Tr(k) + 2.0 * ( Tr(k) - aR_r(k) ) - elseif ( sign(3., aR_r(k) - aL_r(k)) * ( (Tr(k) - aL_r(k)) + (Tr(k) - aR_r(k))) < -abs(aR_r(k) - aL_r(k)) ) then - aR_r(k) = Tr(k) + 2.0 * ( Tr(k) - aL_r(k) ) - endif - enddo + if (continuous) then + call interface_scalar(nk, hl, Tl, Til, 2) + call interface_scalar(nk, hr, Tr, Tir, 2) + call ppm_left_right_edge_values(nk, Tl, Til, aL_l, aR_l) + call ppm_left_right_edge_values(nk, Tr, Tir, aL_r, aR_r) + else + ppoly_r_coeffs_l(:,:) = 0. + ppoly_r_coeffs_r(:,:) = 0. + Tid_l(:,:) = 0. + Tid_r(:,:) = 0. - do k_sublayer = 1, 2*nk+1 + call build_reconstructions_1d( remap_CS, nk, hl, Tl, ppoly_r_coeffs_l, Tid_l, ppoly_r_S_l, iMethod ) + call build_reconstructions_1d( remap_CS, nk, hr, Tr, ppoly_r_coeffs_r, Tid_r, ppoly_r_S_r, iMethod ) + endif + + do k_sublayer = 1, nsurf-1 if (hEff(k_sublayer) == 0.) then Flx(k_sublayer) = 0. else - - klb = KoL(k_sublayer+1) - T_left_bottom = ( 1. - PiL(k_sublayer+1) ) * Til(klb) + PiL(k_sublayer+1) * Til(klb+1) - - klt = KoL(k_sublayer) - T_left_top = ( 1. - PiL(k_sublayer) ) * Til(klt) + PiL(k_sublayer) * Til(klt+1) - - T_left_layer = ppm_ave(PiL(k_sublayer), PiL(k_sublayer+1) + real(klb-klt), & - aL_l(klt), aR_l(klt), Tl(klt)) - - krb = KoR(k_sublayer+1) - T_right_bottom = ( 1. - PiR(k_sublayer+1) ) * Tir(krb) + PiR(k_sublayer+1) * Tir(krb+1) - - krt = KoR(k_sublayer) - T_right_top = ( 1. - PiR(k_sublayer) ) * Tir(krt) + PiR(k_sublayer) * Tir(krt+1) - - T_right_layer = ppm_ave(PiR(k_sublayer), PiR(k_sublayer+1) + real(krb-krt), & - aL_r(krt), aR_r(krt), Tr(krt)) - + if (continuous) then + klb = KoL(k_sublayer+1) + T_left_bottom = ( 1. - PiL(k_sublayer+1) ) * Til(klb) + PiL(k_sublayer+1) * Til(klb+1) + klt = KoL(k_sublayer) + T_left_top = ( 1. - PiL(k_sublayer) ) * Til(klt) + PiL(k_sublayer) * Til(klt+1) + T_left_layer = ppm_ave(PiL(k_sublayer), PiL(k_sublayer+1) + real(klb-klt), & + aL_l(klt), aR_l(klt), Tl(klt)) + + krb = KoR(k_sublayer+1) + T_right_bottom = ( 1. - PiR(k_sublayer+1) ) * Tir(krb) + PiR(k_sublayer+1) * Tir(krb+1) + krt = KoR(k_sublayer) + T_right_top = ( 1. - PiR(k_sublayer) ) * Tir(krt) + PiR(k_sublayer) * Tir(krt+1) + T_right_layer = ppm_ave(PiR(k_sublayer), PiR(k_sublayer+1) + real(krb-krt), & + aL_r(krt), aR_r(krt), Tr(krt)) + else ! Discontinuous reconstruction + klb = KoL(k_sublayer+1) + klt = KoL(k_sublayer) + if (klt .ne. klb) then + call MOM_error(WARNING, "Neutral surfaces span more than one layer") + Flx(k_sublayer) = 0. + cycle + endif + T_left_bottom = evaluation_polynomial( ppoly_r_coeffs_l(klb,:), deg+1, PiL(k_sublayer+1)) + T_left_top = evaluation_polynomial( ppoly_r_coeffs_l(klt,:), deg+1, PiL(k_sublayer)) + T_left_layer = average_value_ppoly(nk, Tl, Tid_l, ppoly_r_coeffs_l, iMethod, klb, & + PiL(k_sublayer), PiL(k_sublayer+1)) + krb = KoR(k_sublayer+1) + krt = KoR(k_sublayer) + if (krt .ne. krb) then + call MOM_error(WARNING, "Neutral surfaces span more than one layer") + Flx(k_sublayer) = 0. + cycle + endif + T_right_bottom = evaluation_polynomial( ppoly_r_coeffs_r(krb,:), deg+1, PiR(k_sublayer+1)) + T_right_top = evaluation_polynomial( ppoly_r_coeffs_r(krt,:), deg+1, PiR(k_sublayer)) + T_right_layer = average_value_ppoly(nk, Tr, Tid_r, ppoly_r_coeffs_r, iMethod, krb, & + PiR(k_sublayer), PiR(k_sublayer+1)) + endif dT_top = T_right_top - T_left_top dT_bottom = T_right_bottom - T_left_bottom dT_ave = 0.5 * ( dT_top + dT_bottom ) @@ -1044,8 +1857,42 @@ subroutine neutral_surface_flux(nk, hl, hr, Tl, Tr, PiL, PiR, KoL, KoR, hEff, Fl end subroutine neutral_surface_flux +!> Discontinuous PPM reconstructions of the left/right edge values within a cell +subroutine ppm_left_right_edge_values(nk, Tl, Ti, aL, aR) + integer, intent(in) :: nk !< Number of levels + real, dimension(nk), intent(in) :: Tl !< Layer tracer (conc, e.g. degC) + real, dimension(nk+1), intent(in) :: Ti !< Interface tracer (conc, e.g. degC) + real, dimension(nk), intent(inout) :: aL !< Left edge value of tracer (conc, e.g. degC) + real, dimension(nk), intent(inout) :: aR !< Right edge value of tracer (conc, e.g. degC) + + integer :: k + ! Setup reconstruction edge values + do k = 1, nk + aL(k) = Ti(k) + aR(k) = Ti(k+1) + if ( signum(1., aR(k) - Tl(k))*signum(1., Tl(k) - aL(k)) <= 0.0 ) then + aL(k) = Tl(k) + aR(k) = Tl(k) + elseif ( sign(3., aR(k) - aL(k)) * ( (Tl(k) - aL(k)) + (Tl(k) - aR(k))) > abs(aR(k) - aL(k)) ) then + aL(k) = Tl(k) + 2.0 * ( Tl(k) - aR(k) ) + elseif ( sign(3., aR(k) - aL(k)) * ( (Tl(k) - aL(k)) + (Tl(k) - aR(k))) < -abs(aR(k) - aL(k)) ) then + aR(k) = Tl(k) + 2.0 * ( Tl(k) - aL(k) ) + endif + enddo +end subroutine ppm_left_right_edge_values + !> Returns true if unit tests of neutral_diffusion functions fail. Otherwise returns false. logical function neutral_diffusion_unit_tests(verbose) + logical, intent(in) :: verbose + + neutral_diffusion_unit_tests = .false. .or. & + ndiff_unit_tests_continuous(verbose) .or. ndiff_unit_tests_discontinuous(verbose) + + +end function neutral_diffusion_unit_tests + +!> Returns true if unit tests of neutral_diffusion functions fail. Otherwise returns false. +logical function ndiff_unit_tests_continuous(verbose) logical, intent(in) :: verbose !< It true, write results to stdout ! Local variables integer, parameter :: nk = 4 @@ -1062,105 +1909,105 @@ logical function neutral_diffusion_unit_tests(verbose) v = verbose - neutral_diffusion_unit_tests = .false. ! Normally return false - write(*,*) '==== MOM_neutral_diffusion: neutral_diffusion_unit_tests =' - - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fv_diff(v,1.,1.,1., 0.,1.,2., 1., 'FV: Straight line on uniform grid') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fv_diff(v,1.,1.,0., 0.,4.,8., 7., 'FV: Vanished right cell') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fv_diff(v,0.,1.,1., 0.,4.,8., 7., 'FV: Vanished left cell') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fv_diff(v,1.,2.,4., 0.,3.,9., 4., 'FV: Stretched grid') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fv_diff(v,2.,0.,2., 0.,1.,2., 0., 'FV: Vanished middle cell') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fv_diff(v,0.,1.,0., 0.,1.,2., 2., 'FV: Vanished on both sides') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fv_diff(v,1.,0.,0., 0.,1.,2., 0., 'FV: Two vanished cell sides') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fv_diff(v,0.,0.,0., 0.,1.,2., 0., 'FV: All vanished cells') - - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fvlsq_slope(v,1.,1.,1., 0.,1.,2., 1., 'LSQ: Straight line on uniform grid') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fvlsq_slope(v,1.,1.,0., 0.,1.,2., 1., 'LSQ: Vanished right cell') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fvlsq_slope(v,0.,1.,1., 0.,1.,2., 1., 'LSQ: Vanished left cell') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fvlsq_slope(v,1.,2.,4., 0.,3.,9., 2., 'LSQ: Stretched grid') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fvlsq_slope(v,1.,0.,1., 0.,1.,2., 2., 'LSQ: Vanished middle cell') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fvlsq_slope(v,0.,1.,0., 0.,1.,2., 0., 'LSQ: Vanished on both sides') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fvlsq_slope(v,1.,0.,0., 0.,1.,2., 0., 'LSQ: Two vanished cell sides') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fvlsq_slope(v,0.,0.,0., 0.,1.,2., 0., 'LSQ: All vanished cells') + ndiff_unit_tests_continuous = .false. ! Normally return false + write(*,*) '==== MOM_neutral_diffusion: ndiff_unit_tests_continuous =' + + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_fv_diff(v,1.,1.,1., 0.,1.,2., 1., 'FV: Straight line on uniform grid') + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_fv_diff(v,1.,1.,0., 0.,4.,8., 7., 'FV: Vanished right cell') + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_fv_diff(v,0.,1.,1., 0.,4.,8., 7., 'FV: Vanished left cell') + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_fv_diff(v,1.,2.,4., 0.,3.,9., 4., 'FV: Stretched grid') + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_fv_diff(v,2.,0.,2., 0.,1.,2., 0., 'FV: Vanished middle cell') + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_fv_diff(v,0.,1.,0., 0.,1.,2., 2., 'FV: Vanished on both sides') + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_fv_diff(v,1.,0.,0., 0.,1.,2., 0., 'FV: Two vanished cell sides') + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_fv_diff(v,0.,0.,0., 0.,1.,2., 0., 'FV: All vanished cells') + + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_fvlsq_slope(v,1.,1.,1., 0.,1.,2., 1., 'LSQ: Straight line on uniform grid') + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_fvlsq_slope(v,1.,1.,0., 0.,1.,2., 1., 'LSQ: Vanished right cell') + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_fvlsq_slope(v,0.,1.,1., 0.,1.,2., 1., 'LSQ: Vanished left cell') + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_fvlsq_slope(v,1.,2.,4., 0.,3.,9., 2., 'LSQ: Stretched grid') + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_fvlsq_slope(v,1.,0.,1., 0.,1.,2., 2., 'LSQ: Vanished middle cell') + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_fvlsq_slope(v,0.,1.,0., 0.,1.,2., 0., 'LSQ: Vanished on both sides') + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_fvlsq_slope(v,1.,0.,0., 0.,1.,2., 0., 'LSQ: Two vanished cell sides') + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_fvlsq_slope(v,0.,0.,0., 0.,1.,2., 0., 'LSQ: All vanished cells') call interface_scalar(4, (/10.,10.,10.,10./), (/24.,18.,12.,6./), Tio, 1) - !neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_data1d(5, Tio, (/27.,21.,15.,9.,3./), 'Linear profile, interface temperatures') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_data1d(v,5, Tio, (/24.,22.5,15.,7.5,6./), 'Linear profile, linear interface temperatures') + !ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_data1d(5, Tio, (/27.,21.,15.,9.,3./), 'Linear profile, interface temperatures') + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_data1d(v,5, Tio, (/24.,22.5,15.,7.5,6./), 'Linear profile, linear interface temperatures') call interface_scalar(4, (/10.,10.,10.,10./), (/24.,18.,12.,6./), Tio, 2) - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_data1d(v,5, Tio, (/24.,22.,15.,8.,6./), 'Linear profile, PPM interface temperatures') - - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(v,-1.0, 0., 1.0, 1.0, 0.5, 'Check mid-point') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(v, 0.0, 0., 1.0, 1.0, 0.0, 'Check bottom') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(v, 0.1, 0., 1.1, 1.0, 0.0, 'Check below') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(v,-1.0, 0., 0.0, 1.0, 1.0, 'Check top') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(v,-1.0, 0., -0.1, 1.0, 1.0, 'Check above') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(v,-1.0, 0., 3.0, 1.0, 0.25, 'Check 1/4') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(v,-3.0, 0., 1.0, 1.0, 0.75, 'Check 3/4') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(v, 1.0, 0., 1.0, 1.0, 0.0, 'Check dRho=0 below') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(v,-1.0, 0., -1.0, 1.0, 1.0, 'Check dRho=0 above') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(v, 0.0, 0., 0.0, 1.0, 0.5, 'Check dRho=0 mid') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(v,-2.0, .5, 5.0, 0.5, 0.5, 'Check dP=0') + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_data1d(v,5, Tio, (/24.,22.,15.,8.,6./), 'Linear profile, PPM interface temperatures') + + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_ifndp(v,-1.0, 0., 1.0, 1.0, 0.5, 'Check mid-point') + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_ifndp(v, 0.0, 0., 1.0, 1.0, 0.0, 'Check bottom') + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_ifndp(v, 0.1, 0., 1.1, 1.0, 0.0, 'Check below') + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_ifndp(v,-1.0, 0., 0.0, 1.0, 1.0, 'Check top') + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_ifndp(v,-1.0, 0., -0.1, 1.0, 1.0, 'Check above') + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_ifndp(v,-1.0, 0., 3.0, 1.0, 0.25, 'Check 1/4') + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_ifndp(v,-3.0, 0., 1.0, 1.0, 0.75, 'Check 3/4') + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_ifndp(v, 1.0, 0., 1.0, 1.0, 0.0, 'Check dRho=0 below') + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_ifndp(v,-1.0, 0., -1.0, 1.0, 1.0, 'Check dRho=0 above') + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_ifndp(v, 0.0, 0., 0.0, 1.0, 0.5, 'Check dRho=0 mid') + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_ifndp(v,-2.0, .5, 5.0, 0.5, 0.5, 'Check dP=0') ! Identical columns - call find_neutral_surface_positions(3, & + call find_neutral_surface_positions_continuous(3, & (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Left positions, T and S (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Right positions, T and S (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS PiLRo, PiRLo, KoL, KoR, hEff) - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(v, 3, KoL, KoR, PiLRo, PiRLo, hEff, & + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, KoL, KoR, PiLRo, PiRLo, hEff, & (/1,1,2,2,3,3,3,3/), & ! KoL (/1,1,2,2,3,3,3,3/), & ! KoR (/0.,0.,0.,0.,0.,0.,1.,1./), & ! pL (/0.,0.,0.,0.,0.,0.,1.,1./), & ! pR (/0.,10.,0.,10.,0.,10.,0./), & ! hEff 'Identical columns') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_data1d(v, 8, & - absolute_positions(3, (/0.,10.,20.,30./), KoL, PiLRo), & + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_data1d(v, 8, & + absolute_positions(3, 8, (/0.,10.,20.,30./), KoL, PiLRo), & (/0.,0.,10.,10.,20.,20.,30.,30./), '... left positions') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_data1d(v, 8, & - absolute_positions(3, (/0.,10.,20.,30./), KoR, PiRLo), & + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_data1d(v, 8, & + absolute_positions(3, 8, (/0.,10.,20.,30./), KoR, PiRLo), & (/0.,0.,10.,10.,20.,20.,30.,30./), '... right positions') - call neutral_surface_flux(3, (/10.,10.,10./), (/10.,10.,10./), & ! nk, hL, hR + call neutral_surface_flux(3, 2*3+2, 2, (/10.,10.,10./), (/10.,10.,10./), & ! nk, hL, hR (/20.,16.,12./), (/20.,16.,12./), & ! Tl, Tr - PiLRo, PiRLo, KoL, KoR, hEff, Flx) - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_data1d(v, 7, Flx, & + PiLRo, PiRLo, KoL, KoR, hEff, Flx, .true.) + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_data1d(v, 7, Flx, & (/0.,0.,0.,0.,0.,0.,0./), 'Identical columns, rho flux (=0)') - call neutral_surface_flux(3, (/10.,10.,10./), (/10.,10.,10./), & ! nk, hL, hR + call neutral_surface_flux(3, 2*3+2, 2, (/10.,10.,10./), (/10.,10.,10./), & ! nk, hL, hR (/-1.,-1.,-1./), (/1.,1.,1./), & ! Sl, Sr - PiLRo, PiRLo, KoL, KoR, hEff, Flx) - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_data1d(v, 7, Flx, & + PiLRo, PiRLo, KoL, KoR, hEff, Flx, .true.) + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_data1d(v, 7, Flx, & (/0.,20.,0.,20.,0.,20.,0./), 'Identical columns, S flux') ! Right column slightly cooler than left - call find_neutral_surface_positions(3, & + call find_neutral_surface_positions_continuous(3, & (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Left positions, T and S (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS (/0.,10.,20.,30./), (/20.,16.,12.,8./), (/0.,0.,0.,0./), & ! Right positions, T and S (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS PiLRo, PiRLo, KoL, KoR, hEff) - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(v, 3, KoL, KoR, PiLRo, PiRLo, hEff, & + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, KoL, KoR, PiLRo, PiRLo, hEff, & (/1,1,2,2,3,3,3,3/), & ! kL (/1,1,1,2,2,3,3,3/), & ! kR (/0.,0.5,0.,0.5,0.,0.5,1.,1./), & ! pL (/0.,0.,0.5,0.,0.5,0.,0.5,1./), & ! pR (/0.,5.,5.,5.,5.,5.,0./), & ! hEff 'Right column slightly cooler') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_data1d(v, 8, & - absolute_positions(3, (/0.,10.,20.,30./), KoL, PiLRo), & + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_data1d(v, 8, & + absolute_positions(3, 8, (/0.,10.,20.,30./), KoL, PiLRo), & (/0.,5.,10.,15.,20.,25.,30.,30./), '... left positions') - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_data1d(v, 8, & - absolute_positions(3, (/0.,10.,20.,30./), KoR, PiRLo), & + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_data1d(v, 8, & + absolute_positions(3, 8, (/0.,10.,20.,30./), KoR, PiRLo), & (/0.,0.,5.,10.,15.,20.,25.,30./), '... right positions') ! Right column slightly warmer than left - call find_neutral_surface_positions(3, & + call find_neutral_surface_positions_continuous(3, & (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Left positions, T and S (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS (/0.,10.,20.,30./), (/24.,20.,16.,12./), (/0.,0.,0.,0./), & ! Right positions, T and S (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS PiLRo, PiRLo, KoL, KoR, hEff) - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(v, 3, KoL, KoR, PiLRo, PiRLo, hEff, & + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, KoL, KoR, PiLRo, PiRLo, hEff, & (/1,1,1,2,2,3,3,3/), & ! kL (/1,1,2,2,3,3,3,3/), & ! kR (/0.,0.,0.5,0.,0.5,0.,0.5,1./), & ! pL @@ -1169,13 +2016,13 @@ logical function neutral_diffusion_unit_tests(verbose) 'Right column slightly warmer') ! Right column somewhat cooler than left - call find_neutral_surface_positions(3, & + call find_neutral_surface_positions_continuous(3, & (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Left positions, T and S (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS (/0.,10.,20.,30./), (/16.,12.,8.,4./), (/0.,0.,0.,0./), & ! Right positions, T and S (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS PiLRo, PiRLo, KoL, KoR, hEff) - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(v, 3, KoL, KoR, PiLRo, PiRLo, hEff, & + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, KoL, KoR, PiLRo, PiRLo, hEff, & (/1,2,2,3,3,3,3,3/), & ! kL (/1,1,1,1,2,2,3,3/), & ! kR (/0.,0.,0.5,0.,0.5,1.,1.,1./), & ! pL @@ -1184,13 +2031,13 @@ logical function neutral_diffusion_unit_tests(verbose) 'Right column somewhat cooler') ! Right column much colder than left with no overlap - call find_neutral_surface_positions(3, & + call find_neutral_surface_positions_continuous(3, & (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Left positions, T and S (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS (/0.,10.,20.,30./), (/9.,7.,5.,3./), (/0.,0.,0.,0./), & ! Right positions, T and S (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS PiLRo, PiRLo, KoL, KoR, hEff) - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(v, 3, KoL, KoR, PiLRo, PiRLo, hEff, & + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, KoL, KoR, PiLRo, PiRLo, hEff, & (/1,2,3,3,3,3,3,3/), & ! kL (/1,1,1,1,1,2,3,3/), & ! kR (/0.,0.,0.,1.,1.,1.,1.,1./), & ! pL @@ -1199,13 +2046,13 @@ logical function neutral_diffusion_unit_tests(verbose) 'Right column much cooler') ! Right column with mixed layer - call find_neutral_surface_positions(3, & + call find_neutral_surface_positions_continuous(3, & (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Left positions, T and S (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS (/0.,10.,20.,30./), (/14.,14.,10.,2./), (/0.,0.,0.,0./), & ! Right positions, T and S (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS PiLRo, PiRLo, KoL, KoR, hEff) - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(v, 3, KoL, KoR, PiLRo, PiRLo, hEff, & + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, KoL, KoR, PiLRo, PiRLo, hEff, & (/1,2,3,3,3,3,3,3/), & ! kL (/1,1,1,1,2,3,3,3/), & ! kR (/0.,0.,0.,0.,0.,1.,1.,1./), & ! pL @@ -1214,13 +2061,13 @@ logical function neutral_diffusion_unit_tests(verbose) 'Right column with mixed layer') ! Identical columns with mixed layer - call find_neutral_surface_positions(3, & + call find_neutral_surface_positions_continuous(3, & (/0.,10.,20.,30./), (/14.,14.,10.,2./), (/0.,0.,0.,0./), & ! Left positions, T and S (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS (/0.,10.,20.,30./), (/14.,14.,10.,2./), (/0.,0.,0.,0./), & ! Right positions, T and S (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS PiLRo, PiRLo, KoL, KoR, hEff) - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(v, 3, KoL, KoR, PiLRo, PiRLo, hEff, & + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, KoL, KoR, PiLRo, PiRLo, hEff, & (/1,1,2,2,3,3,3,3/), & ! kL (/1,1,2,2,3,3,3,3/), & ! kR (/0.,0.,0.,0.,0.,0.,1.,1./), & ! pL @@ -1229,13 +2076,13 @@ logical function neutral_diffusion_unit_tests(verbose) 'Indentical columns with mixed layer') ! Right column with unstable mixed layer - call find_neutral_surface_positions(3, & + call find_neutral_surface_positions_continuous(3, & (/0.,10.,20.,30./), (/14.,14.,10.,2./), (/0.,0.,0.,0./), & ! Left positions, T and S (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS (/0.,10.,20.,30./), (/10.,14.,12.,4./), (/0.,0.,0.,0./), & ! Right positions, T and S (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS PiLRo, PiRLo, KoL, KoR, hEff) - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(v, 3, KoL, KoR, PiLRo, PiRLo, hEff, & + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, KoL, KoR, PiLRo, PiRLo, hEff, & (/1,2,3,3,3,3,3,3/), & ! kL (/1,1,1,2,3,3,3,3/), & ! kR (/0.,0.,0.,0.,0.,0.,.75,1./), & ! pL @@ -1244,13 +2091,13 @@ logical function neutral_diffusion_unit_tests(verbose) 'Right column with unstable mixed layer') ! Left column with unstable mixed layer - call find_neutral_surface_positions(3, & + call find_neutral_surface_positions_continuous(3, & (/0.,10.,20.,30./), (/10.,14.,12.,4./), (/0.,0.,0.,0./), & ! Left positions, T and S (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS (/0.,10.,20.,30./), (/14.,14.,10.,2./), (/0.,0.,0.,0./), & ! Right positions, T and S (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS PiLRo, PiRLo, KoL, KoR, hEff) - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(v, 3, KoL, KoR, PiLRo, PiRLo, hEff, & + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, KoL, KoR, PiLRo, PiRLo, hEff, & (/1,1,1,2,3,3,3,3/), & ! kL (/1,2,3,3,3,3,3,3/), & ! kR (/0.,0.,0.,0.,0.,0.25,1.,1./), & ! pL @@ -1259,13 +2106,13 @@ logical function neutral_diffusion_unit_tests(verbose) 'Left column with unstable mixed layer') ! Two unstable mixed layers - call find_neutral_surface_positions(3, & + call find_neutral_surface_positions_continuous(3, & (/0.,10.,20.,30./), (/8.,12.,10.,2./), (/0.,0.,0.,0./), & ! Left positions, T and S (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS (/0.,10.,20.,30./), (/10.,14.,12.,4./), (/0.,0.,0.,0./), & ! Right positions, T and S (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS PiLRo, PiRLo, KoL, KoR, hEff) - neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(v, 3, KoL, KoR, PiLRo, PiRLo, hEff, & + ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, KoL, KoR, PiLRo, PiRLo, hEff, & (/1,1,1,1,2,3,3,3/), & ! kL (/1,2,3,3,3,3,3,3/), & ! kR (/0.,0.,0.,0.,0.,0.,0.75,1./), & ! pL @@ -1273,9 +2120,181 @@ logical function neutral_diffusion_unit_tests(verbose) (/0.,0.,0.,0.,0.,6.,0./), & ! hEff 'Two unstable mixed layers') - if (.not. neutral_diffusion_unit_tests) write(*,*) 'Pass' + if (.not. ndiff_unit_tests_continuous) write(*,*) 'Pass' -end function neutral_diffusion_unit_tests +end function ndiff_unit_tests_continuous + +logical function ndiff_unit_tests_discontinuous(verbose) + logical, intent(in) :: verbose !< It true, write results to stdout + ! Local variables + integer, parameter :: nk = 3 + integer, parameter :: ns = nk*4 + real, dimension(nk) :: Sl, Sr, Tl, Tr, hl, hr + real, dimension(nk,2) :: TiL, SiL, TiR, SiR + real, dimension(nk+1) :: Pres_l, Pres_R + integer, dimension(ns) :: KoL, KoR + real, dimension(ns) :: PoL, PoR + real, dimension(ns-1) :: hEff, Flx + type(EOS_type), pointer :: EOS ! Structure for linear equation of state + type(remapping_CS), pointer :: remap_CS ! Remapping control structure (PLM) + real, dimension(nk,2) :: poly_T_l, poly_T_r, poly_S, poly_slope ! Linear reconstruction for T + real, dimension(nk,2) :: dRdT, dRdS + integer :: iMethod + integer :: k + logical :: v + + v = verbose + + ndiff_unit_tests_discontinuous = .false. ! Normally return false + write(*,*) '==== MOM_neutral_diffusion: ndiff_unit_tests_discontinuous =' + + ! Unit tests for find_neutral_surface_positions_discontinuous + ! Salinity is 0 for all these tests + Sl(:) = 0. ; Sr(:) = 0. ; poly_S(:,:) = 0. ; SiL(:,:) = 0. ; SiR(:,:) = 0. + dRdT(:,:) = -1. ; dRdS(:,:) = 0. + allocate(remap_CS) + call initialize_remapping( remap_CS, "PLM", boundary_extrapolation = .true. ) + + hL = (/10.,10.,10./) ; hR = (/10.,10.,10./) ; Pres_l(1) = 0. ; Pres_r(1) = 0. + do k = 1,nk ; Pres_l(k+1) = Pres_l(k) + hL(k) ; Pres_r(k+1) = Pres_r(k) + hR(k) ; enddo + ! Identical columns + Tl = (/20.,16.,12./) ; Tr = (/20.,16.,12./) + call build_reconstructions_1d( remap_CS, nk, hL, Tl, poly_T_l, TiL, poly_slope, iMethod ) + call build_reconstructions_1d( remap_CS, nk, hR, Tr, poly_T_r, TiR, poly_slope, iMethod ) + call find_neutral_surface_positions_discontinuous(nk, 1, Pres_l, TiL, SiL, dRdT, dRdS, & + Pres_r, TiR, SiR, dRdT, dRdS, PoL, PoR, KoL, KoR, hEff) + ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, KoL, KoR, PoL, PoR, hEff, & + (/1,1,1,1,2,2,2,2,3,3,3,3/), & ! KoL + (/1,1,1,1,2,2,2,2,3,3,3,3/), & ! KoR + (/0.,0.,1.,1.,0.,0.,1.,1.,0.,0.,1.,1./), & ! pL + (/0.,0.,1.,1.,0.,0.,1.,1.,0.,0.,1.,1./), & ! pR + (/0.,10.,0.,0.,0.,10.,0.,0.,0.,10.,0.,0./), & ! hEff + 'Identical columns') + Tl = (/20.,16.,12./) ; Tr = (/18.,14.,10./) + call build_reconstructions_1d( remap_CS, nk, hL, Tl, poly_T_l, TiL, poly_slope, iMethod ) + call build_reconstructions_1d( remap_CS, nk, hR, Tr, poly_T_r, TiR, poly_slope, iMethod ) + call find_neutral_surface_positions_discontinuous(nk, 1, Pres_l, TiL, SiL, dRdT, dRdS, & + Pres_r, TiR, SiR, dRdT, dRdS, PoL, PoR, KoL, KoR, hEff) + ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, KoL, KoR, PoL, PoR, hEff, & + (/1,1,1,2,2,2,2,3,3,3,3,3/), & ! KoL + (/1,1,1,1,1,2,2,2,2,3,3,3/), & ! KoR + (/0.0, 0.5, 1.0, 0.0, 0.5, 0.5, 1.0, 0.0, 0.5, 0.5, 1.0, 1.0/), & ! pL + (/0.0, 0.0, 0.5, 0.5, 1.0, 0.0, 0.5, 0.5, 1.0, 0.0, 0.5, 1.0/), & ! pR + (/0.0, 5.0, 0.0, 5.0, 0.0, 5.0, 0.0, 5.0, 0.0, 5.0, 0.0/), & ! hEff + 'Right column slightly cooler') + Tl = (/18.,14.,10./) ; Tr = (/20.,16.,12./) ; + call build_reconstructions_1d( remap_CS, nk, hL, Tl, poly_T_l, TiL, poly_slope, iMethod ) + call build_reconstructions_1d( remap_CS, nk, hR, Tr, poly_T_r, TiR, poly_slope, iMethod ) + call find_neutral_surface_positions_discontinuous(nk, 1, Pres_l, TiL, SiL, dRdT, dRdS, & + Pres_r, TiR, SiR, dRdT, dRdS, PoL, PoR, KoL, KoR, hEff) + ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, KoL, KoR, PoL, PoR, hEff, & + (/1,1,1,1,1,2,2,2,2,3,3,3/), & ! KoL + (/1,1,1,2,2,2,2,3,3,3,3,3/), & ! KoR + (/0.0, 0.0, 0.5, 0.5, 1.0, 0.0, 0.5, 0.5, 1.0, 0.0, 0.5, 1.0/), & ! pL + (/0.0, 0.5, 1.0, 0.0, 0.5, 0.5, 1.0, 0.0, 0.5, 0.5, 1.0, 1.0/), & ! pR + (/0.0, 5.0, 0.0, 5.0, 0.0, 5.0, 0.0, 5.0, 0.0, 5.0, 0.0/), & ! hEff + 'Left column slightly cooler') + Tl = (/20.,16.,12./) ; Tr = (/14.,10.,6./) + call build_reconstructions_1d( remap_CS, nk, hL, Tl, poly_T_l, TiL, poly_slope, iMethod ) + call build_reconstructions_1d( remap_CS, nk, hR, Tr, poly_T_r, TiR, poly_slope, iMethod ) + call find_neutral_surface_positions_discontinuous(nk, 1, Pres_l, TiL, SiL, dRdT, dRdS, & + Pres_r, TiR, SiR, dRdT, dRdS, PoL, PoR, KoL, KoR, hEff) + ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, KoL, KoR, PoL, PoR, hEff, & + (/1,1,2,2,2,3,3,3,3,3,3,3/), & ! KoL + (/1,1,1,1,1,1,1,2,2,2,3,3/), & ! KoR + (/0.0, 1.0, 0.0, 0.5, 1.0, 0.0, 0.5, 0.5, 1.0, 1.0, 1.0, 1.0/), & ! pL + (/0.0, 0.0, 0.0, 0.0, 0.5, 0.5, 1.0, 0.0, 0.5, 1.0, 0.0, 1.0/), & ! pR + (/0.0, 0.0, 0.0, 5.0, 0.0, 5.0, 0.0, 5.0, 0.0, 0.0, 0.0/), & ! hEff + 'Right column somewhat cooler') + Tl = (/20.,16.,12./) ; Tr = (/8.,6.,4./) + call build_reconstructions_1d( remap_CS, nk, hL, Tl, poly_T_l, TiL, poly_slope, iMethod ) + call build_reconstructions_1d( remap_CS, nk, hR, Tr, poly_T_r, TiR, poly_slope, iMethod ) + call find_neutral_surface_positions_discontinuous(nk, 1, Pres_l, TiL, SiL, dRdT, dRdS, & + Pres_r, TiR, SiR, dRdT, dRdS, PoL, PoR, KoL, KoR, hEff) + ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, KoL, KoR, PoL, PoR, hEff, & + (/1,1,2,2,3,3,3,3,3,3,3,3/), & ! KoL + (/1,1,1,1,1,1,1,1,2,2,3,3/), & ! KoR + (/0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0/), & ! pL + (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0/), & ! pR + (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/), & ! hEff + 'Right column much cooler') + Tl = (/14.,14.,10./) ; Tr = (/14.,14.,10./) + call build_reconstructions_1d( remap_CS, nk, hL, Tl, poly_T_l, TiL, poly_slope, iMethod ) + call build_reconstructions_1d( remap_CS, nk, hR, Tr, poly_T_r, TiR, poly_slope, iMethod ) + call find_neutral_surface_positions_discontinuous(nk, 1, Pres_l, TiL, SiL, dRdT, dRdS, & + Pres_r, TiR, SiR, dRdT, dRdS, PoL, PoR, KoL, KoR, hEff) + ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, KoL, KoR, PoL, PoR, hEff, & + (/1,1,1,1,2,2,2,2,3,3,3,3/), & ! KoL + (/1,1,1,1,2,2,2,2,3,3,3,3/), & ! KoR + (/0.,0.,1.,1.,0.,0.,1.,1.,0.,0.,1.,1./), & ! pL + (/0.,0.,1.,1.,0.,0.,1.,1.,0.,0.,1.,1./), & ! pR + (/0.,10.,0.,0.,0.,10.,0.,0.,0.,10.,0.,0./), & ! hEff + 'Identical columns with mixed layer') + Tl = (/20.,16.,12./) ; Tr = (/14.,14.,10./) + call build_reconstructions_1d( remap_CS, nk, hL, Tl, poly_T_l, TiL, poly_slope, iMethod ) + call build_reconstructions_1d( remap_CS, nk, hR, Tr, poly_T_r, TiR, poly_slope, iMethod ) + call find_neutral_surface_positions_discontinuous(nk, 1, Pres_l, TiL, SiL, dRdT, dRdS, & + Pres_r, TiR, SiR, dRdT, dRdS, PoL, PoR, KoL, KoR, hEff) + ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, KoL, KoR, PoL, PoR, hEff, & + (/1,1,2,2,2,3,3,3,3,3,3,3/), & ! KoL + (/1,1,1,1,1,1,2,2,2,3,3,3/), & ! KoR + (/0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 1.0/), & ! pL + (/0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.5, 1.0/), & ! pR + (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 5.0, 0.0/), & ! hEff + 'Right column with mixed layer') + Tl = (/14.,14.,6./) ; Tr = (/12.,16.,8./) + call build_reconstructions_1d( remap_CS, nk, hL, Tl, poly_T_l, TiL, poly_slope, iMethod ) + call build_reconstructions_1d( remap_CS, nk, hR, Tr, poly_T_r, TiR, poly_slope, iMethod ) + call find_neutral_surface_positions_discontinuous(nk, 1, Pres_l, TiL, SiL, dRdT, dRdS, & + Pres_r, TiR, SiR, dRdT, dRdS, PoL, PoR, KoL, KoR, hEff) + ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, KoL, KoR, PoL, PoR, hEff, & + (/1,1,2,2,3,3,3,3,3,3,3,3/), & ! KoL + (/1,1,1,1,1,1,2,2,3,3,3,3/), & ! KoR + (/0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, .75, 1.0/), & ! pL + (/0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, .25, 1.0, 1.0/), & ! pR + (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 7.5, 0.0/), & ! hEff + 'Left mixed layer, right unstable mixed layer') + Til(:,1) = (/8.,12.,10./) ; Til(:,2) = (/12.,10.,2./) + Tir(:,1) = (/10.,14.,12./) ; TiR(:,2) = (/14.,12.,4./) + call find_neutral_surface_positions_discontinuous(nk, 1, Pres_l, TiL, SiL, dRdT, dRdS, & + Pres_r, TiR, SiR, dRdT, dRdS, PoL, PoR, KoL, KoR, hEff) + ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, KoL, KoR, PoL, PoR, hEff, & + (/1,1,1,1,1,1,1,2,2,3,3,3/), & ! KoL + (/1,1,2,2,3,3,3,3,3,3,3,3/), & ! KoR + (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, .75, 1.0/), & ! pL + (/0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, .25, .25, 1.0, 1.0/), & ! pR + (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 7.5, 0.0/), & ! hEff + 'Two unstable mixed layers') + deallocate(remap_CS) + + allocate(EOS) + call EOS_manual_init(EOS, form_of_EOS = EOS_LINEAR, dRho_dT = -1., dRho_dS = 2.) + ! Unit tests for refine_nondim_position + ! Tests using Newton's method + ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5,refine_nondim_position( & + 100, 0., 20., 35., -1., 2., 0., 1., 1, (/21., -2./), (/35., 0./), EOS, 0., -1., 1., 0.), & + "Temperature stratified (Newton) ")) + ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5,refine_nondim_position( & + 100, 0., 20., 35., -1., 2., 0., 1., 1, (/20., 0./), (/34., 2./), EOS, 0., -2., 2., 0.), & + "Salinity stratified (Newton) ")) + ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5,refine_nondim_position( & + 100, 0., 20., 35., -1., 2., 0., 1., 1, (/21., -2./), (/34., 2./), EOS, 0., -1., 1., 0.), & + "Temp/Salt stratified (Newton) ")) + ! Tests using Brent's method + ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5,refine_nondim_position( & + 100, 0., 20., 35., -1., 2., 0., 1., 1, (/21., -2./), (/35., 0./), EOS, 0., -1., 1., 0.,force_brent = .true.), & + "Temperature stratified (Brent) ")) + ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5,refine_nondim_position( & + 100, 0., 20., 35., -1., 2., 0., 1., 1, (/20., 0./), (/34., 2./), EOS, 0., -2., 2., 0.,force_brent = .true.), & + "Salinity stratified (Brent) ")) + ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5,refine_nondim_position( & + 100, 0., 20., 35., -1., 2., 0., 1., 1, (/21., -2./), (/34., 2./), EOS, 0., -1., 1., 0.,force_brent = .true.), & + "Temp/Salt stratified (Brent) ")) + deallocate(EOS) + + if (.not. ndiff_unit_tests_discontinuous) write(*,*) 'Pass' + +end function ndiff_unit_tests_discontinuous !> Returns true if a test of fv_diff() fails, and conditionally writes results to stream logical function test_fv_diff(verbose, hkm1, hk, hkp1, Skm1, Sk, Skp1, Ptrue, title) @@ -1438,19 +2457,19 @@ logical function test_data1di(verbose, nk, Po, Ptrue, title) end function test_data1di !> Returns true if output of find_neutral_surface_positions() does not match correct values, and conditionally writes results to stream -logical function test_nsp(verbose, nk, KoL, KoR, pL, pR, hEff, KoL0, KoR0, pL0, pR0, hEff0, title) +logical function test_nsp(verbose, ns, KoL, KoR, pL, pR, hEff, KoL0, KoR0, pL0, pR0, hEff0, title) logical, intent(in) :: verbose !< If true, write results to stdout - integer, intent(in) :: nk !< Number of layers - integer, dimension(2*nk+2), intent(in) :: KoL !< Index of first left interface above neutral surface - integer, dimension(2*nk+2), intent(in) :: KoR !< Index of first right interface above neutral surface - real, dimension(2*nk+2), intent(in) :: pL !< Fractional position of neutral surface within layer KoL of left column - real, dimension(2*nk+2), intent(in) :: pR !< Fractional position of neutral surface within layer KoR of right column - real, dimension(2*nk+1), intent(in) :: hEff !< Effective thickness between two neutral surfaces (Pa) - integer, dimension(2*nk+2), intent(in) :: KoL0 !< Correct value for KoL - integer, dimension(2*nk+2), intent(in) :: KoR0 !< Correct value for KoR - real, dimension(2*nk+2), intent(in) :: pL0 !< Correct value for pL - real, dimension(2*nk+2), intent(in) :: pR0 !< Correct value for pR - real, dimension(2*nk+1), intent(in) :: hEff0 !< Correct value for hEff + integer, intent(in) :: ns !< Number of surfaces + integer, dimension(ns), intent(in) :: KoL !< Index of first left interface above neutral surface + integer, dimension(ns), intent(in) :: KoR !< Index of first right interface above neutral surface + real, dimension(ns), intent(in) :: pL !< Fractional position of neutral surface within layer KoL of left column + real, dimension(ns), intent(in) :: pR !< Fractional position of neutral surface within layer KoR of right column + real, dimension(ns-1), intent(in) :: hEff !< Effective thickness between two neutral surfaces (Pa) + integer, dimension(ns), intent(in) :: KoL0 !< Correct value for KoL + integer, dimension(ns), intent(in) :: KoR0 !< Correct value for KoR + real, dimension(ns), intent(in) :: pL0 !< Correct value for pL + real, dimension(ns), intent(in) :: pR0 !< Correct value for pR + real, dimension(ns-1), intent(in) :: hEff0 !< Correct value for hEff character(len=*), intent(in) :: title !< Title for messages ! Local variables @@ -1458,9 +2477,9 @@ logical function test_nsp(verbose, nk, KoL, KoR, pL, pR, hEff, KoL0, KoR0, pL0, logical :: this_row_failed test_nsp = .false. - do k = 1,2*nk+2 + do k = 1,ns test_nsp = test_nsp .or. compare_nsp_row(KoL(k), KoR(k), pL(k), pR(k), KoL0(k), KoR0(k), pL0(k), pR0(k)) - if (k < 2*nk+2) then + if (k < ns) then if (hEff(k) /= hEff0(k)) test_nsp = .true. endif enddo @@ -1469,7 +2488,7 @@ logical function test_nsp(verbose, nk, KoL, KoR, pL, pR, hEff, KoL0, KoR0, pL0, stdunit = 6 if (test_nsp) stdunit = 0 ! In case of wrong results, write to error stream write(stdunit,'(a)') title - do k = 1,2*nk+2 + do k = 1,ns this_row_failed = compare_nsp_row(KoL(k), KoR(k), pL(k), pR(k), KoL0(k), KoR0(k), pL0(k), pR0(k)) if (this_row_failed) then write(stdunit,10) k,KoL(k),pL(k),KoR(k),pR(k),' <-- WRONG!' @@ -1477,7 +2496,7 @@ logical function test_nsp(verbose, nk, KoL, KoR, pL, pR, hEff, KoL0, KoR0, pL0, else write(stdunit,10) k,KoL(k),pL(k),KoR(k),pR(k) endif - if (k < 2*nk+2) then + if (k < ns) then if (hEff(k) /= hEff0(k)) then write(stdunit,'(i3,8x,"layer hEff =",2(f20.16,a))') k,hEff(k)," .neq. ",hEff0(k),' <-- WRONG!' else @@ -1486,6 +2505,7 @@ logical function test_nsp(verbose, nk, KoL, KoR, pL, pR, hEff, KoL0, KoR0, pL0, endif enddo endif + if (test_nsp) call MOM_error(FATAL,"test_nsp failed") 10 format("ks=",i3," kL=",i3," pL=",f20.16," kR=",i3," pR=",f20.16,a) end function test_nsp @@ -1508,6 +2528,20 @@ logical function compare_nsp_row(KoL, KoR, pL, pR, KoL0, KoR0, pL0, pR0) if (pR /= pR0) compare_nsp_row = .true. end function compare_nsp_row +!> Compares output position from refine_nondim_position with an expected value +logical function test_rnp(expected_pos, test_pos, title) + real, intent(in) :: expected_pos + real, intent(in) :: test_pos + character(len=*), intent(in) :: title + ! Local variables + integer :: stdunit = 6 ! Output to standard error + test_rnp = expected_pos /= test_pos + if (test_rnp) then + write(stdunit,'(A, f20.16, " .neq. ", f20.16, " <-- WRONG")') title, expected_pos, test_pos + else + write(stdunit,'(A, f20.16, " .eq. ", f20.16)') title, expected_pos, test_pos + endif +end function test_rnp !> Deallocates neutral_diffusion control structure subroutine neutral_diffusion_end(CS) type(neutral_diffusion_CS), pointer :: CS