MOM6
Data Types | Functions/Subroutines | Variables
mom_neutral_diffusion Module Reference

Detailed Description

A column-wise toolbox for implementing neutral diffusion.

Data Types

type  neutral_diffusion_cs
 The control structure for the MOM_neutral_diffusion module. More...
 

Functions/Subroutines

logical function, public neutral_diffusion_init (Time, G, param_file, diag, EOS, CS)
 Read parameters and allocate control structure for neutral_diffusion module. More...
 
subroutine, public neutral_diffusion_calc_coeffs (G, GV, h, T, S, CS)
 Calculate remapping factors for u/v columns used to map adjoining columns to a shared coordinate space. More...
 
subroutine, public neutral_diffusion (G, GV, h, Coef_x, Coef_y, dt, Reg, CS)
 Update tracer concentration due to neutral diffusion; layer thickness unchanged by this update. More...
 
subroutine interface_scalar (nk, h, S, Si, i_method, h_neglect)
 Returns interface scalar, Si, for a column of layer values, S. More...
 
real function ppm_edge (hkm1, hk, hkp1, hkp2, Ak, Akp1, Pk, Pkp1, h_neglect)
 Returns the PPM quasi-fourth order edge value at k+1/2 following equation 1.6 in Colella & Woodward, 1984: JCP 54, 174-201. More...
 
real function ppm_ave (xL, xR, aL, aR, aMean)
 Returns the average of a PPM reconstruction between two fractional positions. More...
 
real function signum (a, x)
 A true signum function that returns either -abs(a), when x<0; or abs(a) when x>0; or 0 when x=0. More...
 
subroutine plm_diff (nk, h, S, c_method, b_method, diff)
 Returns PLM slopes for a column where the slopes are the difference in value across each cell. The limiting follows equation 1.8 in Colella & Woodward, 1984: JCP 54, 174-201. More...
 
real function fv_diff (hkm1, hk, hkp1, Skm1, Sk, Skp1)
 Returns the cell-centered second-order finite volume (unlimited PLM) slope using three consecutive cell widths and average values. Slope is returned as a difference across the central cell (i.e. units of scalar S). Discretization follows equation 1.7 in Colella & Woodward, 1984: JCP 54, 174-201. More...
 
real function fvlsq_slope (hkm1, hk, hkp1, Skm1, Sk, Skp1)
 Returns the cell-centered second-order weighted least squares slope using three consecutive cell widths and average values. Slope is returned as a gradient (i.e. units of scalar S over width units). More...
 
subroutine find_neutral_surface_positions_continuous (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. More...
 
subroutine find_neutral_surface_positions_discontinuous (CS, nk, ns, Pres_l, hcol_l, Tl, Sl, dRdT_l, dRdS_l, stable_l, Pres_r, hcol_r, Tr, Sr, dRdT_r, dRdS_r, stable_r, PoL, PoR, KoL, KoR, hEff, ppoly_T_l, ppoly_S_l, ppoly_T_r, ppoly_S_r)
 Higher order version of find_neutral_surface_positions. Returns positions within left/right columns of combined interfaces using intracell reconstructions of T/S. More...
 
real function absolute_position (n, ns, Pint, Karr, NParr, k_surface)
 Converts non-dimensional position within a layer to absolute position (for debugging) More...
 
real function, dimension(ns) absolute_positions (n, ns, Pint, Karr, NParr)
 Converts non-dimensional positions within layers to absolute positions (for debugging) More...
 
subroutine neutral_surface_flux (nk, nsurf, deg, hl, hr, Tl, Tr, PiL, PiR, KoL, KoR, hEff, Flx, continuous, h_neglect, remap_CS, h_neglect_edge)
 Returns a single column of neutral diffusion fluxes of a tracer. More...
 
subroutine neutral_surface_t_eval (nk, ns, k_sub, Ks, Ps, T_mean, T_int, deg, iMethod, T_poly, T_top, T_bot, T_sub, T_top_int, T_bot_int, T_layer)
 Evaluate various parts of the reconstructions to calculate gradient-based flux limter. More...
 
subroutine ppm_left_right_edge_values (nk, Tl, Ti, aL, aR)
 Discontinuous PPM reconstructions of the left/right edge values within a cell. More...
 
logical function, public neutral_diffusion_unit_tests (verbose)
 Returns true if unit tests of neutral_diffusion functions fail. Otherwise returns false. More...
 
logical function ndiff_unit_tests_continuous (verbose)
 Returns true if unit tests of neutral_diffusion functions fail. Otherwise returns false. More...
 
logical function ndiff_unit_tests_discontinuous (verbose)
 
logical function test_fv_diff (verbose, hkm1, hk, hkp1, Skm1, Sk, Skp1, Ptrue, title)
 Returns true if a test of fv_diff() fails, and conditionally writes results to stream. More...
 
logical function test_fvlsq_slope (verbose, hkm1, hk, hkp1, Skm1, Sk, Skp1, Ptrue, title)
 Returns true if a test of fvlsq_slope() fails, and conditionally writes results to stream. More...
 
logical function test_ifndp (verbose, rhoNeg, Pneg, rhoPos, Ppos, Ptrue, title)
 Returns true if a test of interpolate_for_nondim_position() fails, and conditionally writes results to stream. More...
 
logical function test_data1d (verbose, nk, Po, Ptrue, title)
 Returns true if comparison of Po and Ptrue fails, and conditionally writes results to stream. More...
 
logical function test_data1di (verbose, nk, Po, Ptrue, title)
 Returns true if comparison of Po and Ptrue fails, and conditionally writes results to stream. More...
 
logical function test_nsp (verbose, ns, KoL, KoR, pL, pR, hEff, KoL0, KoR0, pL0, pR0, hEff0, title)
 Returns true if output of find_neutral_surface_positions() does not match correct values, and conditionally writes results to stream. More...
 
logical function compare_nsp_row (KoL, KoR, pL, pR, KoL0, KoR0, pL0, pR0)
 Compares a single row, k, of output from find_neutral_surface_positions() More...
 
logical function test_rnp (expected_pos, test_pos, title)
 Compares output position from refine_nondim_position with an expected value. More...
 
subroutine, public neutral_diffusion_end (CS)
 Deallocates neutral_diffusion control structure. More...
 

Variables

character(len=40) mdl = "MOM_neutral_diffusion"
 module name
 

Function/Subroutine Documentation

◆ absolute_position()

real function mom_neutral_diffusion::absolute_position ( integer, intent(in)  n,
integer, intent(in)  ns,
real, dimension(n+1), intent(in)  Pint,
integer, dimension(ns), intent(in)  Karr,
real, dimension(ns), intent(in)  NParr,
integer, intent(in)  k_surface 
)
private

Converts non-dimensional position within a layer to absolute position (for debugging)

Parameters
[in]nNumber of levels
[in]nsNumber of neutral surfaces
[in]pintPosition of interfaces [Pa]
[in]karrIndex of interface above position
[in]nparrNon-dimensional position within layer Karr(:)
[in]k_surfacek-interface to query

Definition at line 1266 of file MOM_neutral_diffusion.F90.

1266  integer, intent(in) :: n !< Number of levels
1267  integer, intent(in) :: ns !< Number of neutral surfaces
1268  real, intent(in) :: Pint(n+1) !< Position of interfaces [Pa]
1269  integer, intent(in) :: Karr(ns) !< Index of interface above position
1270  real, intent(in) :: NParr(ns) !< Non-dimensional position within layer Karr(:)
1271  integer, intent(in) :: k_surface !< k-interface to query
1272  ! Local variables
1273  integer :: k
1274 
1275  k = karr(k_surface)
1276  if (k>n) stop 'absolute_position: k>nk is out of bounds!'
1277  absolute_position = pint(k) + nparr(k_surface) * ( pint(k+1) - pint(k) )
1278 

◆ absolute_positions()

real function, dimension(ns) mom_neutral_diffusion::absolute_positions ( integer, intent(in)  n,
integer, intent(in)  ns,
real, dimension(n+1), intent(in)  Pint,
integer, dimension(ns), intent(in)  Karr,
real, dimension(ns), intent(in)  NParr 
)
private

Converts non-dimensional positions within layers to absolute positions (for debugging)

Parameters
[in]nNumber of levels
[in]nsNumber of neutral surfaces
[in]pintPosition of interface [Pa]
[in]karrIndexes of interfaces about positions
[in]nparrNon-dimensional positions within layers Karr(:)

Definition at line 1283 of file MOM_neutral_diffusion.F90.

1283  integer, intent(in) :: n !< Number of levels
1284  integer, intent(in) :: ns !< Number of neutral surfaces
1285  real, intent(in) :: Pint(n+1) !< Position of interface [Pa]
1286  integer, intent(in) :: Karr(ns) !< Indexes of interfaces about positions
1287  real, intent(in) :: NParr(ns) !< Non-dimensional positions within layers Karr(:)
1288 
1289  real, dimension(ns) :: absolute_positions ! Absolute positions [Pa]
1290 
1291  ! Local variables
1292  integer :: k_surface, k
1293 
1294  do k_surface = 1, ns
1295  absolute_positions(k_surface) = absolute_position(n,ns,pint,karr,nparr,k_surface)
1296  enddo
1297 

◆ compare_nsp_row()

logical function mom_neutral_diffusion::compare_nsp_row ( integer, intent(in)  KoL,
integer, intent(in)  KoR,
real, intent(in)  pL,
real, intent(in)  pR,
integer, intent(in)  KoL0,
integer, intent(in)  KoR0,
real, intent(in)  pL0,
real, intent(in)  pR0 
)
private

Compares a single row, k, of output from find_neutral_surface_positions()

Parameters
[in]kolIndex of first left interface above neutral surface
[in]korIndex of first right interface above neutral surface
[in]plFractional position of neutral surface within layer KoL of left column
[in]prFractional position of neutral surface within layer KoR of right column
[in]kol0Correct value for KoL
[in]kor0Correct value for KoR
[in]pl0Correct value for pL
[in]pr0Correct value for pR

Definition at line 2208 of file MOM_neutral_diffusion.F90.

2208  integer, intent(in) :: KoL !< Index of first left interface above neutral surface
2209  integer, intent(in) :: KoR !< Index of first right interface above neutral surface
2210  real, intent(in) :: pL !< Fractional position of neutral surface within layer KoL of left column
2211  real, intent(in) :: pR !< Fractional position of neutral surface within layer KoR of right column
2212  integer, intent(in) :: KoL0 !< Correct value for KoL
2213  integer, intent(in) :: KoR0 !< Correct value for KoR
2214  real, intent(in) :: pL0 !< Correct value for pL
2215  real, intent(in) :: pR0 !< Correct value for pR
2216 
2217  compare_nsp_row = .false.
2218  if (kol /= kol0) compare_nsp_row = .true.
2219  if (kor /= kor0) compare_nsp_row = .true.
2220  if (pl /= pl0) compare_nsp_row = .true.
2221  if (pr /= pr0) compare_nsp_row = .true.

◆ find_neutral_surface_positions_continuous()

subroutine mom_neutral_diffusion::find_neutral_surface_positions_continuous ( integer, intent(in)  nk,
real, dimension(nk+1), intent(in)  Pl,
real, dimension(nk+1), intent(in)  Tl,
real, dimension(nk+1), intent(in)  Sl,
real, dimension(nk+1), intent(in)  dRdTl,
real, dimension(nk+1), intent(in)  dRdSl,
real, dimension(nk+1), intent(in)  Pr,
real, dimension(nk+1), intent(in)  Tr,
real, dimension(nk+1), intent(in)  Sr,
real, dimension(nk+1), intent(in)  dRdTr,
real, dimension(nk+1), intent(in)  dRdSr,
real, dimension(2*nk+2), intent(inout)  PoL,
real, dimension(2*nk+2), intent(inout)  PoR,
integer, dimension(2*nk+2), intent(inout)  KoL,
integer, dimension(2*nk+2), intent(inout)  KoR,
real, dimension(2*nk+1), intent(inout)  hEff 
)
private

Returns positions within left/right columns of combined interfaces using continuous reconstructions of T/S.

Parameters
[in]nkNumber of levels
[in]plLeft-column interface pressure [Pa]
[in]tlLeft-column interface potential temperature [degC]
[in]slLeft-column interface salinity [ppt]
[in]drdtlLeft-column dRho/dT [kg m-3 degC-1]
[in]drdslLeft-column dRho/dS [kg m-3 ppt-1]
[in]prRight-column interface pressure [Pa]
[in]trRight-column interface potential temperature [degC]
[in]srRight-column interface salinity [ppt]
[in]drdtrLeft-column dRho/dT [kg m-3 degC-1]
[in]drdsrLeft-column dRho/dS [kg m-3 ppt-1]
[in,out]polFractional position of neutral surface within layer KoL of left column
[in,out]porFractional position of neutral surface within layer KoR of right column
[in,out]kolIndex of first left interface above neutral surface
[in,out]korIndex of first right interface above neutral surface
[in,out]heffEffective thickness between two neutral surfaces [Pa]

Definition at line 811 of file MOM_neutral_diffusion.F90.

811  integer, intent(in) :: nk !< Number of levels
812  real, dimension(nk+1), intent(in) :: Pl !< Left-column interface pressure [Pa]
813  real, dimension(nk+1), intent(in) :: Tl !< Left-column interface potential temperature [degC]
814  real, dimension(nk+1), intent(in) :: Sl !< Left-column interface salinity [ppt]
815  real, dimension(nk+1), intent(in) :: dRdTl !< Left-column dRho/dT [kg m-3 degC-1]
816  real, dimension(nk+1), intent(in) :: dRdSl !< Left-column dRho/dS [kg m-3 ppt-1]
817  real, dimension(nk+1), intent(in) :: Pr !< Right-column interface pressure [Pa]
818  real, dimension(nk+1), intent(in) :: Tr !< Right-column interface potential temperature [degC]
819  real, dimension(nk+1), intent(in) :: Sr !< Right-column interface salinity [ppt]
820  real, dimension(nk+1), intent(in) :: dRdTr !< Left-column dRho/dT [kg m-3 degC-1]
821  real, dimension(nk+1), intent(in) :: dRdSr !< Left-column dRho/dS [kg m-3 ppt-1]
822  real, dimension(2*nk+2), intent(inout) :: PoL !< Fractional position of neutral surface within
823  !! layer KoL of left column
824  real, dimension(2*nk+2), intent(inout) :: PoR !< Fractional position of neutral surface within
825  !! layer KoR of right column
826  integer, dimension(2*nk+2), intent(inout) :: KoL !< Index of first left interface above neutral surface
827  integer, dimension(2*nk+2), intent(inout) :: KoR !< Index of first right interface above neutral surface
828  real, dimension(2*nk+1), intent(inout) :: hEff !< Effective thickness between two neutral surfaces [Pa]
829 
830  ! Local variables
831  integer :: ns ! Number of neutral surfaces
832  integer :: k_surface ! Index of neutral surface
833  integer :: kl ! Index of left interface
834  integer :: kr ! Index of right interface
835  real :: dRdT, dRdS ! dRho/dT and dRho/dS for the neutral surface
836  logical :: searching_left_column ! True if searching for the position of a right interface in the left column
837  logical :: searching_right_column ! True if searching for the position of a left interface in the right column
838  logical :: reached_bottom ! True if one of the bottom-most interfaces has been used as the target
839  integer :: krm1, klm1
840  real :: dRho, dRhoTop, dRhoBot, hL, hR
841  integer :: lastK_left, lastK_right
842  real :: lastP_left, lastP_right
843 
844  ns = 2*nk+2
845  ! Initialize variables for the search
846  kr = 1 ; lastk_right = 1 ; lastp_right = 0.
847  kl = 1 ; lastk_left = 1 ; lastp_left = 0.
848  reached_bottom = .false.
849 
850  ! Loop over each neutral surface, working from top to bottom
851  neutral_surfaces: do k_surface = 1, ns
852  klm1 = max(kl-1, 1)
853  if (klm1>nk) stop 'find_neutral_surface_positions(): klm1 went out of bounds!'
854  krm1 = max(kr-1, 1)
855  if (krm1>nk) stop 'find_neutral_surface_positions(): krm1 went out of bounds!'
856 
857  ! Potential density difference, rho(kr) - rho(kl)
858  drho = 0.5 * ( ( drdtr(kr) + drdtl(kl) ) * ( tr(kr) - tl(kl) ) &
859  + ( drdsr(kr) + drdsl(kl) ) * ( sr(kr) - sl(kl) ) )
860  ! Which column has the lighter surface for the current indexes, kr and kl
861  if (.not. reached_bottom) then
862  if (drho < 0.) then
863  searching_left_column = .true.
864  searching_right_column = .false.
865  elseif (drho > 0.) then
866  searching_right_column = .true.
867  searching_left_column = .false.
868  else ! dRho == 0.
869  if (kl + kr == 2) then ! Still at surface
870  searching_left_column = .true.
871  searching_right_column = .false.
872  else ! Not the surface so we simply change direction
873  searching_left_column = .not. searching_left_column
874  searching_right_column = .not. searching_right_column
875  endif
876  endif
877  endif
878 
879  if (searching_left_column) then
880  ! Interpolate for the neutral surface position within the left column, layer klm1
881  ! Potential density difference, rho(kl-1) - rho(kr) (should be negative)
882  drhotop = 0.5 * ( ( drdtl(klm1) + drdtr(kr) ) * ( tl(klm1) - tr(kr) ) &
883  + ( drdsl(klm1) + drdsr(kr) ) * ( sl(klm1) - sr(kr) ) )
884  ! Potential density difference, rho(kl) - rho(kr) (will be positive)
885  drhobot = 0.5 * ( ( drdtl(klm1+1) + drdtr(kr) ) * ( tl(klm1+1) - tr(kr) ) &
886  + ( drdsl(klm1+1) + drdsr(kr) ) * ( sl(klm1+1) - sr(kr) ) )
887 
888  ! Because we are looking left, the right surface, kr, is lighter than klm1+1 and should be denser than klm1
889  ! unless we are still at the top of the left column (kl=1)
890  if (drhotop > 0. .or. kr+kl==2) then
891  pol(k_surface) = 0. ! The right surface is lighter than anything in layer klm1
892  elseif (drhotop >= drhobot) then ! Left layer is unstratified
893  pol(k_surface) = 1.
894  else
895  ! Linearly interpolate for the position between Pl(kl-1) and Pl(kl) where the density difference
896  ! between right and left is zero.
897  pol(k_surface) = interpolate_for_nondim_position( drhotop, pl(klm1), drhobot, pl(klm1+1) )
898  endif
899  if (pol(k_surface)>=1. .and. klm1<nk) then ! >= is really ==, when PoL==1 we point to the bottom of the cell
900  klm1 = klm1 + 1
901  pol(k_surface) = pol(k_surface) - 1.
902  endif
903  if (real(klm1-lastk_left)+(pol(k_surface)-lastp_left)<0.) then
904  pol(k_surface) = lastp_left
905  klm1 = lastk_left
906  endif
907  kol(k_surface) = klm1
908  if (kr <= nk) then
909  por(k_surface) = 0.
910  kor(k_surface) = kr
911  else
912  por(k_surface) = 1.
913  kor(k_surface) = nk
914  endif
915  if (kr <= nk) then
916  kr = kr + 1
917  else
918  reached_bottom = .true.
919  searching_right_column = .true.
920  searching_left_column = .false.
921  endif
922  elseif (searching_right_column) then
923  ! Interpolate for the neutral surface position within the right column, layer krm1
924  ! Potential density difference, rho(kr-1) - rho(kl) (should be negative)
925  drhotop = 0.5 * ( ( drdtr(krm1) + drdtl(kl) ) * ( tr(krm1) - tl(kl) ) &
926  + ( drdsr(krm1) + drdsl(kl) ) * ( sr(krm1) - sl(kl) ) )
927  ! Potential density difference, rho(kr) - rho(kl) (will be positive)
928  drhobot = 0.5 * ( ( drdtr(krm1+1) + drdtl(kl) ) * ( tr(krm1+1) - tl(kl) ) &
929  + ( drdsr(krm1+1) + drdsl(kl) ) * ( sr(krm1+1) - sl(kl) ) )
930 
931  ! Because we are looking right, the left surface, kl, is lighter than krm1+1 and should be denser than krm1
932  ! unless we are still at the top of the right column (kr=1)
933  if (drhotop >= 0. .or. kr+kl==2) then
934  por(k_surface) = 0. ! The left surface is lighter than anything in layer krm1
935  elseif (drhotop >= drhobot) then ! Right layer is unstratified
936  por(k_surface) = 1.
937  else
938  ! Linearly interpolate for the position between Pr(kr-1) and Pr(kr) where the density difference
939  ! between right and left is zero.
940  por(k_surface) = interpolate_for_nondim_position( drhotop, pr(krm1), drhobot, pr(krm1+1) )
941  endif
942  if (por(k_surface)>=1. .and. krm1<nk) then ! >= is really ==, when PoR==1 we point to the bottom of the cell
943  krm1 = krm1 + 1
944  por(k_surface) = por(k_surface) - 1.
945  endif
946  if (real(krm1-lastk_right)+(por(k_surface)-lastp_right)<0.) then
947  por(k_surface) = lastp_right
948  krm1 = lastk_right
949  endif
950  kor(k_surface) = krm1
951  if (kl <= nk) then
952  pol(k_surface) = 0.
953  kol(k_surface) = kl
954  else
955  pol(k_surface) = 1.
956  kol(k_surface) = nk
957  endif
958  if (kl <= nk) then
959  kl = kl + 1
960  else
961  reached_bottom = .true.
962  searching_right_column = .false.
963  searching_left_column = .true.
964  endif
965  else
966  stop 'Else what?'
967  endif
968 
969  lastk_left = kol(k_surface) ; lastp_left = pol(k_surface)
970  lastk_right = kor(k_surface) ; lastp_right = por(k_surface)
971 
972  ! Effective thickness
973  ! NOTE: This would be better expressed in terms of the layers thicknesses rather
974  ! than as differences of position - AJA
975  if (k_surface>1) then
976  hl = absolute_position(nk,ns,pl,kol,pol,k_surface) - absolute_position(nk,ns,pl,kol,pol,k_surface-1)
977  hr = absolute_position(nk,ns,pr,kor,por,k_surface) - absolute_position(nk,ns,pr,kor,por,k_surface-1)
978  if ( hl + hr > 0.) then
979  heff(k_surface-1) = 2. * hl * hr / ( hl + hr ) ! Harmonic mean of layer thicknesses
980  else
981  heff(k_surface-1) = 0.
982  endif
983  endif
984 
985  enddo neutral_surfaces
986 

◆ find_neutral_surface_positions_discontinuous()

subroutine mom_neutral_diffusion::find_neutral_surface_positions_discontinuous ( type(neutral_diffusion_cs), intent(inout)  CS,
integer, intent(in)  nk,
integer, intent(in)  ns,
real, dimension(nk+1), intent(in)  Pres_l,
real, dimension(nk), intent(in)  hcol_l,
real, dimension(nk,2), intent(in)  Tl,
real, dimension(nk,2), intent(in)  Sl,
real, dimension(nk,2), intent(in)  dRdT_l,
real, dimension(nk,2), intent(in)  dRdS_l,
logical, dimension(nk), intent(in)  stable_l,
real, dimension(nk+1), intent(in)  Pres_r,
real, dimension(nk), intent(in)  hcol_r,
real, dimension(nk,2), intent(in)  Tr,
real, dimension(nk,2), intent(in)  Sr,
real, dimension(nk,2), intent(in)  dRdT_r,
real, dimension(nk,2), intent(in)  dRdS_r,
logical, dimension(nk), intent(in)  stable_r,
real, dimension(4*nk), intent(inout)  PoL,
real, dimension(4*nk), intent(inout)  PoR,
integer, dimension(4*nk), intent(inout)  KoL,
integer, dimension(4*nk), intent(inout)  KoR,
real, dimension(4*nk-1), intent(inout)  hEff,
real, dimension(nk,cs%deg+1), intent(in), optional  ppoly_T_l,
real, dimension(nk,cs%deg+1), intent(in), optional  ppoly_S_l,
real, dimension(nk,cs%deg+1), intent(in), optional  ppoly_T_r,
real, dimension(nk,cs%deg+1), intent(in), optional  ppoly_S_r 
)
private

Higher order version of find_neutral_surface_positions. Returns positions within left/right columns of combined interfaces using intracell reconstructions of T/S.

Parameters
[in,out]csNeutral diffusion control structure
[in]nkNumber of levels
[in]nsNumber of neutral surfaces
[in]pres_lLeft-column interface pressure [Pa]
[in]hcol_lLeft-column layer thicknesses
[in]tlLeft-column top interface potential temperature [degC]
[in]slLeft-column top interface salinity [ppt]
[in]drdt_lLeft-column, top interface dRho/dT [kg m-3 degC-1]
[in]drds_lLeft-column, top interface dRho/dS [kg m-3 ppt-1]
[in]stable_lLeft-column, top interface is stable
[in]pres_rRight-column interface pressure [Pa]
[in]hcol_rLeft-column layer thicknesses
[in]trRight-column top interface potential temperature [degC]
[in]srRight-column top interface salinity [ppt]
[in]drdt_rRight-column, top interface dRho/dT [kg m-3 degC-1]
[in]drds_rRight-column, top interface dRho/dS [kg m-3 ppt-1]
[in]stable_rRight-column, top interface is stable
[in,out]polFractional position of neutral surface within layer KoL of left column
[in,out]porFractional position of neutral surface within layer KoR of right column
[in,out]kolIndex of first left interface above neutral surface
[in,out]korIndex of first right interface above neutral surface
[in,out]heffEffective thickness between two neutral surfaces [Pa]
[in]ppoly_t_lLeft-column coefficients of T reconstruction
[in]ppoly_s_lLeft-column coefficients of S reconstruction
[in]ppoly_t_rRight-column coefficients of T reconstruction
[in]ppoly_s_rRight-column coefficients of S reconstruction

Definition at line 994 of file MOM_neutral_diffusion.F90.

994  type(neutral_diffusion_CS), intent(inout) :: CS !< Neutral diffusion control structure
995  integer, intent(in) :: nk !< Number of levels
996  integer, intent(in) :: ns !< Number of neutral surfaces
997  real, dimension(nk+1), intent(in) :: Pres_l !< Left-column interface pressure [Pa]
998  real, dimension(nk), intent(in) :: hcol_l !< Left-column layer thicknesses
999  real, dimension(nk,2), intent(in) :: Tl !< Left-column top interface potential temperature [degC]
1000  real, dimension(nk,2), intent(in) :: Sl !< Left-column top interface salinity [ppt]
1001  real, dimension(nk,2), intent(in) :: dRdT_l !< Left-column, top interface dRho/dT [kg m-3 degC-1]
1002  real, dimension(nk,2), intent(in) :: dRdS_l !< Left-column, top interface dRho/dS [kg m-3 ppt-1]
1003  logical, dimension(nk), intent(in) :: stable_l !< Left-column, top interface is stable
1004  real, dimension(nk+1), intent(in) :: Pres_r !< Right-column interface pressure [Pa]
1005  real, dimension(nk), intent(in) :: hcol_r !< Left-column layer thicknesses
1006  real, dimension(nk,2), intent(in) :: Tr !< Right-column top interface potential temperature [degC]
1007  real, dimension(nk,2), intent(in) :: Sr !< Right-column top interface salinity [ppt]
1008  real, dimension(nk,2), intent(in) :: dRdT_r !< Right-column, top interface dRho/dT [kg m-3 degC-1]
1009  real, dimension(nk,2), intent(in) :: dRdS_r !< Right-column, top interface dRho/dS [kg m-3 ppt-1]
1010  logical, dimension(nk), intent(in) :: stable_r !< Right-column, top interface is stable
1011  real, dimension(4*nk), intent(inout) :: PoL !< Fractional position of neutral surface within
1012  !! layer KoL of left column
1013  real, dimension(4*nk), intent(inout) :: PoR !< Fractional position of neutral surface within
1014  !! layer KoR of right column
1015  integer, dimension(4*nk), intent(inout) :: KoL !< Index of first left interface above neutral surface
1016  integer, dimension(4*nk), intent(inout) :: KoR !< Index of first right interface above neutral surface
1017  real, dimension(4*nk-1), intent(inout) :: hEff !< Effective thickness between two neutral surfaces [Pa]
1018  real, dimension(nk,CS%deg+1), &
1019  optional, intent(in) :: ppoly_T_l !< Left-column coefficients of T reconstruction
1020  real, dimension(nk,CS%deg+1), &
1021  optional, intent(in) :: ppoly_S_l !< Left-column coefficients of S reconstruction
1022  real, dimension(nk,CS%deg+1), &
1023  optional, intent(in) :: ppoly_T_r !< Right-column coefficients of T reconstruction
1024  real, dimension(nk,CS%deg+1), &
1025  optional, intent(in) :: ppoly_S_r !< Right-column coefficients of S reconstruction
1026 
1027  ! Local variables
1028  integer :: k_surface ! Index of neutral surface
1029  integer :: kl_left, kl_right ! Index of layers on the left/right
1030  integer :: ki_left, ki_right ! Index of interfaces on the left/right
1031  logical :: searching_left_column ! True if searching for the position of a right interface in the left column
1032  logical :: searching_right_column ! True if searching for the position of a left interface in the right column
1033  logical :: reached_bottom ! True if one of the bottom-most interfaces has been used as the target
1034  logical :: search_layer
1035  integer :: k, kl_left_0, kl_right_0
1036  real :: dRho, dRhoTop, dRhoBot, hL, hR
1037  integer :: lastK_left, lastK_right
1038  real :: lastP_left, lastP_right
1039  real :: min_bound
1040  real :: T_other, S_other, P_other, dRdT_other, dRdS_other
1041  logical, dimension(nk) :: top_connected_l, top_connected_r
1042  logical, dimension(nk) :: bot_connected_l, bot_connected_r
1043 
1044  top_connected_l(:) = .false. ; top_connected_r(:) = .false.
1045  bot_connected_l(:) = .false. ; bot_connected_r(:) = .false.
1046 
1047  ! Check to make sure that polynomial reconstructions were passed if refine_pos defined)
1048  if (cs%refine_position) then
1049  if (.not. ( present(ppoly_t_l) .and. present(ppoly_s_l) .and. &
1050  present(ppoly_t_r) .and. present(ppoly_s_r) ) ) &
1051  call mom_error(fatal, "fine_neutral_surface_positions_discontinuous: refine_pos is requested, but " //&
1052  "polynomial coefficients not available for T and S")
1053  endif
1054  do k = 1,nk
1055  if (stable_l(k)) then
1056  kl_left = k
1057  kl_left_0 = k
1058  exit
1059  endif
1060  enddo
1061  do k = 1,nk
1062  if (stable_r(k)) then
1063  kl_right = k
1064  kl_right_0 = k
1065  exit
1066  endif
1067  enddo
1068 
1069  ! Initialize variables for the search
1070  ki_right = 1 ; lastk_right = 1 ; lastp_right = -1.
1071  ki_left = 1 ; lastk_left = 1 ; lastp_left = -1.
1072 
1073  reached_bottom = .false.
1074  searching_left_column = .false.
1075  searching_right_column = .false.
1076 
1077  ! Loop over each neutral surface, working from top to bottom
1078  neutral_surfaces: do k_surface = 1, ns
1079  ! Potential density difference, rho(kr) - rho(kl)
1080  drho = 0.5 * ( ( drdt_r(kl_right,ki_right) + drdt_l(kl_left,ki_left) ) * &
1081  ( tr(kl_right,ki_right) - tl(kl_left,ki_left) ) &
1082  + ( drds_r(kl_right,ki_right) + drds_l(kl_left,ki_left) ) * &
1083  ( sr(kl_right,ki_right) - sl(kl_left,ki_left) ) )
1084  if (cs%debug) write(*,'(A,I2,A,E12.4,A,I2,A,I2,A,I2,A,I2)') "k_surface=",k_surface," dRho=",drho, &
1085  "kl_left=",kl_left, " ki_left=",ki_left," kl_right=",kl_right, " ki_right=",ki_right
1086  ! Which column has the lighter surface for the current indexes, kr and kl
1087  if (.not. reached_bottom) then
1088  if (drho < 0.) then
1089  searching_left_column = .true.
1090  searching_right_column = .false.
1091  elseif (drho > 0.) then
1092  searching_right_column = .true.
1093  searching_left_column = .false.
1094  else ! dRho == 0.
1095  if ( ( kl_left == kl_left_0) .and. ( kl_right == kl_right_0 ) .and. &
1096  (ki_left + ki_right == 2) ) then ! Still at surface
1097  searching_left_column = .true.
1098  searching_right_column = .false.
1099  else ! Not the surface so we simply change direction
1100  searching_left_column = .not. searching_left_column
1101  searching_right_column = .not. searching_right_column
1102  endif
1103  endif
1104  endif
1105 
1106  if (searching_left_column) then
1107  ! delta_rho is referenced to the right interface T, S, and P
1108  if (cs%ref_pres>=0.) then
1109  p_other = cs%ref_pres
1110  else
1111  if (ki_right == 1) p_other = pres_r(kl_right)
1112  if (ki_right == 2) p_other = pres_r(kl_right+1)
1113  endif
1114  t_other = tr(kl_right, ki_right)
1115  s_other = sr(kl_right, ki_right)
1116  drdt_other = drdt_r(kl_right, ki_right)
1117  drds_other = drds_r(kl_right, ki_right)
1118  if (cs%refine_position .and. (lastk_left == kl_left)) then
1119  call drho_at_pos(cs%ndiff_aux_CS, t_other, s_other, drdt_other, drds_other, pres_l(kl_left), &
1120  pres_l(kl_left+1), ppoly_t_l(kl_left,:), ppoly_s_l(kl_left,:), lastp_left, drhotop)
1121  else
1122  drhotop = calc_drho(tl(kl_left,1), sl(kl_left,1), drdt_l(kl_left,1), drds_l(kl_left,1), t_other, s_other, &
1123  drdt_other, drds_other)
1124  endif
1125  ! Potential density difference, rho(kl) - rho(kl_right,ki_right) (will be positive)
1126  drhobot = calc_drho(tl(kl_left,2), sl(kl_left,2), drdt_l(kl_left,2), drds_l(kl_left,2), &
1127  t_other, s_other, drdt_other, drds_other)
1128  if (cs%debug) then
1129  write(*,'(A,I2,A,E12.4,A,E12.4,A,E12.4)') "Searching left layer ", kl_left, &
1130  " dRhoTop=", drhotop, " dRhoBot=", drhobot
1131  write(*,'(A,I2,X,I2)') "Searching from right: ", kl_right, ki_right
1132  write(*,*) "Temp/Salt Reference: ", t_other, s_other
1133  write(*,*) "Temp/Salt Top L: ", tl(kl_left,1), sl(kl_left,1)
1134  write(*,*) "Temp/Salt Bot L: ", tl(kl_left,2), sl(kl_left,2)
1135  endif
1136 
1137  ! Set the position within the starting column
1138  por(k_surface) = real(ki_right-1)
1139  kor(k_surface) = kl_right
1140 
1141  ! Set position within the searched column
1142  call search_other_column(drhotop, drhobot, pres_l(kl_left), pres_l(kl_left+1), &
1143  lastp_left, lastk_left, kl_left, kl_left_0, ki_left, &
1144  top_connected_l, bot_connected_l, pol(k_surface), kol(k_surface), search_layer)
1145 
1146  if ( cs%refine_position .and. search_layer ) then
1147  min_bound = 0.
1148  if (k_surface > 1) then
1149  if ( kol(k_surface) == kol(k_surface-1) ) min_bound = pol(k_surface-1)
1150  endif
1151  pol(k_surface) = refine_nondim_position( cs%ndiff_aux_CS, t_other, s_other, drdt_other, drds_other, &
1152  pres_l(kl_left), pres_l(kl_left+1), ppoly_t_l(kl_left,:), ppoly_s_l(kl_left,:), &
1153  drhotop, drhobot, min_bound )
1154  endif
1155  if (pol(k_surface) == 0.) top_connected_l(kol(k_surface)) = .true.
1156  if (pol(k_surface) == 1.) bot_connected_l(kol(k_surface)) = .true.
1157  call increment_interface(nk, kl_right, ki_right, stable_r, reached_bottom, &
1158  searching_right_column, searching_left_column)
1159 
1160  elseif (searching_right_column) then
1161  if (cs%ref_pres>=0.) then
1162  p_other = cs%ref_pres
1163  else
1164  if (ki_left == 1) p_other = pres_l(kl_left)
1165  if (ki_left == 2) p_other = pres_l(kl_left+1)
1166  endif
1167  t_other = tl(kl_left, ki_left)
1168  s_other = sl(kl_left, ki_left)
1169  drdt_other = drdt_l(kl_left, ki_left)
1170  drds_other = drds_l(kl_left, ki_left)
1171  ! Interpolate for the neutral surface position within the right column, layer krm1
1172  ! Potential density difference, rho(kr-1) - rho(kl) (should be negative)
1173 
1174  if (cs%refine_position .and. (lastk_right == kl_right)) then
1175  call drho_at_pos(cs%ndiff_aux_CS, t_other, s_other, drdt_other, drds_other, pres_r(kl_right), &
1176  pres_l(kl_right+1), ppoly_t_r(kl_right,:), ppoly_s_r(kl_right,:), lastp_right, drhotop)
1177  else
1178  drhotop = calc_drho(tr(kl_right,1), sr(kl_right,1), drdt_r(kl_right,1), drds_r(kl_right,1), &
1179  t_other, s_other, drdt_other, drds_other)
1180  endif
1181  drhobot = calc_drho(tr(kl_right,2), sr(kl_right,2), drdt_r(kl_right,2), drds_r(kl_right,2), &
1182  t_other, s_other, drdt_other, drds_other)
1183  if (cs%debug) then
1184  write(*,'(A,I2,A,E12.4,A,E12.4,A,E12.4)') "Searching right layer ", kl_right, &
1185  " dRhoTop=", drhotop, " dRhoBot=", drhobot
1186  write(*,'(A,I2,X,I2)') "Searching from left: ", kl_left, ki_left
1187  write(*,*) "Temp/Salt Reference: ", t_other, s_other
1188  write(*,*) "Temp/Salt Top R: ", tr(kl_right,1), sr(kl_right,1)
1189  write(*,*) "Temp/Salt Bot R: ", tr(kl_right,2), sr(kl_right,2)
1190  endif
1191  ! Set the position within the starting column
1192  pol(k_surface) = real(ki_left-1)
1193  kol(k_surface) = kl_left
1194 
1195  ! Set position within the searched column
1196  call search_other_column(drhotop, drhobot, pres_r(kl_right), pres_r(kl_right+1), lastp_right, lastk_right, &
1197  kl_right, kl_right_0, ki_right, top_connected_r, bot_connected_r, por(k_surface), &
1198  kor(k_surface), search_layer)
1199  if ( cs%refine_position .and. search_layer) then
1200  min_bound = 0.
1201  if (k_surface > 1) then
1202  if ( kor(k_surface) == kor(k_surface-1) ) min_bound = por(k_surface-1)
1203  endif
1204  por(k_surface) = refine_nondim_position(cs%ndiff_aux_CS, t_other, s_other, drdt_other, drds_other, &
1205  pres_r(kl_right), pres_r(kl_right+1), ppoly_t_r(kl_right,:), ppoly_s_r(kl_right,:), &
1206  drhotop, drhobot, min_bound )
1207  endif
1208  if (por(k_surface) == 0.) top_connected_r(kor(k_surface)) = .true.
1209  if (por(k_surface) == 1.) bot_connected_r(kor(k_surface)) = .true.
1210  call increment_interface(nk, kl_left, ki_left, stable_l, reached_bottom, &
1211  searching_left_column, searching_right_column)
1212 
1213  else
1214  stop 'Else what?'
1215  endif
1216  lastk_left = kol(k_surface) ; lastp_left = pol(k_surface)
1217  lastk_right = kor(k_surface) ; lastp_right = por(k_surface)
1218 
1219  if (cs%debug) write(*,'(A,I3,A,ES16.6,A,I2,A,ES16.6)') "KoL:", kol(k_surface), " PoL:", pol(k_surface), &
1220  " KoR:", kor(k_surface), " PoR:", por(k_surface)
1221  ! Effective thickness
1222  if (k_surface>1) then
1223  ! This is useful as a check to make sure that positions are monotonically increasing
1224  hl = absolute_position(nk,ns,pres_l,kol,pol,k_surface) - absolute_position(nk,ns,pres_l,kol,pol,k_surface-1)
1225  hr = absolute_position(nk,ns,pres_r,kor,por,k_surface) - absolute_position(nk,ns,pres_r,kor,por,k_surface-1)
1226  ! In the case of a layer being unstably stratified, may get a negative thickness. Set the previous position
1227  ! to the current location
1228  if ( hl<0. .or. hr<0. ) then
1229  heff(k_surface-1) = 0.
1230  call mom_error(warning, "hL or hR is negative")
1231  elseif ( hl > 0. .and. hr > 0.) then
1232  hl = (pol(k_surface) - pol(k_surface-1))*hcol_l(kol(k_surface))
1233  hr = (por(k_surface) - por(k_surface-1))*hcol_r(kor(k_surface))
1234  heff(k_surface-1) = 2. * ( (hl * hr) / ( hl + hr ) )! Harmonic mean
1235  else
1236  heff(k_surface-1) = 0.
1237  endif
1238  endif
1239  enddo neutral_surfaces
1240  if (cs%debug) then
1241  write (*,*) "==========Start Neutral Surfaces=========="
1242  do k = 1,ns-1
1243  if (heff(k)>0.) then
1244  kl_left = kol(k)
1245  kl_right = kor(k)
1246  write (*,'(A,I3,X,ES16.6,X,I3,X,ES16.6)') "Top surface KoL, PoL, KoR, PoR: ", kl_left, pol(k), kl_right, por(k)
1247  call check_neutral_positions(cs%ndiff_aux_CS, pres_l(kl_left), pres_l(kl_left+1), pres_r(kl_right), &
1248  pres_r(kl_right+1), pol(k), por(k), ppoly_t_l(kl_left,:), ppoly_t_r(kl_right,:), &
1249  ppoly_s_l(kl_left,:), ppoly_s_r(kl_right,:))
1250  kl_left = kol(k+1)
1251  kl_right = kor(k+1)
1252  write (*,'(A,I3,X,ES16.6,X,I3,X,ES16.6)') "Bot surface KoL, PoL, KoR, PoR: ", kl_left, pol(k+1), kl_right, por(k)
1253  call check_neutral_positions(cs%ndiff_aux_CS, pres_l(kl_left), pres_l(kl_left+1), pres_r(kl_right), &
1254  pres_r(kl_right+1), pol(k), por(k), ppoly_t_l(kl_left,:), ppoly_t_r(kl_right,:), &
1255  ppoly_s_l(kl_left,:), ppoly_s_r(kl_right,:))
1256  endif
1257  enddo
1258  write(*,'(A,E16.6)') "Total thickness of sublayers: ", sum(heff)
1259  write(*,*) "==========End Neutral Surfaces=========="
1260  endif
1261 

◆ fv_diff()

real function mom_neutral_diffusion::fv_diff ( real, intent(in)  hkm1,
real, intent(in)  hk,
real, intent(in)  hkp1,
real, intent(in)  Skm1,
real, intent(in)  Sk,
real, intent(in)  Skp1 
)
private

Returns the cell-centered second-order finite volume (unlimited PLM) slope using three consecutive cell widths and average values. Slope is returned as a difference across the central cell (i.e. units of scalar S). Discretization follows equation 1.7 in Colella & Woodward, 1984: JCP 54, 174-201.

Parameters
[in]hkm1Left cell width
[in]hkCenter cell width
[in]hkp1Right cell width
[in]skm1Left cell average value
[in]skCenter cell average value
[in]skp1Right cell average value

Definition at line 754 of file MOM_neutral_diffusion.F90.

754  real, intent(in) :: hkm1 !< Left cell width
755  real, intent(in) :: hk !< Center cell width
756  real, intent(in) :: hkp1 !< Right cell width
757  real, intent(in) :: Skm1 !< Left cell average value
758  real, intent(in) :: Sk !< Center cell average value
759  real, intent(in) :: Skp1 !< Right cell average value
760 
761  ! Local variables
762  real :: h_sum, hp, hm
763 
764  h_sum = ( hkm1 + hkp1 ) + hk
765  if (h_sum /= 0.) h_sum = 1./ h_sum
766  hm = hkm1 + hk
767  if (hm /= 0.) hm = 1./ hm
768  hp = hkp1 + hk
769  if (hp /= 0.) hp = 1./ hp
770  fv_diff = ( hk * h_sum ) * &
771  ( ( 2. * hkm1 + hk ) * hp * ( skp1 - sk ) &
772  + ( 2. * hkp1 + hk ) * hm * ( sk - skm1 ) )

◆ fvlsq_slope()

real function mom_neutral_diffusion::fvlsq_slope ( real, intent(in)  hkm1,
real, intent(in)  hk,
real, intent(in)  hkp1,
real, intent(in)  Skm1,
real, intent(in)  Sk,
real, intent(in)  Skp1 
)
private

Returns the cell-centered second-order weighted least squares slope using three consecutive cell widths and average values. Slope is returned as a gradient (i.e. units of scalar S over width units).

Parameters
[in]hkm1Left cell width
[in]hkCenter cell width
[in]hkp1Right cell width
[in]skm1Left cell average value
[in]skCenter cell average value
[in]skp1Right cell average value

Definition at line 780 of file MOM_neutral_diffusion.F90.

780  real, intent(in) :: hkm1 !< Left cell width
781  real, intent(in) :: hk !< Center cell width
782  real, intent(in) :: hkp1 !< Right cell width
783  real, intent(in) :: Skm1 !< Left cell average value
784  real, intent(in) :: Sk !< Center cell average value
785  real, intent(in) :: Skp1 !< Right cell average value
786 
787  ! Local variables
788  real :: xkm1, xkp1
789  real :: h_sum, hx_sum, hxsq_sum, hxy_sum, hy_sum, det
790 
791  xkm1 = -0.5 * ( hk + hkm1 )
792  xkp1 = 0.5 * ( hk + hkp1 )
793  h_sum = ( hkm1 + hkp1 ) + hk
794  hx_sum = hkm1*xkm1 + hkp1*xkp1
795  hxsq_sum = hkm1*(xkm1**2) + hkp1*(xkp1**2)
796  hxy_sum = hkm1*xkm1*skm1 + hkp1*xkp1*skp1
797  hy_sum = ( hkm1*skm1 + hkp1*skp1 ) + hk*sk
798  det = h_sum * hxsq_sum - hx_sum**2
799  if (det /= 0.) then
800  !a = ( hxsq_sum * hy_sum - hx_sum*hxy_sum ) / det ! a would be mean of straight line fit
801  fvlsq_slope = ( h_sum * hxy_sum - hx_sum*hy_sum ) / det ! Gradient of straight line fit
802  else
803  fvlsq_slope = 0. ! Adcroft's reciprocal rule
804  endif

◆ interface_scalar()

subroutine mom_neutral_diffusion::interface_scalar ( integer, intent(in)  nk,
real, dimension(nk), intent(in)  h,
real, dimension(nk), intent(in)  S,
real, dimension(nk+1), intent(inout)  Si,
integer, intent(in)  i_method,
real, intent(in)  h_neglect 
)
private

Returns interface scalar, Si, for a column of layer values, S.

Parameters
[in]nkNumber of levels
[in]hLayer thickness [H ~> m or kg m-2]
[in]sLayer scalar (conc, e.g. ppt)
[in,out]siInterface scalar (conc, e.g. ppt)
[in]i_method=1 use average of PLM edges =2 use continuous PPM edge interpolation
[in]h_neglectA negligibly small thickness [H ~> m or kg m-2]

Definition at line 571 of file MOM_neutral_diffusion.F90.

571  integer, intent(in) :: nk !< Number of levels
572  real, dimension(nk), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
573  real, dimension(nk), intent(in) :: S !< Layer scalar (conc, e.g. ppt)
574  real, dimension(nk+1), intent(inout) :: Si !< Interface scalar (conc, e.g. ppt)
575  integer, intent(in) :: i_method !< =1 use average of PLM edges
576  !! =2 use continuous PPM edge interpolation
577  real, intent(in) :: h_neglect !< A negligibly small thickness [H ~> m or kg m-2]
578  ! Local variables
579  integer :: k, km2, kp1
580  real, dimension(nk) :: diff
581  real :: Sb, Sa
582 
583  call plm_diff(nk, h, s, 2, 1, diff)
584  si(1) = s(1) - 0.5 * diff(1)
585  if (i_method==1) then
586  do k = 2, nk
587  ! Average of the two edge values (will be bounded and,
588  ! when slopes are unlimited, notionally second-order accurate)
589  sa = s(k-1) + 0.5 * diff(k-1) ! Lower edge value of a PLM reconstruction for layer above
590  sb = s(k) - 0.5 * diff(k) ! Upper edge value of a PLM reconstruction for layer below
591  si(k) = 0.5 * ( sa + sb )
592  enddo
593  elseif (i_method==2) then
594  do k = 2, nk
595  ! PPM quasi-fourth order interpolation for edge values following
596  ! equation 1.6 in Colella & Woodward, 1984: JCP 54, 174-201.
597  km2 = max(1, k-2)
598  kp1 = min(nk, k+1)
599  si(k) = ppm_edge(h(km2), h(k-1), h(k), h(kp1), s(k-1), s(k), diff(k-1), diff(k), h_neglect)
600  enddo
601  endif
602  si(nk+1) = s(nk) + 0.5 * diff(nk)
603 

◆ ndiff_unit_tests_continuous()

logical function mom_neutral_diffusion::ndiff_unit_tests_continuous ( logical, intent(in)  verbose)
private

Returns true if unit tests of neutral_diffusion functions fail. Otherwise returns false.

Parameters
[in]verboseIf true, write results to stdout

Definition at line 1517 of file MOM_neutral_diffusion.F90.

1517  logical, intent(in) :: verbose !< If true, write results to stdout
1518  ! Local variables
1519  integer, parameter :: nk = 4
1520  real, dimension(nk+1) :: TiL, TiR1, TiR2, TiR4, Tio ! Test interface temperatures
1521  real, dimension(nk) :: TL ! Test layer temperatures
1522  real, dimension(nk+1) :: SiL ! Test interface salinities
1523  real, dimension(nk+1) :: PiL, PiR4 ! Test interface positions
1524  real, dimension(2*nk+2) :: PiLRo, PiRLo ! Test positions
1525  integer, dimension(2*nk+2) :: KoL, KoR ! Test indexes
1526  real, dimension(2*nk+1) :: hEff ! Test positions
1527  real, dimension(2*nk+1) :: Flx ! Test flux
1528  integer :: k
1529  logical :: v
1530  real :: h_neglect, h_neglect_edge
1531 
1532  h_neglect_edge = 1.0e-10 ; h_neglect = 1.0e-30
1533 
1534  v = verbose
1535 
1536  ndiff_unit_tests_continuous = .false. ! Normally return false
1537  write(*,*) '==== MOM_neutral_diffusion: ndiff_unit_tests_continuous ='
1538 
1539  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
1540  test_fv_diff(v,1.,1.,1., 0.,1.,2., 1., 'FV: Straight line on uniform grid')
1541  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
1542  test_fv_diff(v,1.,1.,0., 0.,4.,8., 7., 'FV: Vanished right cell')
1543  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
1544  test_fv_diff(v,0.,1.,1., 0.,4.,8., 7., 'FV: Vanished left cell')
1545  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
1546  test_fv_diff(v,1.,2.,4., 0.,3.,9., 4., 'FV: Stretched grid')
1547  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
1548  test_fv_diff(v,2.,0.,2., 0.,1.,2., 0., 'FV: Vanished middle cell')
1549  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
1550  test_fv_diff(v,0.,1.,0., 0.,1.,2., 2., 'FV: Vanished on both sides')
1551  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
1552  test_fv_diff(v,1.,0.,0., 0.,1.,2., 0., 'FV: Two vanished cell sides')
1553  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
1554  test_fv_diff(v,0.,0.,0., 0.,1.,2., 0., 'FV: All vanished cells')
1555 
1556  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
1557  test_fvlsq_slope(v,1.,1.,1., 0.,1.,2., 1., 'LSQ: Straight line on uniform grid')
1558  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
1559  test_fvlsq_slope(v,1.,1.,0., 0.,1.,2., 1., 'LSQ: Vanished right cell')
1560  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
1561  test_fvlsq_slope(v,0.,1.,1., 0.,1.,2., 1., 'LSQ: Vanished left cell')
1562  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
1563  test_fvlsq_slope(v,1.,2.,4., 0.,3.,9., 2., 'LSQ: Stretched grid')
1564  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
1565  test_fvlsq_slope(v,1.,0.,1., 0.,1.,2., 2., 'LSQ: Vanished middle cell')
1566  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
1567  test_fvlsq_slope(v,0.,1.,0., 0.,1.,2., 0., 'LSQ: Vanished on both sides')
1568  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
1569  test_fvlsq_slope(v,1.,0.,0., 0.,1.,2., 0., 'LSQ: Two vanished cell sides')
1570  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
1571  test_fvlsq_slope(v,0.,0.,0., 0.,1.,2., 0., 'LSQ: All vanished cells')
1572 
1573  call interface_scalar(4, (/10.,10.,10.,10./), (/24.,18.,12.,6./), tio, 1, h_neglect)
1574  !ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
1575  ! test_data1d(5, Tio, (/27.,21.,15.,9.,3./), 'Linear profile, interface temperatures')
1576  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
1577  test_data1d(v,5, tio, (/24.,22.5,15.,7.5,6./), 'Linear profile, linear interface temperatures')
1578  call interface_scalar(4, (/10.,10.,10.,10./), (/24.,18.,12.,6./), tio, 2, h_neglect)
1579  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
1580  test_data1d(v,5, tio, (/24.,22.,15.,8.,6./), 'Linear profile, PPM interface temperatures')
1581 
1582  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
1583  test_ifndp(v,-1.0, 0., 1.0, 1.0, 0.5, 'Check mid-point')
1584  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
1585  test_ifndp(v, 0.0, 0., 1.0, 1.0, 0.0, 'Check bottom')
1586  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
1587  test_ifndp(v, 0.1, 0., 1.1, 1.0, 0.0, 'Check below')
1588  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
1589  test_ifndp(v,-1.0, 0., 0.0, 1.0, 1.0, 'Check top')
1590  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
1591  test_ifndp(v,-1.0, 0., -0.1, 1.0, 1.0, 'Check above')
1592  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
1593  test_ifndp(v,-1.0, 0., 3.0, 1.0, 0.25, 'Check 1/4')
1594  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
1595  test_ifndp(v,-3.0, 0., 1.0, 1.0, 0.75, 'Check 3/4')
1596  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
1597  test_ifndp(v, 1.0, 0., 1.0, 1.0, 0.0, 'Check dRho=0 below')
1598  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
1599  test_ifndp(v,-1.0, 0., -1.0, 1.0, 1.0, 'Check dRho=0 above')
1600  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
1601  test_ifndp(v, 0.0, 0., 0.0, 1.0, 0.5, 'Check dRho=0 mid')
1602  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
1603  test_ifndp(v,-2.0, .5, 5.0, 0.5, 0.5, 'Check dP=0')
1604 
1605  ! Identical columns
1606  call find_neutral_surface_positions_continuous(3, &
1607  (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Left positions, T and S
1608  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
1609  (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Right positions, T and S
1610  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
1611  pilro, pirlo, kol, kor, heff)
1612  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, kol, kor, pilro, pirlo, heff, &
1613  (/1,1,2,2,3,3,3,3/), & ! KoL
1614  (/1,1,2,2,3,3,3,3/), & ! KoR
1615  (/0.,0.,0.,0.,0.,0.,1.,1./), & ! pL
1616  (/0.,0.,0.,0.,0.,0.,1.,1./), & ! pR
1617  (/0.,10.,0.,10.,0.,10.,0./), & ! hEff
1618  'Identical columns')
1619  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_data1d(v, 8, &
1620  absolute_positions(3, 8, (/0.,10.,20.,30./), kol, pilro), &
1621  (/0.,0.,10.,10.,20.,20.,30.,30./), '... left positions')
1622  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_data1d(v, 8, &
1623  absolute_positions(3, 8, (/0.,10.,20.,30./), kor, pirlo), &
1624  (/0.,0.,10.,10.,20.,20.,30.,30./), '... right positions')
1625  call neutral_surface_flux(3, 2*3+2, 2, (/10.,10.,10./), (/10.,10.,10./), & ! nk, hL, hR
1626  (/20.,16.,12./), (/20.,16.,12./), & ! Tl, Tr
1627  pilro, pirlo, kol, kor, heff, flx, .true., &
1628  h_neglect, h_neglect_edge=h_neglect_edge)
1629  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_data1d(v, 7, flx, &
1630  (/0.,0.,0.,0.,0.,0.,0./), 'Identical columns, rho flux (=0)')
1631  call neutral_surface_flux(3, 2*3+2, 2, (/10.,10.,10./), (/10.,10.,10./), & ! nk, hL, hR
1632  (/-1.,-1.,-1./), (/1.,1.,1./), & ! Sl, Sr
1633  pilro, pirlo, kol, kor, heff, flx, .true., &
1634  h_neglect, h_neglect_edge=h_neglect_edge)
1635  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_data1d(v, 7, flx, &
1636  (/0.,20.,0.,20.,0.,20.,0./), 'Identical columns, S flux')
1637 
1638  ! Right column slightly cooler than left
1639  call find_neutral_surface_positions_continuous(3, &
1640  (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Left positions, T and S
1641  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
1642  (/0.,10.,20.,30./), (/20.,16.,12.,8./), (/0.,0.,0.,0./), & ! Right positions, T and S
1643  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
1644  pilro, pirlo, kol, kor, heff)
1645  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, kol, kor, pilro, pirlo, heff, &
1646  (/1,1,2,2,3,3,3,3/), & ! kL
1647  (/1,1,1,2,2,3,3,3/), & ! kR
1648  (/0.,0.5,0.,0.5,0.,0.5,1.,1./), & ! pL
1649  (/0.,0.,0.5,0.,0.5,0.,0.5,1./), & ! pR
1650  (/0.,5.,5.,5.,5.,5.,0./), & ! hEff
1651  'Right column slightly cooler')
1652  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_data1d(v, 8, &
1653  absolute_positions(3, 8, (/0.,10.,20.,30./), kol, pilro), &
1654  (/0.,5.,10.,15.,20.,25.,30.,30./), '... left positions')
1655  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_data1d(v, 8, &
1656  absolute_positions(3, 8, (/0.,10.,20.,30./), kor, pirlo), &
1657  (/0.,0.,5.,10.,15.,20.,25.,30./), '... right positions')
1658 
1659  ! Right column slightly warmer than left
1660  call find_neutral_surface_positions_continuous(3, &
1661  (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Left positions, T and S
1662  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
1663  (/0.,10.,20.,30./), (/24.,20.,16.,12./), (/0.,0.,0.,0./), & ! Right positions, T and S
1664  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
1665  pilro, pirlo, kol, kor, heff)
1666  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, kol, kor, pilro, pirlo, heff, &
1667  (/1,1,1,2,2,3,3,3/), & ! kL
1668  (/1,1,2,2,3,3,3,3/), & ! kR
1669  (/0.,0.,0.5,0.,0.5,0.,0.5,1./), & ! pL
1670  (/0.,0.5,0.,0.5,0.,0.5,1.,1./), & ! pR
1671  (/0.,5.,5.,5.,5.,5.,0./), & ! hEff
1672  'Right column slightly warmer')
1673 
1674  ! Right column somewhat cooler than left
1675  call find_neutral_surface_positions_continuous(3, &
1676  (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Left positions, T and S
1677  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
1678  (/0.,10.,20.,30./), (/16.,12.,8.,4./), (/0.,0.,0.,0./), & ! Right positions, T and S
1679  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
1680  pilro, pirlo, kol, kor, heff)
1681  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, kol, kor, pilro, pirlo, heff, &
1682  (/1,2,2,3,3,3,3,3/), & ! kL
1683  (/1,1,1,1,2,2,3,3/), & ! kR
1684  (/0.,0.,0.5,0.,0.5,1.,1.,1./), & ! pL
1685  (/0.,0.,0.,0.5,0.,0.5,0.,1./), & ! pR
1686  (/0.,0.,5.,5.,5.,0.,0./), & ! hEff
1687  'Right column somewhat cooler')
1688 
1689  ! Right column much colder than left with no overlap
1690  call find_neutral_surface_positions_continuous(3, &
1691  (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Left positions, T and S
1692  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
1693  (/0.,10.,20.,30./), (/9.,7.,5.,3./), (/0.,0.,0.,0./), & ! Right positions, T and S
1694  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
1695  pilro, pirlo, kol, kor, heff)
1696  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, kol, kor, pilro, pirlo, heff, &
1697  (/1,2,3,3,3,3,3,3/), & ! kL
1698  (/1,1,1,1,1,2,3,3/), & ! kR
1699  (/0.,0.,0.,1.,1.,1.,1.,1./), & ! pL
1700  (/0.,0.,0.,0.,0.,0.,0.,1./), & ! pR
1701  (/0.,0.,0.,0.,0.,0.,0./), & ! hEff
1702  'Right column much cooler')
1703 
1704  ! Right column with mixed layer
1705  call find_neutral_surface_positions_continuous(3, &
1706  (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Left positions, T and S
1707  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
1708  (/0.,10.,20.,30./), (/14.,14.,10.,2./), (/0.,0.,0.,0./), & ! Right positions, T and S
1709  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
1710  pilro, pirlo, kol, kor, heff)
1711  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, kol, kor, pilro, pirlo, heff, &
1712  (/1,2,3,3,3,3,3,3/), & ! kL
1713  (/1,1,1,1,2,3,3,3/), & ! kR
1714  (/0.,0.,0.,0.,0.,1.,1.,1./), & ! pL
1715  (/0.,0.,0.,0.,0.,0.,0.,1./), & ! pR
1716  (/0.,0.,0.,0.,10.,0.,0./), & ! hEff
1717  'Right column with mixed layer')
1718 
1719  ! Identical columns with mixed layer
1720  call find_neutral_surface_positions_continuous(3, &
1721  (/0.,10.,20.,30./), (/14.,14.,10.,2./), (/0.,0.,0.,0./), & ! Left positions, T and S
1722  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
1723  (/0.,10.,20.,30./), (/14.,14.,10.,2./), (/0.,0.,0.,0./), & ! Right positions, T and S
1724  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
1725  pilro, pirlo, kol, kor, heff)
1726  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, kol, kor, pilro, pirlo, heff, &
1727  (/1,1,2,2,3,3,3,3/), & ! kL
1728  (/1,1,2,2,3,3,3,3/), & ! kR
1729  (/0.,0.,0.,0.,0.,0.,1.,1./), & ! pL
1730  (/0.,0.,0.,0.,0.,0.,1.,1./), & ! pR
1731  (/0.,10.,0.,10.,0.,10.,0./), & ! hEff
1732  'Indentical columns with mixed layer')
1733 
1734  ! Right column with unstable mixed layer
1735  call find_neutral_surface_positions_continuous(3, &
1736  (/0.,10.,20.,30./), (/14.,14.,10.,2./), (/0.,0.,0.,0./), & ! Left positions, T and S
1737  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
1738  (/0.,10.,20.,30./), (/10.,14.,12.,4./), (/0.,0.,0.,0./), & ! Right positions, T and S
1739  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
1740  pilro, pirlo, kol, kor, heff)
1741  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, kol, kor, pilro, pirlo, heff, &
1742  (/1,2,3,3,3,3,3,3/), & ! kL
1743  (/1,1,1,2,3,3,3,3/), & ! kR
1744  (/0.,0.,0.,0.,0.,0.,.75,1./), & ! pL
1745  (/0.,0.,0.,0.,0.,0.25,1.,1./), & ! pR
1746  (/0.,0.,0.,0.,0.,7.5,0./), & ! hEff
1747  'Right column with unstable mixed layer')
1748 
1749  ! Left column with unstable mixed layer
1750  call find_neutral_surface_positions_continuous(3, &
1751  (/0.,10.,20.,30./), (/10.,14.,12.,4./), (/0.,0.,0.,0./), & ! Left positions, T and S
1752  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
1753  (/0.,10.,20.,30./), (/14.,14.,10.,2./), (/0.,0.,0.,0./), & ! Right positions, T and S
1754  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
1755  pilro, pirlo, kol, kor, heff)
1756  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, kol, kor, pilro, pirlo, heff, &
1757  (/1,1,1,2,3,3,3,3/), & ! kL
1758  (/1,2,3,3,3,3,3,3/), & ! kR
1759  (/0.,0.,0.,0.,0.,0.25,1.,1./), & ! pL
1760  (/0.,0.,0.,0.,0.,0.,.75,1./), & ! pR
1761  (/0.,0.,0.,0.,0.,7.5,0./), & ! hEff
1762  'Left column with unstable mixed layer')
1763 
1764  ! Two unstable mixed layers
1765  call find_neutral_surface_positions_continuous(3, &
1766  (/0.,10.,20.,30./), (/8.,12.,10.,2./), (/0.,0.,0.,0./), & ! Left positions, T and S
1767  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
1768  (/0.,10.,20.,30./), (/10.,14.,12.,4./), (/0.,0.,0.,0./), & ! Right positions, T and S
1769  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
1770  pilro, pirlo, kol, kor, heff)
1771  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, kol, kor, pilro, pirlo, heff, &
1772  (/1,1,1,1,2,3,3,3/), & ! kL
1773  (/1,2,3,3,3,3,3,3/), & ! kR
1774  (/0.,0.,0.,0.,0.,0.,0.75,1./), & ! pL
1775  (/0.,0.,0.,0.5,0.5,0.5,1.,1./), & ! pR
1776  (/0.,0.,0.,0.,0.,6.,0./), & ! hEff
1777  'Two unstable mixed layers')
1778 
1779  if (.not. ndiff_unit_tests_continuous) write(*,*) 'Pass'
1780 

◆ ndiff_unit_tests_discontinuous()

logical function mom_neutral_diffusion::ndiff_unit_tests_discontinuous ( logical, intent(in)  verbose)
private
Parameters
[in]verboseIt true, write results to stdout

Definition at line 1784 of file MOM_neutral_diffusion.F90.

1784  logical, intent(in) :: verbose !< It true, write results to stdout
1785  ! Local variables
1786  integer, parameter :: nk = 3
1787  integer, parameter :: ns = nk*4
1788  real, dimension(nk) :: Sl, Sr, Tl, Tr, hl, hr
1789  real, dimension(nk,2) :: TiL, SiL, TiR, SiR
1790  real, dimension(nk+1) :: Pres_l, Pres_R
1791  integer, dimension(ns) :: KoL, KoR
1792  real, dimension(ns) :: PoL, PoR
1793  real, dimension(ns-1) :: hEff, Flx
1794  type(neutral_diffusion_CS) :: CS !< Neutral diffusion control structure
1795  type(EOS_type), pointer :: EOS !< Structure for linear equation of state
1796  type(remapping_CS), pointer :: remap_CS !< Remapping control structure (PLM)
1797  real, dimension(nk,2) :: poly_T_l, poly_T_r, poly_S, poly_slope ! Linear reconstruction for T
1798  real, dimension(nk,2) :: dRdT, dRdS
1799  logical, dimension(nk) :: stable_l, stable_r
1800  integer :: iMethod
1801  integer :: ns_l, ns_r
1802  real :: h_neglect, h_neglect_edge
1803  integer :: k
1804  logical :: v
1805 
1806  v = verbose
1807  ndiff_unit_tests_discontinuous = .false. ! Normally return false
1808  write(*,*) '==== MOM_neutral_diffusion: ndiff_unit_tests_discontinuous ='
1809 
1810  h_neglect = 1.0e-30 ; h_neglect_edge = 1.0e-10
1811 
1812  ! Unit tests for find_neutral_surface_positions_discontinuous
1813  ! Salinity is 0 for all these tests
1814  sl(:) = 0. ; sr(:) = 0. ; poly_s(:,:) = 0. ; sil(:,:) = 0. ; sir(:,:) = 0.
1815  drdt(:,:) = -1. ; drds(:,:) = 0.
1816 
1817  ! Intialize any control structures needed for unit tests
1818  cs%refine_position = .false.
1819  cs%ref_pres = -1.
1820  allocate(remap_cs)
1821  call initialize_remapping( remap_cs, "PLM", boundary_extrapolation = .true. )
1822 
1823  hl = (/10.,10.,10./) ; hr = (/10.,10.,10./) ; pres_l(1) = 0. ; pres_r(1) = 0.
1824  do k = 1,nk ; pres_l(k+1) = pres_l(k) + hl(k) ; pres_r(k+1) = pres_r(k) + hr(k) ; enddo
1825  ! Identical columns
1826  tl = (/20.,16.,12./) ; tr = (/20.,16.,12./)
1827  call build_reconstructions_1d( remap_cs, nk, hl, tl, poly_t_l, til, poly_slope, imethod, h_neglect, h_neglect_edge )
1828  call build_reconstructions_1d( remap_cs, nk, hr, tr, poly_t_r, tir, poly_slope, imethod, h_neglect, h_neglect_edge )
1829  call mark_unstable_cells( nk, drdt, drds, til, sil, stable_l, ns_l )
1830  call mark_unstable_cells( nk, drdt, drds, tir, sir, stable_r, ns_r )
1831  call find_neutral_surface_positions_discontinuous(cs, nk, ns_l+ns_r, pres_l, hl, til, sil, drdt, drds, stable_l, &
1832  pres_r, hr, tir, sir, drdt, drds, stable_r, pol, por, kol, kor, heff)
1833  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, kol, kor, pol, por, heff, &
1834  (/1,1,1,1,2,2,2,2,3,3,3,3/), & ! KoL
1835  (/1,1,1,1,2,2,2,2,3,3,3,3/), & ! KoR
1836  (/0.,0.,1.,1.,0.,0.,1.,1.,0.,0.,1.,1./), & ! pL
1837  (/0.,0.,1.,1.,0.,0.,1.,1.,0.,0.,1.,1./), & ! pR
1838  (/0.,10.,0.,0.,0.,10.,0.,0.,0.,10.,0.,0./), & ! hEff
1839  'Identical columns')
1840  tl = (/20.,16.,12./) ; tr = (/18.,14.,10./)
1841  call build_reconstructions_1d( remap_cs, nk, hl, tl, poly_t_l, til, poly_slope, imethod, h_neglect, h_neglect_edge )
1842  call build_reconstructions_1d( remap_cs, nk, hr, tr, poly_t_r, tir, poly_slope, imethod, h_neglect, h_neglect_edge )
1843  call mark_unstable_cells( nk, drdt, drds, til, sil, stable_l, ns_l )
1844  call mark_unstable_cells( nk, drdt, drds, tir, sir, stable_r, ns_r )
1845  call find_neutral_surface_positions_discontinuous(cs, nk, ns_l+ns_r, pres_l, hl, til, sil, drdt, drds, stable_l, &
1846  pres_r, hr, tir, sir, drdt, drds, stable_r, pol, por, kol, kor, heff)
1847  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, kol, kor, pol, por, heff, &
1848  (/1,1,1,2,2,2,2,3,3,3,3,3/), & ! KoL
1849  (/1,1,1,1,1,2,2,2,2,3,3,3/), & ! KoR
1850  (/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
1851  (/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
1852  (/0.0, 5.0, 0.0, 5.0, 0.0, 5.0, 0.0, 5.0, 0.0, 5.0, 0.0/), & ! hEff
1853  'Right column slightly cooler')
1854  tl = (/18.,14.,10./) ; tr = (/20.,16.,12./)
1855  call build_reconstructions_1d( remap_cs, nk, hl, tl, poly_t_l, til, poly_slope, imethod, h_neglect, h_neglect_edge )
1856  call build_reconstructions_1d( remap_cs, nk, hr, tr, poly_t_r, tir, poly_slope, imethod, h_neglect, h_neglect_edge )
1857  call mark_unstable_cells( nk, drdt, drds, til, sil, stable_l, ns_l )
1858  call mark_unstable_cells( nk, drdt, drds, tir, sir, stable_r, ns_r )
1859  call find_neutral_surface_positions_discontinuous(cs, nk, ns_l+ns_r, pres_l, hl, til, sil, drdt, drds, stable_l, &
1860  pres_r, hr, tir, sir, drdt, drds, stable_r, pol, por, kol, kor, heff)
1861  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, kol, kor, pol, por, heff, &
1862  (/1,1,1,1,1,2,2,2,2,3,3,3/), & ! KoL
1863  (/1,1,1,2,2,2,2,3,3,3,3,3/), & ! KoR
1864  (/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
1865  (/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
1866  (/0.0, 5.0, 0.0, 5.0, 0.0, 5.0, 0.0, 5.0, 0.0, 5.0, 0.0/), & ! hEff
1867  'Left column slightly cooler')
1868  tl = (/20.,16.,12./) ; tr = (/14.,10.,6./)
1869  call build_reconstructions_1d( remap_cs, nk, hl, tl, poly_t_l, til, poly_slope, imethod, h_neglect, h_neglect_edge )
1870  call build_reconstructions_1d( remap_cs, nk, hr, tr, poly_t_r, tir, poly_slope, imethod, h_neglect, h_neglect_edge )
1871  call mark_unstable_cells( nk, drdt, drds, til, sil, stable_l, ns_l )
1872  call mark_unstable_cells( nk, drdt, drds, tir, sir, stable_r, ns_r )
1873  call find_neutral_surface_positions_discontinuous(cs, nk, ns_l+ns_r, pres_l, hl, til, sil, drdt, drds, stable_l, &
1874  pres_r, hr, tir, sir, drdt, drds, stable_r, pol, por, kol, kor, heff)
1875  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, kol, kor, pol, por, heff, &
1876  (/1,1,2,2,2,3,3,3,3,3,3,3/), & ! KoL
1877  (/1,1,1,1,1,1,1,2,2,2,3,3/), & ! KoR
1878  (/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
1879  (/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
1880  (/0.0, 0.0, 0.0, 5.0, 0.0, 5.0, 0.0, 5.0, 0.0, 0.0, 0.0/), & ! hEff
1881  'Right column somewhat cooler')
1882  tl = (/20.,16.,12./) ; tr = (/8.,6.,4./)
1883  call build_reconstructions_1d( remap_cs, nk, hl, tl, poly_t_l, til, poly_slope, imethod, h_neglect, h_neglect_edge )
1884  call build_reconstructions_1d( remap_cs, nk, hr, tr, poly_t_r, tir, poly_slope, imethod, h_neglect, h_neglect_edge )
1885  call mark_unstable_cells( nk, drdt, drds, til, sil, stable_l, ns_l )
1886  call mark_unstable_cells( nk, drdt, drds, tir, sir, stable_r, ns_r )
1887  call find_neutral_surface_positions_discontinuous(cs, nk, ns_l+ns_r, pres_l, hl, til, sil, drdt, drds, stable_l, &
1888  pres_r, hr, tir, sir, drdt, drds, stable_r, pol, por, kol, kor, heff)
1889  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, kol, kor, pol, por, heff, &
1890  (/1,1,2,2,3,3,3,3,3,3,3,3/), & ! KoL
1891  (/1,1,1,1,1,1,1,1,2,2,3,3/), & ! KoR
1892  (/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
1893  (/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
1894  (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/), & ! hEff
1895  'Right column much cooler')
1896  tl = (/14.,14.,10./) ; tr = (/14.,14.,10./)
1897  call build_reconstructions_1d( remap_cs, nk, hl, tl, poly_t_l, til, poly_slope, imethod, h_neglect, h_neglect_edge )
1898  call build_reconstructions_1d( remap_cs, nk, hr, tr, poly_t_r, tir, poly_slope, imethod, h_neglect, h_neglect_edge )
1899  call mark_unstable_cells( nk, drdt, drds, til, sil, stable_l, ns_l )
1900  call mark_unstable_cells( nk, drdt, drds, tir, sir, stable_r, ns_r )
1901  call find_neutral_surface_positions_discontinuous(cs, nk, ns_l+ns_r, pres_l, hl, til, sil, drdt, drds, stable_l, &
1902  pres_r, hr, tir, sir, drdt, drds, stable_r, pol, por, kol, kor, heff)
1903  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, kol, kor, pol, por, heff, &
1904  (/1,1,1,1,2,2,2,2,3,3,3,3/), & ! KoL
1905  (/1,1,1,1,2,2,2,2,3,3,3,3/), & ! KoR
1906  (/0.,0.,1.,1.,0.,0.,1.,1.,0.,0.,1.,1./), & ! pL
1907  (/0.,0.,1.,1.,0.,0.,1.,1.,0.,0.,1.,1./), & ! pR
1908  (/0.,10.,0.,0.,0.,10.,0.,0.,0.,10.,0.,0./), & ! hEff
1909  'Identical columns with mixed layer')
1910  tl = (/20.,16.,12./) ; tr = (/14.,14.,10./)
1911  call build_reconstructions_1d( remap_cs, nk, hl, tl, poly_t_l, til, poly_slope, imethod, h_neglect, h_neglect_edge )
1912  call build_reconstructions_1d( remap_cs, nk, hr, tr, poly_t_r, tir, poly_slope, imethod, h_neglect, h_neglect_edge )
1913  call mark_unstable_cells( nk, drdt, drds, til, sil, stable_l, ns_l )
1914  call mark_unstable_cells( nk, drdt, drds, tir, sir, stable_r, ns_r )
1915  call find_neutral_surface_positions_discontinuous(cs, nk, ns_l+ns_r, pres_l, hl, til, sil, drdt, drds, stable_l, &
1916  pres_r, hr, tir, sir, drdt, drds, stable_r, pol, por, kol, kor, heff)
1917  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, kol, kor, pol, por, heff, &
1918  (/1,1,2,2,2,3,3,3,3,3,3,3/), & ! KoL
1919  (/1,1,1,1,1,1,2,2,2,3,3,3/), & ! KoR
1920  (/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
1921  (/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
1922  (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 5.0, 0.0/), & ! hEff
1923  'Right column with mixed layer')
1924  tl = (/14.,14.,6./) ; tr = (/12.,16.,8./)
1925  call build_reconstructions_1d( remap_cs, nk, hl, tl, poly_t_l, til, poly_slope, imethod, h_neglect, h_neglect_edge )
1926  call build_reconstructions_1d( remap_cs, nk, hr, tr, poly_t_r, tir, poly_slope, imethod, h_neglect, h_neglect_edge )
1927  call mark_unstable_cells( nk, drdt, drds, til, sil, stable_l, ns_l )
1928  call mark_unstable_cells( nk, drdt, drds, tir, sir, stable_r, ns_r )
1929  call find_neutral_surface_positions_discontinuous(cs, nk, ns_l+ns_r, pres_l, hl, til, sil, drdt, drds, stable_l, &
1930  pres_r, hr, tir, sir, drdt, drds, stable_r, pol, por, kol, kor, heff)
1931  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 10, kol, kor, pol, por, heff, &
1932  (/1,1,1,1,2,2,2,3,3,3/), & ! KoL
1933  (/2,2,2,3,3,3,3,3,3,3/), & ! KoR
1934  (/0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, .75, 1.0/), & ! pL
1935  (/0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, .25, 1.0, 1.0/), & ! pR
1936  (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 7.5, 0.0/), & ! hEff
1937  'Left mixed layer, right unstable mixed layer')
1938 
1939  tl = (/10.,11.,6./) ; tr = (/12.,13.,8./)
1940  til(:,1) = (/8.,12.,10./) ; til(:,2) = (/12.,10.,2./)
1941  tir(:,1) = (/10.,14.,12./) ; tir(:,2) = (/14.,12.,4./)
1942  call mark_unstable_cells( nk, drdt, drds, til, sil, stable_l, ns_l )
1943  call mark_unstable_cells( nk, drdt, drds, tir, sir, stable_r, ns_r )
1944  call find_neutral_surface_positions_discontinuous(cs, nk, ns_l+ns_r, pres_l, hl, til, sil, drdt, drds, stable_l, &
1945  pres_r, hr, tir, sir, drdt, drds, stable_r, pol, por, kol, kor, heff)
1946  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 8, kol, kor, pol, por, heff, &
1947  (/2,2,2,2,2,3,3,3/), & ! KoL
1948  (/2,2,2,3,3,3,3,3/), & ! KoR
1949  (/0.0, 0.0, 0.0, 0.0, 1.0, 0.0, .75, 1.0/), & ! pL
1950  (/0.0, 1.0, 1.0, 0.0, .25, .25, 1.0, 1.0/), & ! pR
1951  (/0.0, 0.0, 0.0, 4.0, 0.0, 7.5, 0.0/), & ! hEff
1952  'Two unstable mixed layers')
1953  deallocate(remap_cs)
1954 
1955  allocate(eos)
1956  call eos_manual_init(eos, form_of_eos = eos_linear, drho_dt = -1., drho_ds = 2.)
1957  ! Unit tests for refine_nondim_position
1958  allocate(cs%ndiff_aux_CS)
1959  call set_ndiff_aux_params(cs%ndiff_aux_CS, deg = 1, max_iter = 10, drho_tol = 0., xtol = 0., eos = eos)
1960  ! Tests using Newton's method
1961  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5,refine_nondim_position( &
1962  cs%ndiff_aux_CS, 20., 35., -1., 2., 0., 1., (/21., -2./), (/35., 0./), -1., 1., 0.), &
1963  "Temperature stratified (Newton) "))
1964  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5,refine_nondim_position( &
1965  cs%ndiff_aux_CS, 20., 35., -1., 2., 0., 1., (/20., 0./), (/34., 2./), -2., 2., 0.), &
1966  "Salinity stratified (Newton) "))
1967  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5,refine_nondim_position( &
1968  cs%ndiff_aux_CS, 20., 35., -1., 2., 0., 1., (/21., -2./), (/34., 2./), -1., 1., 0.), &
1969  "Temp/Salt stratified (Newton) "))
1970  call set_ndiff_aux_params(cs%ndiff_aux_CS, force_brent = .true.)
1971  ! Tests using Brent's method
1972  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5,refine_nondim_position( &
1973  cs%ndiff_aux_CS, 20., 35., -1., 2., 0., 1., (/21., -2./), (/35., 0./), -1., 1., 0.), &
1974  "Temperature stratified (Brent) "))
1975  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5,refine_nondim_position( &
1976  cs%ndiff_aux_CS, 20., 35., -1., 2., 0., 1., (/20., 0./), (/34., 2./), -2., 2., 0.), &
1977  "Salinity stratified (Brent) "))
1978  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5,refine_nondim_position( &
1979  cs%ndiff_aux_CS, 20., 35., -1., 2., 0., 1., (/21., -2./), (/34., 2./), -1., 1., 0.), &
1980  "Temp/Salt stratified (Brent) "))
1981  deallocate(eos)
1982 
1983  if (.not. ndiff_unit_tests_discontinuous) write(*,*) 'Pass'
1984 

◆ neutral_diffusion()

subroutine, public mom_neutral_diffusion::neutral_diffusion ( type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h,
real, dimension( g %isdb: g %iedb, g %jsd: g %jed), intent(in)  Coef_x,
real, dimension( g %isd: g %ied, g %jsdb: g %jedb), intent(in)  Coef_y,
real, intent(in)  dt,
type(tracer_registry_type), pointer  Reg,
type(neutral_diffusion_cs), pointer  CS 
)

Update tracer concentration due to neutral diffusion; layer thickness unchanged by this update.

Parameters
[in]gOcean grid structure
[in]gvocean vertical grid structure
[in]hLayer thickness [H ~> m or kg m-2]
[in]coef_xdt * Kh * dy / dx at u-points [m2]
[in]coef_ydt * Kh * dx / dy at v-points [m2]
[in]dtTracer time step * I_numitts (I_numitts in tracer_hordiff)
regTracer registry
csNeutral diffusion control structure

Definition at line 411 of file MOM_neutral_diffusion.F90.

411  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
412  type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
413  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
414  real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: Coef_x !< dt * Kh * dy / dx at u-points [m2]
415  real, dimension(SZI_(G),SZJB_(G)), intent(in) :: Coef_y !< dt * Kh * dx / dy at v-points [m2]
416  real, intent(in) :: dt !< Tracer time step * I_numitts
417  !! (I_numitts in tracer_hordiff)
418  type(tracer_registry_type), pointer :: Reg !< Tracer registry
419  type(neutral_diffusion_CS), pointer :: CS !< Neutral diffusion control structure
420 
421  ! Local variables
422  real, dimension(SZIB_(G),SZJ_(G),CS%nsurf-1) :: uFlx ! Zonal flux of tracer [H conc ~> m conc or conc kg m-2]
423  real, dimension(SZI_(G),SZJB_(G),CS%nsurf-1) :: vFlx ! Meridional flux of tracer
424  ! [H conc ~> m conc or conc kg m-2]
425  real, dimension(SZI_(G),SZJ_(G),G%ke) :: tendency ! tendency array for diagn
426  real, dimension(SZI_(G),SZJ_(G)) :: tendency_2d ! depth integrated content tendency for diagn
427  real, dimension(SZIB_(G),SZJ_(G)) :: trans_x_2d ! depth integrated diffusive tracer x-transport diagn
428  real, dimension(SZI_(G),SZJB_(G)) :: trans_y_2d ! depth integrated diffusive tracer y-transport diagn
429  real, dimension(G%ke) :: dTracer ! change in tracer concentration due to ndiffusion
430 
431  type(tracer_type), pointer :: Tracer => null() ! Pointer to the current tracer
432 
433  integer :: i, j, k, m, ks, nk
434  real :: Idt
435  real :: h_neglect, h_neglect_edge
436 
437  !### Try replacing both of these with GV%H_subroundoff
438  h_neglect_edge = gv%m_to_H*1.0e-10
439  h_neglect = gv%m_to_H*1.0e-30
440 
441  nk = gv%ke
442 
443  do m = 1,reg%ntr ! Loop over tracer registry
444 
445  tracer => reg%Tr(m)
446 
447  ! for diagnostics
448  if (tracer%id_dfxy_conc > 0 .or. tracer%id_dfxy_cont > 0 .or. tracer%id_dfxy_cont_2d > 0 .or. &
449  tracer%id_dfx_2d > 0 .or. tracer%id_dfy_2d > 0) then
450  idt = 1.0/dt
451  tendency(:,:,:) = 0.0
452  endif
453 
454  uflx(:,:,:) = 0.
455  vflx(:,:,:) = 0.
456 
457  ! x-flux
458  do j = g%jsc,g%jec ; do i = g%isc-1,g%iec
459  if (g%mask2dCu(i,j)>0.) then
460  call neutral_surface_flux(nk, cs%nsurf, cs%deg, h(i,j,:), h(i+1,j,:), &
461  tracer%t(i,j,:), tracer%t(i+1,j,:), &
462  cs%uPoL(i,j,:), cs%uPoR(i,j,:), &
463  cs%uKoL(i,j,:), cs%uKoR(i,j,:), &
464  cs%uhEff(i,j,:), uflx(i,j,:), &
465  cs%continuous_reconstruction, h_neglect, cs%remap_CS, h_neglect_edge)
466  endif
467  enddo ; enddo
468 
469  ! y-flux
470  do j = g%jsc-1,g%jec ; do i = g%isc,g%iec
471  if (g%mask2dCv(i,j)>0.) then
472  call neutral_surface_flux(nk, cs%nsurf, cs%deg, h(i,j,:), h(i,j+1,:), &
473  tracer%t(i,j,:), tracer%t(i,j+1,:), &
474  cs%vPoL(i,j,:), cs%vPoR(i,j,:), &
475  cs%vKoL(i,j,:), cs%vKoR(i,j,:), &
476  cs%vhEff(i,j,:), vflx(i,j,:), &
477  cs%continuous_reconstruction, h_neglect, cs%remap_CS, h_neglect_edge)
478  endif
479  enddo ; enddo
480 
481  ! Update the tracer concentration from divergence of neutral diffusive flux components
482  do j = g%jsc,g%jec ; do i = g%isc,g%iec
483  if (g%mask2dT(i,j)>0.) then
484 
485  dtracer(:) = 0.
486  do ks = 1,cs%nsurf-1
487  k = cs%uKoL(i,j,ks)
488  dtracer(k) = dtracer(k) + coef_x(i,j) * uflx(i,j,ks)
489  k = cs%uKoR(i-1,j,ks)
490  dtracer(k) = dtracer(k) - coef_x(i-1,j) * uflx(i-1,j,ks)
491  k = cs%vKoL(i,j,ks)
492  dtracer(k) = dtracer(k) + coef_y(i,j) * vflx(i,j,ks)
493  k = cs%vKoR(i,j-1,ks)
494  dtracer(k) = dtracer(k) - coef_y(i,j-1) * vflx(i,j-1,ks)
495  enddo
496  do k = 1, gv%ke
497  tracer%t(i,j,k) = tracer%t(i,j,k) + dtracer(k) * &
498  ( g%IareaT(i,j) / ( h(i,j,k) + gv%H_subroundoff ) )
499  enddo
500 
501  if (tracer%id_dfxy_conc > 0 .or. tracer%id_dfxy_cont > 0 .or. tracer%id_dfxy_cont_2d > 0 ) then
502  do k = 1, gv%ke
503  tendency(i,j,k) = dtracer(k) * g%IareaT(i,j) * idt
504  enddo
505  endif
506 
507  endif
508  enddo ; enddo
509 
510  ! Diagnose vertically summed zonal flux, giving zonal tracer transport from ndiff.
511  ! Note sign corresponds to downgradient flux convention.
512  if (tracer%id_dfx_2d > 0) then
513  do j = g%jsc,g%jec ; do i = g%isc-1,g%iec
514  trans_x_2d(i,j) = 0.
515  if (g%mask2dCu(i,j)>0.) then
516  do ks = 1,cs%nsurf-1
517  trans_x_2d(i,j) = trans_x_2d(i,j) - coef_x(i,j) * uflx(i,j,ks)
518  enddo
519  trans_x_2d(i,j) = trans_x_2d(i,j) * idt
520  endif
521  enddo ; enddo
522  call post_data(tracer%id_dfx_2d, trans_x_2d(:,:), cs%diag)
523  endif
524 
525  ! Diagnose vertically summed merid flux, giving meridional tracer transport from ndiff.
526  ! Note sign corresponds to downgradient flux convention.
527  if (tracer%id_dfy_2d > 0) then
528  do j = g%jsc-1,g%jec ; do i = g%isc,g%iec
529  trans_y_2d(i,j) = 0.
530  if (g%mask2dCv(i,j)>0.) then
531  do ks = 1,cs%nsurf-1
532  trans_y_2d(i,j) = trans_y_2d(i,j) - coef_y(i,j) * vflx(i,j,ks)
533  enddo
534  trans_y_2d(i,j) = trans_y_2d(i,j) * idt
535  endif
536  enddo ; enddo
537  call post_data(tracer%id_dfy_2d, trans_y_2d(:,:), cs%diag)
538  endif
539 
540  ! post tendency of tracer content
541  if (tracer%id_dfxy_cont > 0) then
542  call post_data(tracer%id_dfxy_cont, tendency(:,:,:), cs%diag)
543  endif
544 
545  ! post depth summed tendency for tracer content
546  if (tracer%id_dfxy_cont_2d > 0) then
547  tendency_2d(:,:) = 0.
548  do j = g%jsc,g%jec ; do i = g%isc,g%iec
549  do k = 1, gv%ke
550  tendency_2d(i,j) = tendency_2d(i,j) + tendency(i,j,k)
551  enddo
552  enddo ; enddo
553  call post_data(tracer%id_dfxy_cont_2d, tendency_2d(:,:), cs%diag)
554  endif
555 
556  ! post tendency of tracer concentration; this step must be
557  ! done after posting tracer content tendency, since we alter
558  ! the tendency array.
559  if (tracer%id_dfxy_conc > 0) then
560  do k = 1, gv%ke ; do j = g%jsc,g%jec ; do i = g%isc,g%iec
561  tendency(i,j,k) = tendency(i,j,k) / ( h(i,j,k) + gv%H_subroundoff )
562  enddo ; enddo ; enddo
563  call post_data(tracer%id_dfxy_conc, tendency, cs%diag)
564  endif
565  enddo ! Loop over tracer registry
566 

◆ neutral_diffusion_calc_coeffs()

subroutine, public mom_neutral_diffusion::neutral_diffusion_calc_coeffs ( type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  T,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  S,
type(neutral_diffusion_cs), pointer  CS 
)

Calculate remapping factors for u/v columns used to map adjoining columns to a shared coordinate space.

Parameters
[in]gOcean grid structure
[in]gvocean vertical grid structure
[in]hLayer thickness [H ~> m or kg m-2]
[in]tPotential temperature [degC]
[in]sSalinity [ppt]
csNeutral diffusion control structure

Definition at line 234 of file MOM_neutral_diffusion.F90.

234  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
235  type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
236  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
237  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: T !< Potential temperature [degC]
238  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: S !< Salinity [ppt]
239  type(neutral_diffusion_CS), pointer :: CS !< Neutral diffusion control structure
240 
241  ! Local variables
242  integer :: i, j, k
243  ! Variables used for reconstructions
244  real, dimension(SZK_(G),2) :: ppoly_r_S ! Reconstruction slopes
245  real, dimension(SZI_(G), SZJ_(G)) :: hEff_sum
246  integer :: iMethod
247  real, dimension(SZI_(G)) :: ref_pres ! Reference pressure used to calculate alpha/beta
248  real :: h_neglect, h_neglect_edge
249  real :: pa_to_H
250 
251  pa_to_h = 1. / gv%H_to_pa
252 
253  !### Try replacing both of these with GV%H_subroundoff
254  if (gv%Boussinesq) then
255  h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
256  else
257  h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10
258  endif
259 
260  ! If doing along isopycnal diffusion (as opposed to neutral diffusion, set the reference pressure)
261  if (cs%ref_pres>=0.) then
262  ref_pres(:) = cs%ref_pres
263  endif
264 
265  if (cs%continuous_reconstruction) then
266  cs%dRdT(:,:,:) = 0.
267  cs%dRdS(:,:,:) = 0.
268  else
269  cs%T_i(:,:,:,:) = 0.
270  cs%S_i(:,:,:,:) = 0.
271  cs%dRdT_i(:,:,:,:) = 0.
272  cs%dRdS_i(:,:,:,:) = 0.
273  cs%ns(:,:) = 0.
274  cs%stable_cell(:,:,:) = .true.
275  endif
276 
277  ! Calculate pressure at interfaces and layer averaged alpha/beta
278  cs%Pint(:,:,1) = 0.
279  do k=1,g%ke ; do j=g%jsc-1, g%jec+1 ; do i=g%isc-1,g%iec+1
280  cs%Pint(i,j,k+1) = cs%Pint(i,j,k) + h(i,j,k)*gv%H_to_Pa
281  enddo ; enddo ; enddo
282 
283  do j = g%jsc-1, g%jec+1
284  ! Interpolate state to interface
285  do i = g%isc-1, g%iec+1
286  if (cs%continuous_reconstruction) then
287  call interface_scalar(g%ke, h(i,j,:), t(i,j,:), cs%Tint(i,j,:), 2, h_neglect)
288  call interface_scalar(g%ke, h(i,j,:), s(i,j,:), cs%Sint(i,j,:), 2, h_neglect)
289  else
290  call build_reconstructions_1d( cs%remap_CS, g%ke, h(i,j,:), t(i,j,:), cs%ppoly_coeffs_T(i,j,:,:), &
291  cs%T_i(i,j,:,:), ppoly_r_s, imethod, h_neglect, h_neglect_edge )
292  call build_reconstructions_1d( cs%remap_CS, g%ke, h(i,j,:), s(i,j,:), cs%ppoly_coeffs_S(i,j,:,:), &
293  cs%S_i(i,j,:,:), ppoly_r_s, imethod, h_neglect, h_neglect_edge )
294  endif
295  enddo
296 
297  ! Continuous reconstruction
298  if (cs%continuous_reconstruction) then
299  do k = 1, g%ke+1
300  if (cs%ref_pres<0) ref_pres(:) = cs%Pint(:,j,k)
301  call calculate_density_derivs(cs%Tint(:,j,k), cs%Sint(:,j,k), ref_pres, &
302  cs%dRdT(:,j,k), cs%dRdS(:,j,k), g%isc-1, g%iec-g%isc+3, cs%EOS)
303  enddo
304  else ! Discontinuous reconstruction
305  do k = 1, g%ke
306  if (cs%ref_pres<0) ref_pres(:) = cs%Pint(:,j,k)
307  call calculate_density_derivs(cs%T_i(:,j,k,1), cs%S_i(:,j,k,1), ref_pres, &
308  cs%dRdT_i(:,j,k,1), cs%dRdS_i(:,j,k,1), g%isc-1, g%iec-g%isc+3, cs%EOS)
309  if (cs%ref_pres<0) then
310  ref_pres(:) = cs%Pint(:,j,k+1)
311  endif
312  call calculate_density_derivs(cs%T_i(:,j,k,2), cs%S_i(:,j,k,2), ref_pres, &
313  cs%dRdT_i(:,j,k,2), cs%dRdS_i(:,j,k,2), g%isc-1, g%iec-g%isc+3, cs%EOS)
314  enddo
315  endif
316  enddo
317 
318  if (.not. cs%continuous_reconstruction) then
319  do j = g%jsc-1, g%jec+1 ; do i = g%isc-1, g%iec+1
320  call mark_unstable_cells( g%ke, cs%dRdT_i(i,j,:,:), cs%dRdS_i(i,j,:,:), cs%T_i(i,j,:,:), cs%S_i(i,j,:,:), &
321  cs%stable_cell(i,j,:), cs%ns(i,j) )
322  enddo ; enddo
323  endif
324 
325  cs%uhEff(:,:,:) = 0.
326  cs%vhEff(:,:,:) = 0.
327  cs%uPoL(:,:,:) = 0.
328  cs%vPoL(:,:,:) = 0.
329  cs%uPoR(:,:,:) = 0.
330  cs%vPoR(:,:,:) = 0.
331  cs%uKoL(:,:,:) = 1
332  cs%vKoL(:,:,:) = 1
333  cs%uKoR(:,:,:) = 1
334  cs%vKoR(:,:,:) = 1
335 
336  ! Neutral surface factors at U points
337  do j = g%jsc, g%jec ; do i = g%isc-1, g%iec
338  if (g%mask2dCu(i,j) > 0.) then
339  if (cs%continuous_reconstruction) then
340  call find_neutral_surface_positions_continuous(g%ke, &
341  cs%Pint(i,j,:), cs%Tint(i,j,:), cs%Sint(i,j,:), cs%dRdT(i,j,:), cs%dRdS(i,j,:), &
342  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,:), &
343  cs%uPoL(i,j,:), cs%uPoR(i,j,:), cs%uKoL(i,j,:), cs%uKoR(i,j,:), cs%uhEff(i,j,:) )
344  else
345  call find_neutral_surface_positions_discontinuous(cs, g%ke, cs%ns(i,j)+cs%ns(i+1,j), &
346  cs%Pint(i,j,:), h(i,j,:), cs%T_i(i,j,:,:), cs%S_i(i,j,:,:), &
347  cs%dRdT_i(i,j,:,:), cs%dRdS_i(i,j,:,:), cs%stable_cell(i,j,:), &
348  cs%Pint(i+1,j,:), h(i+1,j,:), cs%T_i(i+1,j,:,:), cs%S_i(i+1,j,:,:), &
349  cs%dRdT_i(i+1,j,:,:), cs%dRdS_i(i+1,j,:,:), cs%stable_cell(i+1,j,:), &
350  cs%uPoL(i,j,:), cs%uPoR(i,j,:), cs%uKoL(i,j,:), cs%uKoR(i,j,:), cs%uhEff(i,j,:), &
351  cs%ppoly_coeffs_T(i,j,:,:), cs%ppoly_coeffs_S(i,j,:,:), &
352  cs%ppoly_coeffs_T(i+1,j,:,:), cs%ppoly_coeffs_S(i+1,j,:,:))
353  endif
354  endif
355  enddo ; enddo
356 
357  ! Neutral surface factors at V points
358  do j = g%jsc-1, g%jec ; do i = g%isc, g%iec
359  if (g%mask2dCv(i,j) > 0.) then
360  if (cs%continuous_reconstruction) then
361  call find_neutral_surface_positions_continuous(g%ke, &
362  cs%Pint(i,j,:), cs%Tint(i,j,:), cs%Sint(i,j,:), cs%dRdT(i,j,:), cs%dRdS(i,j,:), &
363  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,:), &
364  cs%vPoL(i,j,:), cs%vPoR(i,j,:), cs%vKoL(i,j,:), cs%vKoR(i,j,:), cs%vhEff(i,j,:) )
365  else
366  call find_neutral_surface_positions_discontinuous(cs, g%ke, cs%ns(i,j)+cs%ns(i,j+1), &
367  cs%Pint(i,j,:), h(i,j,:), cs%T_i(i,j,:,:), cs%S_i(i,j,:,:), &
368  cs%dRdT_i(i,j,:,:), cs%dRdS_i(i,j,:,:), cs%stable_cell(i,j,:), &
369  cs%Pint(i,j+1,:), h(i,j+1,:), cs%T_i(i,j+1,:,:), cs%S_i(i,j+1,:,:), &
370  cs%dRdT_i(i,j+1,:,:), cs%dRdS_i(i,j+1,:,:), cs%stable_cell(i,j+1,:), &
371  cs%vPoL(i,j,:), cs%vPoR(i,j,:), cs%vKoL(i,j,:), cs%vKoR(i,j,:), cs%vhEff(i,j,:), &
372  cs%ppoly_coeffs_T(i,j,:,:), cs%ppoly_coeffs_S(i,j,:,:), &
373  cs%ppoly_coeffs_T(i,j+1,:,:), cs%ppoly_coeffs_S(i,j+1,:,:))
374 
375  endif
376  endif
377  enddo ; enddo
378 
379  ! Continuous reconstructions calculate hEff as the difference between the pressures of the
380  ! neutral surfaces which need to be reconverted to thickness units. The discontinuous version
381  ! calculates hEff from the fraction of the nondimensional fraction of the layer occupied by
382  ! the... (Please finish this thought. -RWH)
383  if (cs%continuous_reconstruction) then
384  do k = 1, cs%nsurf-1 ; do j = g%jsc, g%jec ; do i = g%isc-1, g%iec
385  if (g%mask2dCu(i,j) > 0.) cs%uhEff(i,j,k) = cs%uhEff(i,j,k) * pa_to_h
386  enddo ; enddo ; enddo
387  do k = 1, cs%nsurf-1 ; do j = g%jsc-1, g%jec ; do i = g%isc, g%iec
388  if (g%mask2dCv(i,j) > 0.) cs%vhEff(i,j,k) = cs%vhEff(i,j,k) * pa_to_h
389  enddo ; enddo ; enddo
390  endif
391 
392  if (cs%id_uhEff_2d>0) then
393  heff_sum(:,:) = 0.
394  do k = 1,cs%nsurf-1 ; do j=g%jsc,g%jec ; do i=g%isc-1,g%iec
395  heff_sum(i,j) = heff_sum(i,j) + cs%uhEff(i,j,k)
396  enddo ; enddo ; enddo
397  call post_data(cs%id_uhEff_2d, heff_sum, cs%diag)
398  endif
399  if (cs%id_vhEff_2d>0) then
400  heff_sum(:,:) = 0.
401  do k = 1,cs%nsurf-1 ; do j=g%jsc-1,g%jec ; do i=g%isc,g%iec
402  heff_sum(i,j) = heff_sum(i,j) + cs%vhEff(i,j,k)
403  enddo ; enddo ; enddo
404  call post_data(cs%id_vhEff_2d, heff_sum, cs%diag)
405  endif
406 

◆ neutral_diffusion_end()

subroutine, public mom_neutral_diffusion::neutral_diffusion_end ( type(neutral_diffusion_cs), pointer  CS)

Deallocates neutral_diffusion control structure.

Parameters
csNeutral diffusion control structure

Definition at line 2240 of file MOM_neutral_diffusion.F90.

2240  type(neutral_diffusion_CS), pointer :: CS !< Neutral diffusion control structure
2241 
2242  if (associated(cs)) deallocate(cs)
2243 

◆ neutral_diffusion_init()

logical function, public mom_neutral_diffusion::neutral_diffusion_init ( type(time_type), intent(in), target  Time,
type(ocean_grid_type), intent(in)  G,
type(param_file_type), intent(in)  param_file,
type(diag_ctrl), intent(inout), target  diag,
type(eos_type), intent(in), target  EOS,
type(neutral_diffusion_cs), pointer  CS 
)

Read parameters and allocate control structure for neutral_diffusion module.

Parameters
[in]timeTime structure
[in]gGrid structure
[in,out]diagDiagnostics control structure
[in]param_fileParameter file structure
[in]eosEquation of state
csNeutral diffusion control structure

Definition at line 101 of file MOM_neutral_diffusion.F90.

101  type(time_type), target, intent(in) :: Time !< Time structure
102  type(ocean_grid_type), intent(in) :: G !< Grid structure
103  type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure
104  type(param_file_type), intent(in) :: param_file !< Parameter file structure
105  type(EOS_type), target, intent(in) :: EOS !< Equation of state
106  type(neutral_diffusion_CS), pointer :: CS !< Neutral diffusion control structure
107 
108  ! Local variables
109  character(len=256) :: mesg ! Message for error messages.
110  character(len=80) :: string ! Temporary strings
111  logical :: boundary_extrap
112  ! For refine_pos
113  integer :: max_iter
114  real :: drho_tol, xtol, ref_pres
115 
116  if (associated(cs)) then
117  call mom_error(fatal, "neutral_diffusion_init called with associated control structure.")
118  return
119  endif
120 
121  ! Log this module and master switch for turning it on/off
122  call log_version(param_file, mdl, version, &
123  "This module implements neutral diffusion of tracers")
124  call get_param(param_file, mdl, "USE_NEUTRAL_DIFFUSION", neutral_diffusion_init, &
125  "If true, enables the neutral diffusion module.", &
126  default=.false.)
127 
128  if (.not.neutral_diffusion_init) then
129  return
130  endif
131 
132  allocate(cs)
133  allocate(cs%ndiff_aux_CS)
134  cs%diag => diag
135  cs%EOS => eos
136  ! call openParameterBlock(param_file,'NEUTRAL_DIFF')
137 
138  ! Read all relevant parameters and write them to the model log.
139  call get_param(param_file, mdl, "NDIFF_CONTINUOUS", cs%continuous_reconstruction, &
140  "If true, uses a continuous reconstruction of T and S when "//&
141  "finding neutral surfaces along which diffusion will happen. "//&
142  "If false, a PPM discontinuous reconstruction of T and S "//&
143  "is done which results in a higher order routine but exacts "//&
144  "a higher computational cost.", default=.true.)
145  call get_param(param_file, mdl, "NDIFF_REF_PRES", cs%ref_pres, &
146  "The reference pressure (Pa) used for the derivatives of "//&
147  "the equation of state. If negative (default), local "//&
148  "pressure is used.", &
149  default = -1.)
150  ! Initialize and configure remapping
151  if (cs%continuous_reconstruction .eqv. .false.) then
152  call get_param(param_file, mdl, "NDIFF_BOUNDARY_EXTRAP", boundary_extrap, &
153  "Uses a rootfinding approach to find the position of a "//&
154  "neutral surface within a layer taking into account the "//&
155  "nonlinearity of the equation of state and the "//&
156  "polynomial reconstructions of T/S.", &
157  default=.false.)
158  call get_param(param_file, mdl, "NDIFF_REMAPPING_SCHEME", string, &
159  "This sets the reconstruction scheme used "//&
160  "for vertical remapping for all variables. "//&
161  "It can be one of the following schemes: "//&
162  trim(remappingschemesdoc), default=remappingdefaultscheme)
163  call initialize_remapping( cs%remap_CS, string, boundary_extrapolation = boundary_extrap )
164  call extract_member_remapping_cs(cs%remap_CS, degree=cs%deg)
165  call get_param(param_file, mdl, "NDIFF_REFINE_POSITION", cs%refine_position, &
166  "Uses a rootfinding approach to find the position of a "//&
167  "neutral surface within a layer taking into account the "//&
168  "nonlinearity of the equation of state and the "//&
169  "polynomial reconstructions of T/S.", &
170  default=.false.)
171  if (cs%refine_position) then
172  call get_param(param_file, mdl, "NDIFF_DRHO_TOL", drho_tol, &
173  "Sets the convergence criterion for finding the neutral "//&
174  "position within a layer in kg m-3.", &
175  default=1.e-10)
176  call get_param(param_file, mdl, "NDIFF_X_TOL", xtol, &
177  "Sets the convergence criterion for a change in nondim "//&
178  "position within a layer.", &
179  default=0.)
180  call get_param(param_file, mdl, "NDIFF_MAX_ITER", max_iter, &
181  "The maximum number of iterations to be done before "//&
182  "exiting the iterative loop to find the neutral surface", &
183  default=10)
184  call set_ndiff_aux_params(cs%ndiff_aux_CS, max_iter = max_iter, drho_tol = drho_tol, xtol = xtol)
185  endif
186  call get_param(param_file, mdl, "NDIFF_DEBUG", cs%debug, &
187  "Turns on verbose output for discontinuous neutral "//&
188  "diffusion routines.", &
189  default = .false.)
190  call set_ndiff_aux_params(cs%ndiff_aux_CS, deg=cs%deg, ref_pres = cs%ref_pres, eos = eos, debug = cs%debug)
191  endif
192 
193 ! call get_param(param_file, mdl, "KHTR", CS%KhTr, &
194 ! "The background along-isopycnal tracer diffusivity.", &
195 ! units="m2 s-1", default=0.0)
196 ! call closeParameterBlock(param_file)
197  if (cs%continuous_reconstruction) then
198  cs%nsurf = 2*g%ke+2 ! Continuous reconstruction means that every interface has two connections
199  allocate(cs%dRdT(szi_(g),szj_(g),szk_(g)+1)) ; cs%dRdT(:,:,:) = 0.
200  allocate(cs%dRdS(szi_(g),szj_(g),szk_(g)+1)) ; cs%dRdS(:,:,:) = 0.
201  else
202  cs%nsurf = 4*g%ke ! Discontinuous means that every interface has four connections
203  allocate(cs%T_i(szi_(g),szj_(g),szk_(g),2)) ; cs%T_i(:,:,:,:) = 0.
204  allocate(cs%S_i(szi_(g),szj_(g),szk_(g),2)) ; cs%S_i(:,:,:,:) = 0.
205  allocate(cs%dRdT_i(szi_(g),szj_(g),szk_(g),2)) ; cs%dRdT_i(:,:,:,:) = 0.
206  allocate(cs%dRdS_i(szi_(g),szj_(g),szk_(g),2)) ; cs%dRdS_i(:,:,:,:) = 0.
207  allocate(cs%ppoly_coeffs_T(szi_(g),szj_(g),szk_(g),cs%deg+1)) ; cs%ppoly_coeffs_T(:,:,:,:) = 0.
208  allocate(cs%ppoly_coeffs_S(szi_(g),szj_(g),szk_(g),cs%deg+1)) ; cs%ppoly_coeffs_S(:,:,:,:) = 0.
209  allocate(cs%ns(szi_(g),szj_(g))) ; cs%ns(:,:) = 0.
210  endif
211  ! T-points
212  allocate(cs%Tint(szi_(g),szj_(g),szk_(g)+1)) ; cs%Tint(:,:,:) = 0.
213  allocate(cs%Sint(szi_(g),szj_(g),szk_(g)+1)) ; cs%Sint(:,:,:) = 0.
214  allocate(cs%Pint(szi_(g),szj_(g),szk_(g)+1)) ; cs%Pint(:,:,:) = 0.
215  allocate(cs%stable_cell(szi_(g),szj_(g),szk_(g))) ; cs%stable_cell(:,:,:) = .true.
216  ! U-points
217  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.
218  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.
219  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
220  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
221  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
222  ! V-points
223  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.
224  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.
225  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
226  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
227  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
228 

◆ neutral_diffusion_unit_tests()

logical function, public mom_neutral_diffusion::neutral_diffusion_unit_tests ( logical, intent(in)  verbose)

Returns true if unit tests of neutral_diffusion functions fail. Otherwise returns false.

Parameters
[in]verboseIf true, write results to stdout

Definition at line 1507 of file MOM_neutral_diffusion.F90.

1507  logical, intent(in) :: verbose !< If true, write results to stdout
1508 
1509  neutral_diffusion_unit_tests = .false. .or. &
1510  ndiff_unit_tests_continuous(verbose) .or. ndiff_unit_tests_discontinuous(verbose)
1511 
1512 

◆ neutral_surface_flux()

subroutine mom_neutral_diffusion::neutral_surface_flux ( integer, intent(in)  nk,
integer, intent(in)  nsurf,
integer, intent(in)  deg,
real, dimension(nk), intent(in)  hl,
real, dimension(nk), intent(in)  hr,
real, dimension(nk), intent(in)  Tl,
real, dimension(nk), intent(in)  Tr,
real, dimension(nsurf), intent(in)  PiL,
real, dimension(nsurf), intent(in)  PiR,
integer, dimension(nsurf), intent(in)  KoL,
integer, dimension(nsurf), intent(in)  KoR,
real, dimension(nsurf-1), intent(in)  hEff,
real, dimension(nsurf-1), intent(inout)  Flx,
logical, intent(in)  continuous,
real, intent(in)  h_neglect,
type(remapping_cs), intent(in), optional  remap_CS,
real, intent(in), optional  h_neglect_edge 
)
private

Returns a single column of neutral diffusion fluxes of a tracer.

Parameters
[in]nkNumber of levels
[in]nsurfNumber of neutral surfaces
[in]degDegree of polynomial reconstructions
[in]hlLeft-column layer thickness [Pa]
[in]hrRight-column layer thickness [Pa]
[in]tlLeft-column layer tracer (conc, e.g. degC)
[in]trRight-column layer tracer (conc, e.g. degC)
[in]pilFractional position of neutral surface within layer KoL of left column
[in]pirFractional position of neutral surface within layer KoR of right column
[in]kolIndex of first left interface above neutral surface
[in]korIndex of first right interface above neutral surface
[in]heffEffective thickness between two neutral surfaces [Pa]
[in,out]flxFlux of tracer between pairs of neutral layers (conc H)
[in]continuousTrue if using continuous reconstruction
[in]h_neglectA negligibly small width for the purpose of cell reconstructions in the same units as h0.
[in]remap_csRemapping control structure used to create sublayers
[in]h_neglect_edgeA negligibly small width for the purpose of edge value calculations in the same units as h0.

Definition at line 1303 of file MOM_neutral_diffusion.F90.

1303  integer, intent(in) :: nk !< Number of levels
1304  integer, intent(in) :: nsurf !< Number of neutral surfaces
1305  integer, intent(in) :: deg !< Degree of polynomial reconstructions
1306  real, dimension(nk), intent(in) :: hl !< Left-column layer thickness [Pa]
1307  real, dimension(nk), intent(in) :: hr !< Right-column layer thickness [Pa]
1308  real, dimension(nk), intent(in) :: Tl !< Left-column layer tracer (conc, e.g. degC)
1309  real, dimension(nk), intent(in) :: Tr !< Right-column layer tracer (conc, e.g. degC)
1310  real, dimension(nsurf), intent(in) :: PiL !< Fractional position of neutral surface
1311  !! within layer KoL of left column
1312  real, dimension(nsurf), intent(in) :: PiR !< Fractional position of neutral surface
1313  !! within layer KoR of right column
1314  integer, dimension(nsurf), intent(in) :: KoL !< Index of first left interface above neutral surface
1315  integer, dimension(nsurf), intent(in) :: KoR !< Index of first right interface above neutral surface
1316  real, dimension(nsurf-1), intent(in) :: hEff !< Effective thickness between two neutral surfaces [Pa]
1317  real, dimension(nsurf-1), intent(inout) :: Flx !< Flux of tracer between pairs of neutral layers (conc H)
1318  logical, intent(in) :: continuous !< True if using continuous reconstruction
1319  real, intent(in) :: h_neglect !< A negligibly small width for the
1320  !! purpose of cell reconstructions
1321  !! in the same units as h0.
1322  type(remapping_CS), optional, intent(in) :: remap_CS !< Remapping control structure used
1323  !! to create sublayers
1324  real, optional, intent(in) :: h_neglect_edge !< A negligibly small width
1325  !! for the purpose of edge value calculations
1326  !! in the same units as h0.
1327  ! Local variables
1328  integer :: k_sublayer, klb, klt, krb, krt, k
1329  real :: T_right_top, T_right_bottom, T_right_layer, T_right_sub, T_right_top_int, T_right_bot_int
1330  real :: T_left_top, T_left_bottom, T_left_layer, T_left_sub, T_left_top_int, T_left_bot_int
1331  real :: dT_top, dT_bottom, dT_layer, dT_ave, dT_sublayer, dT_top_int, dT_bot_int
1332  real, dimension(nk+1) :: Til !< Left-column interface tracer (conc, e.g. degC)
1333  real, dimension(nk+1) :: Tir !< Right-column interface tracer (conc, e.g. degC)
1334  real, dimension(nk) :: aL_l !< Left-column left edge value of tracer (conc, e.g. degC)
1335  real, dimension(nk) :: aR_l !< Left-column right edge value of tracer (conc, e.g. degC)
1336  real, dimension(nk) :: aL_r !< Right-column left edge value of tracer (conc, e.g. degC)
1337  real, dimension(nk) :: aR_r !< Right-column right edge value of tracer (conc, e.g. degC)
1338  ! Discontinuous reconstruction
1339  integer :: iMethod
1340  real, dimension(nk,2) :: Tid_l !< Left-column interface tracer (conc, e.g. degC)
1341  real, dimension(nk,2) :: Tid_r !< Right-column interface tracer (conc, e.g. degC)
1342  real, dimension(nk,deg+1) :: ppoly_r_coeffs_l
1343  real, dimension(nk,deg+1) :: ppoly_r_coeffs_r
1344  real, dimension(nk,deg+1) :: ppoly_r_S_l
1345  real, dimension(nk,deg+1) :: ppoly_r_S_r
1346  logical :: down_flux
1347  ! Setup reconstruction edge values
1348  if (continuous) then
1349  call interface_scalar(nk, hl, tl, til, 2, h_neglect)
1350  call interface_scalar(nk, hr, tr, tir, 2, h_neglect)
1351  call ppm_left_right_edge_values(nk, tl, til, al_l, ar_l)
1352  call ppm_left_right_edge_values(nk, tr, tir, al_r, ar_r)
1353  else
1354  ppoly_r_coeffs_l(:,:) = 0.
1355  ppoly_r_coeffs_r(:,:) = 0.
1356  tid_l(:,:) = 0.
1357  tid_r(:,:) = 0.
1358 
1359  call build_reconstructions_1d( remap_cs, nk, hl, tl, ppoly_r_coeffs_l, tid_l, &
1360  ppoly_r_s_l, imethod, h_neglect, h_neglect_edge )
1361  call build_reconstructions_1d( remap_cs, nk, hr, tr, ppoly_r_coeffs_r, tid_r, &
1362  ppoly_r_s_r, imethod, h_neglect, h_neglect_edge )
1363  endif
1364 
1365  do k_sublayer = 1, nsurf-1
1366  if (heff(k_sublayer) == 0.) then
1367  flx(k_sublayer) = 0.
1368  else
1369  if (continuous) then
1370  klb = kol(k_sublayer+1)
1371  t_left_bottom = ( 1. - pil(k_sublayer+1) ) * til(klb) + pil(k_sublayer+1) * til(klb+1)
1372  klt = kol(k_sublayer)
1373  t_left_top = ( 1. - pil(k_sublayer) ) * til(klt) + pil(k_sublayer) * til(klt+1)
1374  t_left_layer = ppm_ave(pil(k_sublayer), pil(k_sublayer+1) + real(klb-klt), &
1375  al_l(klt), ar_l(klt), tl(klt))
1376 
1377  krb = kor(k_sublayer+1)
1378  t_right_bottom = ( 1. - pir(k_sublayer+1) ) * tir(krb) + pir(k_sublayer+1) * tir(krb+1)
1379  krt = kor(k_sublayer)
1380  t_right_top = ( 1. - pir(k_sublayer) ) * tir(krt) + pir(k_sublayer) * tir(krt+1)
1381  t_right_layer = ppm_ave(pir(k_sublayer), pir(k_sublayer+1) + real(krb-krt), &
1382  al_r(krt), ar_r(krt), tr(krt))
1383  dt_top = t_right_top - t_left_top
1384  dt_bottom = t_right_bottom - t_left_bottom
1385  dt_ave = 0.5 * ( dt_top + dt_bottom )
1386  dt_layer = t_right_layer - t_left_layer
1387  if (signum(1.,dt_top) * signum(1.,dt_bottom) <= 0. .or. signum(1.,dt_ave) * signum(1.,dt_layer) <= 0.) then
1388  dt_ave = 0.
1389  else
1390  dt_ave = dt_layer
1391  endif
1392  flx(k_sublayer) = dt_ave * heff(k_sublayer)
1393  else ! Discontinuous reconstruction
1394  ! Calculate tracer values on left and right side of the neutral surface
1395  call neutral_surface_t_eval(nk, nsurf, k_sublayer, kol, pil, tl, tid_l, deg, imethod, &
1396  ppoly_r_coeffs_l, t_left_top, t_left_bottom, t_left_sub, &
1397  t_left_top_int, t_left_bot_int, t_left_layer)
1398  call neutral_surface_t_eval(nk, nsurf, k_sublayer, kor, pir, tr, tid_r, deg, imethod, &
1399  ppoly_r_coeffs_r, t_right_top, t_right_bottom, t_right_sub, &
1400  t_right_top_int, t_right_bot_int, t_right_layer)
1401 
1402  dt_top = t_right_top - t_left_top
1403  dt_bottom = t_right_bottom - t_left_bottom
1404  dt_sublayer = t_right_sub - t_left_sub
1405  dt_top_int = t_right_top_int - t_left_top_int
1406  dt_bot_int = t_right_bot_int - t_left_bot_int
1407  ! Enforcing the below criterion incorrectly zero out fluxes
1408  !dT_layer = T_right_layer - T_left_layer
1409 
1410  down_flux = dt_top <= 0. .and. dt_bottom <= 0. .and. &
1411  dt_sublayer <= 0. .and. dt_top_int <= 0. .and. &
1412  dt_bot_int <= 0.
1413  down_flux = down_flux .or. &
1414  (dt_top >= 0. .and. dt_bottom >= 0. .and. &
1415  dt_sublayer >= 0. .and. dt_top_int >= 0. .and. &
1416  dt_bot_int >= 0.)
1417  if (down_flux) then
1418  flx(k_sublayer) = dt_sublayer * heff(k_sublayer)
1419  else
1420  flx(k_sublayer) = 0.
1421  endif
1422  endif
1423  endif
1424  enddo
1425 

◆ neutral_surface_t_eval()

subroutine mom_neutral_diffusion::neutral_surface_t_eval ( integer, intent(in)  nk,
integer, intent(in)  ns,
integer, intent(in)  k_sub,
integer, dimension(ns), intent(in)  Ks,
real, dimension(ns), intent(in)  Ps,
real, dimension(nk), intent(in)  T_mean,
real, dimension(nk,2), intent(in)  T_int,
integer, intent(in)  deg,
integer, intent(in)  iMethod,
real, dimension(nk,deg+1), intent(in)  T_poly,
real, intent(out)  T_top,
real, intent(out)  T_bot,
real, intent(out)  T_sub,
real, intent(out)  T_top_int,
real, intent(out)  T_bot_int,
real, intent(out)  T_layer 
)
private

Evaluate various parts of the reconstructions to calculate gradient-based flux limter.

Parameters
[in]nkNumber of cell everages
[in]nsNumber of neutral surfaces
[in]k_subIndex of current neutral layer
[in]ksList of the layers associated with each neutral surface
[in]psList of the positions within a layer of each surface
[in]t_meanCell average of tracer
[in]t_intCell interface values of tracer from reconstruction
[in]degDegree of reconstruction polynomial (e.g. 1 is linear)
[in]imethodMethod of integration to use
[in]t_polyCoefficients of polynomial reconstructions
[out]t_topTracer value at top (across discontinuity if necessary)
[out]t_botTracer value at bottom (across discontinuity if necessary)
[out]t_subAverage of the tracer value over the sublayer
[out]t_top_intTracer value at top interface of neutral layer
[out]t_bot_intTracer value at bottom interface of neutral layer
[out]t_layerCell-average that the the reconstruction belongs to

Definition at line 1431 of file MOM_neutral_diffusion.F90.

1431  integer, intent(in ) :: nk !< Number of cell everages
1432  integer, intent(in ) :: ns !< Number of neutral surfaces
1433  integer, intent(in ) :: k_sub !< Index of current neutral layer
1434  integer, dimension(ns), intent(in ) :: Ks !< List of the layers associated with each neutral surface
1435  real, dimension(ns), intent(in ) :: Ps !< List of the positions within a layer of each surface
1436  real, dimension(nk), intent(in ) :: T_mean !< Cell average of tracer
1437  real, dimension(nk,2), intent(in ) :: T_int !< Cell interface values of tracer from reconstruction
1438  integer, intent(in ) :: deg !< Degree of reconstruction polynomial (e.g. 1 is linear)
1439  integer, intent(in ) :: iMethod !< Method of integration to use
1440  real, dimension(nk,deg+1), intent(in ) :: T_poly !< Coefficients of polynomial reconstructions
1441  real, intent( out) :: T_top !< Tracer value at top (across discontinuity if necessary)
1442  real, intent( out) :: T_bot !< Tracer value at bottom (across discontinuity if necessary)
1443  real, intent( out) :: T_sub !< Average of the tracer value over the sublayer
1444  real, intent( out) :: T_top_int !< Tracer value at top interface of neutral layer
1445  real, intent( out) :: T_bot_int !< Tracer value at bottom interface of neutral layer
1446  real, intent( out) :: T_layer !< Cell-average that the the reconstruction belongs to
1447 
1448  integer :: kl, ks_top, ks_bot
1449 
1450  ks_top = k_sub
1451  ks_bot = k_sub + 1
1452  if ( ks(ks_top) /= ks(ks_bot) ) then
1453  call mom_error(fatal, "Neutral surfaces span more than one layer")
1454  endif
1455  kl = ks(k_sub)
1456  ! First if the neutral surfaces spans the entirety of a cell, then do not search across the discontinuity
1457  if ( (ps(ks_top) == 0.) .and. (ps(ks_bot) == 1.)) then
1458  t_top = t_int(kl,1)
1459  t_bot = t_int(kl,2)
1460  else
1461  ! Search across potential discontinuity at top
1462  if ( (kl > 1) .and. (ps(ks_top) == 0.) ) then
1463  t_top = t_int(kl-1,2)
1464  else
1465  t_top = evaluation_polynomial( t_poly(kl,:), deg+1, ps(ks_top) )
1466  endif
1467  ! Search across potential discontinuity at bottom
1468  if ( (kl < nk) .and. (ps(ks_bot) == 1.) ) then
1469  t_bot = t_int(kl+1,1)
1470  else
1471  t_bot = evaluation_polynomial( t_poly(kl,:), deg+1, ps(ks_bot) )
1472  endif
1473  endif
1474  t_sub = average_value_ppoly(nk, t_mean, t_int, t_poly, imethod, kl, ps(ks_top), ps(ks_bot))
1475  t_top_int = evaluation_polynomial( t_poly(kl,:), deg+1, ps(ks_top))
1476  t_bot_int = evaluation_polynomial( t_poly(kl,:), deg+1, ps(ks_bot))
1477  t_layer = t_mean(kl)
1478 

◆ plm_diff()

subroutine mom_neutral_diffusion::plm_diff ( integer, intent(in)  nk,
real, dimension(nk), intent(in)  h,
real, dimension(nk), intent(in)  S,
integer, intent(in)  c_method,
integer, intent(in)  b_method,
real, dimension(nk), intent(inout)  diff 
)
private

Returns PLM slopes for a column where the slopes are the difference in value across each cell. The limiting follows equation 1.8 in Colella & Woodward, 1984: JCP 54, 174-201.

Parameters
[in]nkNumber of levels
[in]hLayer thickness [H ~> m or kg m-2]
[in]sLayer salinity (conc, e.g. ppt)
[in]c_methodMethod to use for the centered difference
[in]b_method=1, use PCM in first/last cell, =2 uses linear extrapolation
[in,out]diffScalar difference across layer (conc, e.g. ppt) determined by the following values for c_method:
  1. Second order finite difference (not recommended)
  2. Second order finite volume (used in original PPM)
  3. Finite-volume weighted least squares linear fit
Todo:
The use of c_method to choose a scheme is inefficient and should eventually be moved up the call tree.

Definition at line 687 of file MOM_neutral_diffusion.F90.

687  integer, intent(in) :: nk !< Number of levels
688  real, dimension(nk), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
689  real, dimension(nk), intent(in) :: S !< Layer salinity (conc, e.g. ppt)
690  integer, intent(in) :: c_method !< Method to use for the centered difference
691  integer, intent(in) :: b_method !< =1, use PCM in first/last cell, =2 uses linear extrapolation
692  real, dimension(nk), intent(inout) :: diff !< Scalar difference across layer (conc, e.g. ppt)
693  !! determined by the following values for c_method:
694  !! 1. Second order finite difference (not recommended)
695  !! 2. Second order finite volume (used in original PPM)
696  !! 3. Finite-volume weighted least squares linear fit
697  !! \todo The use of c_method to choose a scheme is inefficient
698  !! and should eventually be moved up the call tree.
699 
700  ! Local variables
701  integer :: k
702  real :: hkm1, hk, hkp1, Skm1, Sk, Skp1, diff_l, diff_r, diff_c
703 
704  do k = 2, nk-1
705  hkm1 = h(k-1)
706  hk = h(k)
707  hkp1 = h(k+1)
708 
709  if ( ( hkp1 + hk ) * ( hkm1 + hk ) > 0.) then
710  skm1 = s(k-1)
711  sk = s(k)
712  skp1 = s(k+1)
713  if (c_method==1) then
714  ! Simple centered diff (from White)
715  if ( hk + 0.5 * (hkm1 + hkp1) /= 0. ) then
716  diff_c = ( skp1 - skm1 ) * ( hk / ( hk + 0.5 * (hkm1 + hkp1) ) )
717  else
718  diff_c = 0.
719  endif
720  elseif (c_method==2) then
721  ! Second order accurate centered FV slope (from Colella and Woodward, JCP 1984)
722  diff_c = fv_diff(hkm1, hk, hkp1, skm1, sk, skp1)
723  elseif (c_method==3) then
724  ! Second order accurate finite-volume least squares slope
725  diff_c = hk * fvlsq_slope(hkm1, hk, hkp1, skm1, sk, skp1)
726  endif
727  ! Limit centered slope by twice the side differenced slopes
728  diff_l = 2. * ( sk - skm1 )
729  diff_r = 2. * ( skp1 - sk )
730  if ( signum(1., diff_l) * signum(1., diff_r) <= 0. ) then
731  diff(k) = 0. ! PCM for local extrema
732  else
733  diff(k) = sign( min( abs(diff_l), abs(diff_c), abs(diff_r) ), diff_c )
734  endif
735  else
736  diff(k) = 0. ! PCM next to vanished layers
737  endif
738  enddo
739  if (b_method==1) then ! PCM for top and bottom layer
740  diff(1) = 0.
741  diff(nk) = 0.
742  elseif (b_method==2) then ! Linear extrapolation for top and bottom interfaces
743  diff(1) = ( s(2) - s(1) ) * 2. * ( h(1) / ( h(1) + h(2) ) )
744  diff(nk) = s(nk) - s(nk-1) * 2. * ( h(nk) / ( h(nk-1) + h(nk) ) )
745  endif
746 

◆ ppm_ave()

real function mom_neutral_diffusion::ppm_ave ( real, intent(in)  xL,
real, intent(in)  xR,
real, intent(in)  aL,
real, intent(in)  aR,
real, intent(in)  aMean 
)
private

Returns the average of a PPM reconstruction between two fractional positions.

Parameters
[in]xlFraction position of left bound (0,1)
[in]xrFraction position of right bound (0,1)
[in]alLeft edge scalar value, at x=0
[in]arRight edge scalar value, at x=1
[in]ameanAverage scalar value of cell

Definition at line 649 of file MOM_neutral_diffusion.F90.

649  real, intent(in) :: xL !< Fraction position of left bound (0,1)
650  real, intent(in) :: xR !< Fraction position of right bound (0,1)
651  real, intent(in) :: aL !< Left edge scalar value, at x=0
652  real, intent(in) :: aR !< Right edge scalar value, at x=1
653  real, intent(in) :: aMean !< Average scalar value of cell
654 
655  ! Local variables
656  real :: dx, xave, a6, a6o3
657 
658  dx = xr - xl
659  xave = 0.5 * ( xr + xl )
660  a6o3 = 2. * amean - ( al + ar ) ! a6 / 3.
661  a6 = 3. * a6o3
662 
663  if (dx<0.) then
664  stop 'ppm_ave: dx<0 should not happend!'
665  elseif (dx>1.) then
666  stop 'ppm_ave: dx>1 should not happend!'
667  elseif (dx==0.) then
668  ppm_ave = al + ( ar - al ) * xr + a6 * xr * ( 1. - xr )
669  else
670  ppm_ave = ( al + xave * ( ( ar - al ) + a6 ) ) - a6o3 * ( xr**2 + xr * xl + xl**2 )
671  endif

◆ ppm_edge()

real function mom_neutral_diffusion::ppm_edge ( real, intent(in)  hkm1,
real, intent(in)  hk,
real, intent(in)  hkp1,
real, intent(in)  hkp2,
real, intent(in)  Ak,
real, intent(in)  Akp1,
real, intent(in)  Pk,
real, intent(in)  Pkp1,
real, intent(in)  h_neglect 
)
private

Returns the PPM quasi-fourth order edge value at k+1/2 following equation 1.6 in Colella & Woodward, 1984: JCP 54, 174-201.

Parameters
[in]hkm1Width of cell k-1
[in]hkWidth of cell k
[in]hkp1Width of cell k+1
[in]hkp2Width of cell k+2
[in]akAverage scalar value of cell k
[in]akp1Average scalar value of cell k+1
[in]pkPLM slope for cell k
[in]pkp1PLM slope for cell k+1
[in]h_neglectA negligibly small thickness [H ~> m or kg m-2]

Definition at line 609 of file MOM_neutral_diffusion.F90.

609  real, intent(in) :: hkm1 !< Width of cell k-1
610  real, intent(in) :: hk !< Width of cell k
611  real, intent(in) :: hkp1 !< Width of cell k+1
612  real, intent(in) :: hkp2 !< Width of cell k+2
613  real, intent(in) :: Ak !< Average scalar value of cell k
614  real, intent(in) :: Akp1 !< Average scalar value of cell k+1
615  real, intent(in) :: Pk !< PLM slope for cell k
616  real, intent(in) :: Pkp1 !< PLM slope for cell k+1
617  real, intent(in) :: h_neglect !< A negligibly small thickness [H ~> m or kg m-2]
618 
619  ! Local variables
620  real :: R_hk_hkp1, R_2hk_hkp1, R_hk_2hkp1, f1, f2, f3, f4
621 
622  r_hk_hkp1 = hk + hkp1
623  if (r_hk_hkp1 <= 0.) then
624  ppm_edge = 0.5 * ( ak + akp1 )
625  return
626  endif
627  r_hk_hkp1 = 1. / r_hk_hkp1
628  if (hk<hkp1) then
629  ppm_edge = ak + ( hk * r_hk_hkp1 ) * ( akp1 - ak )
630  else
631  ppm_edge = akp1 + ( hkp1 * r_hk_hkp1 ) * ( ak - akp1 )
632  endif
633 
634  r_2hk_hkp1 = 1. / ( ( 2. * hk + hkp1 ) + h_neglect )
635  r_hk_2hkp1 = 1. / ( ( hk + 2. * hkp1 ) + h_neglect )
636  f1 = 1./ ( ( hk + hkp1) + ( hkm1 + hkp2 ) )
637  f2 = 2. * ( hkp1 * hk ) * r_hk_hkp1 * &
638  ( ( hkm1 + hk ) * r_2hk_hkp1 - ( hkp2 + hkp1 ) * r_hk_2hkp1 )
639  f3 = hk * ( hkm1 + hk ) * r_2hk_hkp1
640  f4 = hkp1 * ( hkp1 + hkp2 ) * r_hk_2hkp1
641 
642  ppm_edge = ppm_edge + f1 * ( f2 * ( akp1 - ak ) - ( f3 * pkp1 - f4 * pk ) )
643 

◆ ppm_left_right_edge_values()

subroutine mom_neutral_diffusion::ppm_left_right_edge_values ( integer, intent(in)  nk,
real, dimension(nk), intent(in)  Tl,
real, dimension(nk+1), intent(in)  Ti,
real, dimension(nk), intent(inout)  aL,
real, dimension(nk), intent(inout)  aR 
)
private

Discontinuous PPM reconstructions of the left/right edge values within a cell.

Parameters
[in]nkNumber of levels
[in]tlLayer tracer (conc, e.g. degC)
[in]tiInterface tracer (conc, e.g. degC)
[in,out]alLeft edge value of tracer (conc, e.g. degC)
[in,out]arRight edge value of tracer (conc, e.g. degC)

Definition at line 1483 of file MOM_neutral_diffusion.F90.

1483  integer, intent(in) :: nk !< Number of levels
1484  real, dimension(nk), intent(in) :: Tl !< Layer tracer (conc, e.g. degC)
1485  real, dimension(nk+1), intent(in) :: Ti !< Interface tracer (conc, e.g. degC)
1486  real, dimension(nk), intent(inout) :: aL !< Left edge value of tracer (conc, e.g. degC)
1487  real, dimension(nk), intent(inout) :: aR !< Right edge value of tracer (conc, e.g. degC)
1488 
1489  integer :: k
1490  ! Setup reconstruction edge values
1491  do k = 1, nk
1492  al(k) = ti(k)
1493  ar(k) = ti(k+1)
1494  if ( signum(1., ar(k) - tl(k))*signum(1., tl(k) - al(k)) <= 0.0 ) then
1495  al(k) = tl(k)
1496  ar(k) = tl(k)
1497  elseif ( sign(3., ar(k) - al(k)) * ( (tl(k) - al(k)) + (tl(k) - ar(k))) > abs(ar(k) - al(k)) ) then
1498  al(k) = tl(k) + 2.0 * ( tl(k) - ar(k) )
1499  elseif ( sign(3., ar(k) - al(k)) * ( (tl(k) - al(k)) + (tl(k) - ar(k))) < -abs(ar(k) - al(k)) ) then
1500  ar(k) = tl(k) + 2.0 * ( tl(k) - al(k) )
1501  endif
1502  enddo

◆ signum()

real function mom_neutral_diffusion::signum ( real, intent(in)  a,
real, intent(in)  x 
)
private

A true signum function that returns either -abs(a), when x<0; or abs(a) when x>0; or 0 when x=0.

Parameters
[in]aThe magnitude argument
[in]xThe sign (or zero) argument

Definition at line 676 of file MOM_neutral_diffusion.F90.

676  real, intent(in) :: a !< The magnitude argument
677  real, intent(in) :: x !< The sign (or zero) argument
678 
679  signum = sign(a,x)
680  if (x==0.) signum = 0.
681 

◆ test_data1d()

logical function mom_neutral_diffusion::test_data1d ( logical, intent(in)  verbose,
integer, intent(in)  nk,
real, dimension(nk), intent(in)  Po,
real, dimension(nk), intent(in)  Ptrue,
character(len=*), intent(in)  title 
)
private

Returns true if comparison of Po and Ptrue fails, and conditionally writes results to stream.

Parameters
[in]verboseIf true, write results to stdout
[in]nkNumber of layers
[in]poCalculated answer
[in]ptrueTrue answer
[in]titleTitle for messages

Definition at line 2085 of file MOM_neutral_diffusion.F90.

2085  logical, intent(in) :: verbose !< If true, write results to stdout
2086  integer, intent(in) :: nk !< Number of layers
2087  real, dimension(nk), intent(in) :: Po !< Calculated answer
2088  real, dimension(nk), intent(in) :: Ptrue !< True answer
2089  character(len=*), intent(in) :: title !< Title for messages
2090 
2091  ! Local variables
2092  integer :: k, stdunit
2093 
2094  test_data1d = .false.
2095  do k = 1,nk
2096  if (po(k) /= ptrue(k)) test_data1d = .true.
2097  enddo
2098 
2099  if (test_data1d .or. verbose) then
2100  stdunit = 6
2101  if (test_data1d) stdunit = 0 ! In case of wrong results, write to error stream
2102  write(stdunit,'(a)') title
2103  do k = 1,nk
2104  if (po(k) /= ptrue(k)) then
2105  test_data1d = .true.
2106  write(stdunit,'(a,i2,2(x,a,f20.16),x,a,1pe22.15,x,a)') &
2107  'k=',k,'Po=',po(k),'Ptrue=',ptrue(k),'err=',po(k)-ptrue(k),'WRONG!'
2108  else
2109  if (verbose) &
2110  write(stdunit,'(a,i2,2(x,a,f20.16),x,a,1pe22.15)') &
2111  'k=',k,'Po=',po(k),'Ptrue=',ptrue(k),'err=',po(k)-ptrue(k)
2112  endif
2113  enddo
2114  endif
2115 

◆ test_data1di()

logical function mom_neutral_diffusion::test_data1di ( logical, intent(in)  verbose,
integer, intent(in)  nk,
integer, dimension(nk), intent(in)  Po,
integer, dimension(nk), intent(in)  Ptrue,
character(len=*), intent(in)  title 
)
private

Returns true if comparison of Po and Ptrue fails, and conditionally writes results to stream.

Parameters
[in]verboseIf true, write results to stdout
[in]nkNumber of layers
[in]poCalculated answer
[in]ptrueTrue answer
[in]titleTitle for messages

Definition at line 2120 of file MOM_neutral_diffusion.F90.

2120  logical, intent(in) :: verbose !< If true, write results to stdout
2121  integer, intent(in) :: nk !< Number of layers
2122  integer, dimension(nk), intent(in) :: Po !< Calculated answer
2123  integer, dimension(nk), intent(in) :: Ptrue !< True answer
2124  character(len=*), intent(in) :: title !< Title for messages
2125 
2126  ! Local variables
2127  integer :: k, stdunit
2128 
2129  test_data1di = .false.
2130  do k = 1,nk
2131  if (po(k) /= ptrue(k)) test_data1di = .true.
2132  enddo
2133 
2134  if (test_data1di .or. verbose) then
2135  stdunit = 6
2136  if (test_data1di) stdunit = 0 ! In case of wrong results, write to error stream
2137  write(stdunit,'(a)') title
2138  do k = 1,nk
2139  if (po(k) /= ptrue(k)) then
2140  test_data1di = .true.
2141  write(stdunit,'(a,i2,2(x,a,i5),x,a)') 'k=',k,'Io=',po(k),'Itrue=',ptrue(k),'WRONG!'
2142  else
2143  if (verbose) &
2144  write(stdunit,'(a,i2,2(x,a,i5))') 'k=',k,'Io=',po(k),'Itrue=',ptrue(k)
2145  endif
2146  enddo
2147  endif
2148 

◆ test_fv_diff()

logical function mom_neutral_diffusion::test_fv_diff ( logical, intent(in)  verbose,
real, intent(in)  hkm1,
real, intent(in)  hk,
real, intent(in)  hkp1,
real, intent(in)  Skm1,
real, intent(in)  Sk,
real, intent(in)  Skp1,
real, intent(in)  Ptrue,
character(len=*), intent(in)  title 
)
private

Returns true if a test of fv_diff() fails, and conditionally writes results to stream.

Parameters
[in]verboseIf true, write results to stdout
[in]hkm1Left cell width
[in]hkCenter cell width
[in]hkp1Right cell width
[in]skm1Left cell average value
[in]skCenter cell average value
[in]skp1Right cell average value
[in]ptrueTrue answer [Pa]
[in]titleTitle for messages

Definition at line 1989 of file MOM_neutral_diffusion.F90.

1989  logical, intent(in) :: verbose !< If true, write results to stdout
1990  real, intent(in) :: hkm1 !< Left cell width
1991  real, intent(in) :: hk !< Center cell width
1992  real, intent(in) :: hkp1 !< Right cell width
1993  real, intent(in) :: Skm1 !< Left cell average value
1994  real, intent(in) :: Sk !< Center cell average value
1995  real, intent(in) :: Skp1 !< Right cell average value
1996  real, intent(in) :: Ptrue !< True answer [Pa]
1997  character(len=*), intent(in) :: title !< Title for messages
1998 
1999  ! Local variables
2000  integer :: stdunit
2001  real :: Pret
2002 
2003  pret = fv_diff(hkm1, hk, hkp1, skm1, sk, skp1)
2004  test_fv_diff = (pret /= ptrue)
2005 
2006  if (test_fv_diff .or. verbose) then
2007  stdunit = 6
2008  if (test_fv_diff) stdunit = 0 ! In case of wrong results, write to error stream
2009  write(stdunit,'(a)') title
2010  if (test_fv_diff) then
2011  write(stdunit,'(2(x,a,f20.16),x,a)') 'pRet=',pret,'pTrue=',ptrue,'WRONG!'
2012  else
2013  write(stdunit,'(2(x,a,f20.16))') 'pRet=',pret,'pTrue=',ptrue
2014  endif
2015  endif
2016 

◆ test_fvlsq_slope()

logical function mom_neutral_diffusion::test_fvlsq_slope ( logical, intent(in)  verbose,
real, intent(in)  hkm1,
real, intent(in)  hk,
real, intent(in)  hkp1,
real, intent(in)  Skm1,
real, intent(in)  Sk,
real, intent(in)  Skp1,
real, intent(in)  Ptrue,
character(len=*), intent(in)  title 
)
private

Returns true if a test of fvlsq_slope() fails, and conditionally writes results to stream.

Parameters
[in]verboseIf true, write results to stdout
[in]hkm1Left cell width
[in]hkCenter cell width
[in]hkp1Right cell width
[in]skm1Left cell average value
[in]skCenter cell average value
[in]skp1Right cell average value
[in]ptrueTrue answer [Pa]
[in]titleTitle for messages

Definition at line 2021 of file MOM_neutral_diffusion.F90.

2021  logical, intent(in) :: verbose !< If true, write results to stdout
2022  real, intent(in) :: hkm1 !< Left cell width
2023  real, intent(in) :: hk !< Center cell width
2024  real, intent(in) :: hkp1 !< Right cell width
2025  real, intent(in) :: Skm1 !< Left cell average value
2026  real, intent(in) :: Sk !< Center cell average value
2027  real, intent(in) :: Skp1 !< Right cell average value
2028  real, intent(in) :: Ptrue !< True answer [Pa]
2029  character(len=*), intent(in) :: title !< Title for messages
2030 
2031  ! Local variables
2032  integer :: stdunit
2033  real :: Pret
2034 
2035  pret = fvlsq_slope(hkm1, hk, hkp1, skm1, sk, skp1)
2036  test_fvlsq_slope = (pret /= ptrue)
2037 
2038  if (test_fvlsq_slope .or. verbose) then
2039  stdunit = 6
2040  if (test_fvlsq_slope) stdunit = 0 ! In case of wrong results, write to error stream
2041  write(stdunit,'(a)') title
2042  if (test_fvlsq_slope) then
2043  write(stdunit,'(2(x,a,f20.16),x,a)') 'pRet=',pret,'pTrue=',ptrue,'WRONG!'
2044  else
2045  write(stdunit,'(2(x,a,f20.16))') 'pRet=',pret,'pTrue=',ptrue
2046  endif
2047  endif
2048 

◆ test_ifndp()

logical function mom_neutral_diffusion::test_ifndp ( logical, intent(in)  verbose,
real, intent(in)  rhoNeg,
real, intent(in)  Pneg,
real, intent(in)  rhoPos,
real, intent(in)  Ppos,
real, intent(in)  Ptrue,
character(len=*), intent(in)  title 
)
private

Returns true if a test of interpolate_for_nondim_position() fails, and conditionally writes results to stream.

Parameters
[in]verboseIf true, write results to stdout
[in]rhonegLighter density [kg m-3]
[in]pnegInterface position of lighter density [Pa]
[in]rhoposHeavier density [kg m-3]
[in]pposInterface position of heavier density [Pa]
[in]ptrueTrue answer [Pa]
[in]titleTitle for messages

Definition at line 2053 of file MOM_neutral_diffusion.F90.

2053  logical, intent(in) :: verbose !< If true, write results to stdout
2054  real, intent(in) :: rhoNeg !< Lighter density [kg m-3]
2055  real, intent(in) :: Pneg !< Interface position of lighter density [Pa]
2056  real, intent(in) :: rhoPos !< Heavier density [kg m-3]
2057  real, intent(in) :: Ppos !< Interface position of heavier density [Pa]
2058  real, intent(in) :: Ptrue !< True answer [Pa]
2059  character(len=*), intent(in) :: title !< Title for messages
2060 
2061  ! Local variables
2062  integer :: stdunit
2063  real :: Pret
2064 
2065  pret = interpolate_for_nondim_position(rhoneg, pneg, rhopos, ppos)
2066  test_ifndp = (pret /= ptrue)
2067 
2068  if (test_ifndp .or. verbose) then
2069  stdunit = 6
2070  if (test_ifndp) stdunit = 0 ! In case of wrong results, write to error stream
2071  write(stdunit,'(a)') title
2072  if (test_ifndp) then
2073  write(stdunit,'(4(x,a,f20.16),2(x,a,1pe22.15),x,a)') &
2074  'r1=',rhoneg,'p1=',pneg,'r2=',rhopos,'p2=',ppos,'pRet=',pret,'pTrue=',ptrue,'WRONG!'
2075  else
2076  write(stdunit,'(4(x,a,f20.16),2(x,a,1pe22.15))') &
2077  'r1=',rhoneg,'p1=',pneg,'r2=',rhopos,'p2=',ppos,'pRet=',pret,'pTrue=',ptrue
2078  endif
2079  endif
2080 

◆ test_nsp()

logical function mom_neutral_diffusion::test_nsp ( logical, intent(in)  verbose,
integer, intent(in)  ns,
integer, dimension(ns), intent(in)  KoL,
integer, dimension(ns), intent(in)  KoR,
real, dimension(ns), intent(in)  pL,
real, dimension(ns), intent(in)  pR,
real, dimension(ns-1), intent(in)  hEff,
integer, dimension(ns), intent(in)  KoL0,
integer, dimension(ns), intent(in)  KoR0,
real, dimension(ns), intent(in)  pL0,
real, dimension(ns), intent(in)  pR0,
real, dimension(ns-1), intent(in)  hEff0,
character(len=*), intent(in)  title 
)
private

Returns true if output of find_neutral_surface_positions() does not match correct values, and conditionally writes results to stream.

Parameters
[in]verboseIf true, write results to stdout
[in]nsNumber of surfaces
[in]kolIndex of first left interface above neutral surface
[in]korIndex of first right interface above neutral surface
[in]plFractional position of neutral surface within layer KoL of left column
[in]prFractional position of neutral surface within layer KoR of right column
[in]heffEffective thickness between two neutral surfaces [Pa]
[in]kol0Correct value for KoL
[in]kor0Correct value for KoR
[in]pl0Correct value for pL
[in]pr0Correct value for pR
[in]heff0Correct value for hEff
[in]titleTitle for messages

Definition at line 2154 of file MOM_neutral_diffusion.F90.

2154  logical, intent(in) :: verbose !< If true, write results to stdout
2155  integer, intent(in) :: ns !< Number of surfaces
2156  integer, dimension(ns), intent(in) :: KoL !< Index of first left interface above neutral surface
2157  integer, dimension(ns), intent(in) :: KoR !< Index of first right interface above neutral surface
2158  real, dimension(ns), intent(in) :: pL !< Fractional position of neutral surface within layer KoL of left column
2159  real, dimension(ns), intent(in) :: pR !< Fractional position of neutral surface within layer KoR of right column
2160  real, dimension(ns-1), intent(in) :: hEff !< Effective thickness between two neutral surfaces [Pa]
2161  integer, dimension(ns), intent(in) :: KoL0 !< Correct value for KoL
2162  integer, dimension(ns), intent(in) :: KoR0 !< Correct value for KoR
2163  real, dimension(ns), intent(in) :: pL0 !< Correct value for pL
2164  real, dimension(ns), intent(in) :: pR0 !< Correct value for pR
2165  real, dimension(ns-1), intent(in) :: hEff0 !< Correct value for hEff
2166  character(len=*), intent(in) :: title !< Title for messages
2167 
2168  ! Local variables
2169  integer :: k, stdunit
2170  logical :: this_row_failed
2171 
2172  test_nsp = .false.
2173  do k = 1,ns
2174  test_nsp = test_nsp .or. compare_nsp_row(kol(k), kor(k), pl(k), pr(k), kol0(k), kor0(k), pl0(k), pr0(k))
2175  if (k < ns) then
2176  if (heff(k) /= heff0(k)) test_nsp = .true.
2177  endif
2178  enddo
2179 
2180  if (test_nsp .or. verbose) then
2181  stdunit = 6
2182  if (test_nsp) stdunit = 0 ! In case of wrong results, write to error stream
2183  write(stdunit,'(a)') title
2184  do k = 1,ns
2185  this_row_failed = compare_nsp_row(kol(k), kor(k), pl(k), pr(k), kol0(k), kor0(k), pl0(k), pr0(k))
2186  if (this_row_failed) then
2187  write(stdunit,10) k,kol(k),pl(k),kor(k),pr(k),' <-- WRONG!'
2188  write(stdunit,10) k,kol0(k),pl0(k),kor0(k),pr0(k),' <-- should be this'
2189  else
2190  write(stdunit,10) k,kol(k),pl(k),kor(k),pr(k)
2191  endif
2192  if (k < ns) then
2193  if (heff(k) /= heff0(k)) then
2194  write(stdunit,'(i3,8x,"layer hEff =",2(f20.16,a))') k,heff(k)," .neq. ",heff0(k),' <-- WRONG!'
2195  else
2196  write(stdunit,'(i3,8x,"layer hEff =",f20.16)') k,heff(k)
2197  endif
2198  endif
2199  enddo
2200  endif
2201  if (test_nsp) call mom_error(fatal,"test_nsp failed")
2202 
2203 10 format("ks=",i3," kL=",i3," pL=",f20.16," kR=",i3," pR=",f20.16,a)

◆ test_rnp()

logical function mom_neutral_diffusion::test_rnp ( real, intent(in)  expected_pos,
real, intent(in)  test_pos,
character(len=*), intent(in)  title 
)
private

Compares output position from refine_nondim_position with an expected value.

Parameters
[in]expected_posThe expected position
[in]test_posThe position returned by the code
[in]titleA label for this test

Definition at line 2226 of file MOM_neutral_diffusion.F90.

2226  real, intent(in) :: expected_pos !< The expected position
2227  real, intent(in) :: test_pos !< The position returned by the code
2228  character(len=*), intent(in) :: title !< A label for this test
2229  ! Local variables
2230  integer :: stdunit = 6 ! Output to standard error
2231  test_rnp = expected_pos /= test_pos
2232  if (test_rnp) then
2233  write(stdunit,'(A, f20.16, " .neq. ", f20.16, " <-- WRONG")') title, expected_pos, test_pos
2234  else
2235  write(stdunit,'(A, f20.16, " == ", f20.16)') title, expected_pos, test_pos
2236  endif