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
 

Functions/Subroutines

logical function, public neutral_diffusion_init (Time, G, param_file, diag, CS)
 Read parameters and allocate control structure for neutral_diffusion module. More...
 
subroutine, public neutral_diffusion_diag_init (Time, G, diag, C_p, Reg, CS)
 Diagnostic handles for neutral diffusion tendencies. More...
 
subroutine, public neutral_diffusion_calc_coeffs (G, GV, h, T, S, EOS, 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, Tracer, m, dt, name, CS)
 Update tracer concentration due to neutral diffusion; layer thickness unchanged by this update. More...
 
subroutine interface_scalar (nk, h, S, Si, i_method)
 Returns interface scalar, Si, for a column of layer values, S. More...
 
real function ppm_edge (hkm1, hk, hkp1, hkp2, Ak, Akp1, Pk, Pkp1)
 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 (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. More...
 
real function absolute_position (n, Pint, Karr, NParr, k_surface)
 Converts non-dimensional position within a layer to absolute position (for debugging) More...
 
real function, dimension(2 *n+2) absolute_positions (n, Pint, Karr, NParr)
 Converts non-dimensional positions within layers to absolute positions (for debugging) More...
 
real function interpolate_for_nondim_position (dRhoNeg, Pneg, dRhoPos, Ppos)
 Returns the non-dimensional position between Pneg and Ppos where the interpolated density difference equals zero. The result is always bounded to be between 0 and 1. More...
 
subroutine neutral_surface_flux (nk, hl, hr, Tl, Tr, PiL, PiR, KoL, KoR, hEff, Flx)
 Returns a single column of neutral diffusion fluxes of a tracer. 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 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, nk, 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...
 
subroutine, public neutral_diffusion_end (CS)
 Deallocates neutral_diffusion control structure. More...
 

Variables

character(len=40) mdl = "MOM_neutral_diffusion"
 
logical, parameter debug_this_module = .false.
 

Function/Subroutine Documentation

◆ absolute_position()

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

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

Parameters
[in]nNumber of levels
[in]pintPosition of interfaces (Pa)
[in]karrIndex of interface above position
[in]nparrNon-dimensional position within layer Karr(:)

Definition at line 892 of file MOM_neutral_diffusion.F90.

Referenced by absolute_positions(), and find_neutral_surface_positions().

892  integer, intent(in) :: n !< Number of levels
893  real, intent(in) :: pint(n+1) !< Position of interfaces (Pa)
894  integer, intent(in) :: karr(2*n+2) !< Index of interface above position
895  real, intent(in) :: nparr(2*n+2) !< Non-dimensional position within layer Karr(:)
896 
897  ! Local variables
898  integer :: k_surface, k
899 
900  k = karr(k_surface)
901  if (k>n) stop 'absolute_position: k>nk is out of bounds!'
902  absolute_position = pint(k) + nparr(k_surface) * ( pint(k+1) - pint(k) )
903 
Here is the caller graph for this function:

◆ absolute_positions()

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

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

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

Definition at line 908 of file MOM_neutral_diffusion.F90.

References absolute_position().

Referenced by neutral_diffusion_unit_tests().

908  integer, intent(in) :: n !< Number of levels
909  real, intent(in) :: pint(n+1) !< Position of interface (Pa)
910  integer, intent(in) :: karr(2*n+2) !< Indexes of interfaces about positions
911  real, intent(in) :: nparr(2*n+2) !< Non-dimensional positions within layers Karr(:)
912 
913  real, dimension(2*n+2) :: absolute_positions ! Absolute positions (Pa)
914 
915  ! Local variables
916  integer :: k_surface, k
917 
918  do k_surface = 1, 2*n+2
919  absolute_positions(k_surface) = absolute_position(n,pint,karr,nparr,k_surface)
920  enddo
921 
Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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 1495 of file MOM_neutral_diffusion.F90.

Referenced by test_nsp().

1495  integer, intent(in) :: kol !< Index of first left interface above neutral surface
1496  integer, intent(in) :: kor !< Index of first right interface above neutral surface
1497  real, intent(in) :: pl !< Fractional position of neutral surface within layer KoL of left column
1498  real, intent(in) :: pr !< Fractional position of neutral surface within layer KoR of right column
1499  integer, intent(in) :: kol0 !< Correct value for KoL
1500  integer, intent(in) :: kor0 !< Correct value for KoR
1501  real, intent(in) :: pl0 !< Correct value for pL
1502  real, intent(in) :: pr0 !< Correct value for pR
1503 
1504  compare_nsp_row = .false.
1505  if (kol /= kol0) compare_nsp_row = .true.
1506  if (kor /= kor0) compare_nsp_row = .true.
1507  if (pl /= pl0) compare_nsp_row = .true.
1508  if (pr /= pr0) compare_nsp_row = .true.
Here is the caller graph for this function:

◆ find_neutral_surface_positions()

subroutine mom_neutral_diffusion::find_neutral_surface_positions ( 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.

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/m3/degC)
[in]drdslLeft-column dRho/dS (kg/m3/ppt)
[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/m3/degC)
[in]drdsrLeft-column dRho/dS (kg/m3/ppt)
[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 716 of file MOM_neutral_diffusion.F90.

References absolute_position(), and interpolate_for_nondim_position().

Referenced by neutral_diffusion_calc_coeffs(), and neutral_diffusion_unit_tests().

716  integer, intent(in) :: nk !< Number of levels
717  real, dimension(nk+1), intent(in) :: pl !< Left-column interface pressure (Pa)
718  real, dimension(nk+1), intent(in) :: tl !< Left-column interface potential temperature (degC)
719  real, dimension(nk+1), intent(in) :: sl !< Left-column interface salinity (ppt)
720  real, dimension(nk+1), intent(in) :: drdtl !< Left-column dRho/dT (kg/m3/degC)
721  real, dimension(nk+1), intent(in) :: drdsl !< Left-column dRho/dS (kg/m3/ppt)
722  real, dimension(nk+1), intent(in) :: pr !< Right-column interface pressure (Pa)
723  real, dimension(nk+1), intent(in) :: tr !< Right-column interface potential temperature (degC)
724  real, dimension(nk+1), intent(in) :: sr !< Right-column interface salinity (ppt)
725  real, dimension(nk+1), intent(in) :: drdtr !< Left-column dRho/dT (kg/m3/degC)
726  real, dimension(nk+1), intent(in) :: drdsr !< Left-column dRho/dS (kg/m3/ppt)
727  real, dimension(2*nk+2), intent(inout) :: pol !< Fractional position of neutral surface within layer KoL of left column
728  real, dimension(2*nk+2), intent(inout) :: por !< Fractional position of neutral surface within layer KoR of right column
729  integer, dimension(2*nk+2), intent(inout) :: kol !< Index of first left interface above neutral surface
730  integer, dimension(2*nk+2), intent(inout) :: kor !< Index of first right interface above neutral surface
731  real, dimension(2*nk+1), intent(inout) :: heff !< Effective thickness between two neutral surfaces (Pa)
732 
733  ! Local variables
734  integer :: k_surface ! Index of neutral surface
735  integer :: kl ! Index of left interface
736  integer :: kr ! Index of right interface
737  real :: drdt, drds ! dRho/dT and dRho/dS for the neutral surface
738  logical :: searching_left_column ! True if searching for the position of a right interface in the left column
739  logical :: searching_right_column ! True if searching for the position of a left interface in the right column
740  logical :: reached_bottom ! True if one of the bottom-most interfaces has been used as the target
741  integer :: krm1, klm1
742  real :: drho, drhotop, drhobot, hl, hr
743  integer :: lastk_left, lastk_right
744  real :: lastp_left, lastp_right
745 
746  ! Initialize variables for the search
747  kr = 1 ; lastk_right = 1 ; lastp_right = 0.
748  kl = 1 ; lastk_left = 1 ; lastp_left = 0.
749  reached_bottom = .false.
750 
751  ! Loop over each neutral surface, working from top to bottom
752  neutral_surfaces: do k_surface = 1, 2*nk+2
753  klm1 = max(kl-1, 1)
754  if (klm1>nk) stop 'find_neutral_surface_positions(): klm1 went out of bounds!'
755  krm1 = max(kr-1, 1)
756  if (krm1>nk) stop 'find_neutral_surface_positions(): krm1 went out of bounds!'
757 
758  ! Potential density difference, rho(kr) - rho(kl)
759  drho = 0.5 * ( ( drdtr(kr) + drdtl(kl) ) * ( tr(kr) - tl(kl) ) &
760  + ( drdsr(kr) + drdsl(kl) ) * ( sr(kr) - sl(kl) ) )
761  ! Which column has the lighter surface for the current indexes, kr and kl
762  if (.not. reached_bottom) then
763  if (drho < 0.) then
764  searching_left_column = .true.
765  searching_right_column = .false.
766  elseif (drho > 0.) then
767  searching_right_column = .true.
768  searching_left_column = .false.
769  else ! dRho == 0.
770  if (kl + kr == 2) then ! Still at surface
771  searching_left_column = .true.
772  searching_right_column = .false.
773  else ! Not the surface so we simply change direction
774  searching_left_column = .not. searching_left_column
775  searching_right_column = .not. searching_right_column
776  endif
777  endif
778  endif
779 
780  if (searching_left_column) then
781  ! Interpolate for the neutral surface position within the left column, layer klm1
782  ! Potential density difference, rho(kl-1) - rho(kr) (should be negative)
783  drhotop = 0.5 * ( ( drdtl(klm1) + drdtr(kr) ) * ( tl(klm1) - tr(kr) ) &
784  + ( drdsl(klm1) + drdsr(kr) ) * ( sl(klm1) - sr(kr) ) )
785  ! Potential density difference, rho(kl) - rho(kr) (will be positive)
786  drhobot = 0.5 * ( ( drdtl(klm1+1) + drdtr(kr) ) * ( tl(klm1+1) - tr(kr) ) &
787  + ( drdsl(klm1+1) + drdsr(kr) ) * ( sl(klm1+1) - sr(kr) ) )
788 
789  ! Because we are looking left, the right surface, kr, is lighter than klm1+1 and should be denser than klm1
790  ! unless we are still at the top of the left column (kl=1)
791  if (drhotop > 0. .or. kr+kl==2) then
792  pol(k_surface) = 0. ! The right surface is lighter than anything in layer klm1
793  elseif (drhotop >= drhobot) then ! Left layer is unstratified
794  pol(k_surface) = 1.
795  else
796  ! Linearly interpolate for the position between Pl(kl-1) and Pl(kl) where the density difference
797  ! between right and left is zero.
798  pol(k_surface) = interpolate_for_nondim_position( drhotop, pl(klm1), drhobot, pl(klm1+1) )
799  endif
800  if (pol(k_surface)>=1. .and. klm1<nk) then ! >= is really ==, when PoL==1 we point to the bottom of the cell
801  klm1 = klm1 + 1
802  pol(k_surface) = pol(k_surface) - 1.
803  endif
804  if (real(klm1-lastk_left)+(pol(k_surface)-lastp_left)<0.) then
805  pol(k_surface) = lastp_left
806  klm1 = lastk_left
807  endif
808  kol(k_surface) = klm1
809  if (kr <= nk) then
810  por(k_surface) = 0.
811  kor(k_surface) = kr
812  else
813  por(k_surface) = 1.
814  kor(k_surface) = nk
815  endif
816  if (kr <= nk) then
817  kr = kr + 1
818  else
819  reached_bottom = .true.
820  searching_right_column = .true.
821  searching_left_column = .false.
822  endif
823  elseif (searching_right_column) then
824  ! Interpolate for the neutral surface position within the right column, layer krm1
825  ! Potential density difference, rho(kr-1) - rho(kl) (should be negative)
826  drhotop = 0.5 * ( ( drdtr(krm1) + drdtl(kl) ) * ( tr(krm1) - tl(kl) ) &
827  + ( drdsr(krm1) + drdsl(kl) ) * ( sr(krm1) - sl(kl) ) )
828  ! Potential density difference, rho(kr) - rho(kl) (will be positive)
829  drhobot = 0.5 * ( ( drdtr(krm1+1) + drdtl(kl) ) * ( tr(krm1+1) - tl(kl) ) &
830  + ( drdsr(krm1+1) + drdsl(kl) ) * ( sr(krm1+1) - sl(kl) ) )
831 
832  ! Because we are looking right, the left surface, kl, is lighter than krm1+1 and should be denser than krm1
833  ! unless we are still at the top of the right column (kr=1)
834  if (drhotop >= 0. .or. kr+kl==2) then
835  por(k_surface) = 0. ! The left surface is lighter than anything in layer krm1
836  elseif (drhotop >= drhobot) then ! Right layer is unstratified
837  por(k_surface) = 1.
838  else
839  ! Linearly interpolate for the position between Pr(kr-1) and Pr(kr) where the density difference
840  ! between right and left is zero.
841  por(k_surface) = interpolate_for_nondim_position( drhotop, pr(krm1), drhobot, pr(krm1+1) )
842  endif
843  if (por(k_surface)>=1. .and. krm1<nk) then ! >= is really ==, when PoR==1 we point to the bottom of the cell
844  krm1 = krm1 + 1
845  por(k_surface) = por(k_surface) - 1.
846  endif
847  if (real(krm1-lastk_right)+(por(k_surface)-lastp_right)<0.) then
848  por(k_surface) = lastp_right
849  krm1 = lastk_right
850  endif
851  kor(k_surface) = krm1
852  if (kl <= nk) then
853  pol(k_surface) = 0.
854  kol(k_surface) = kl
855  else
856  pol(k_surface) = 1.
857  kol(k_surface) = nk
858  endif
859  if (kl <= nk) then
860  kl = kl + 1
861  else
862  reached_bottom = .true.
863  searching_right_column = .false.
864  searching_left_column = .true.
865  endif
866  else
867  stop 'Else what?'
868  endif
869 
870  lastk_left = kol(k_surface) ; lastp_left = pol(k_surface)
871  lastk_right = kor(k_surface) ; lastp_right = por(k_surface)
872 
873  ! Effective thickness
874  ! NOTE: This would be better expressed in terms of the layers thicknesses rather
875  ! than as differences of position - AJA
876  if (k_surface>1) then
877  hl = absolute_position(nk,pl,kol,pol,k_surface) - absolute_position(nk,pl,kol,pol,k_surface-1)
878  hr = absolute_position(nk,pr,kor,por,k_surface) - absolute_position(nk,pr,kor,por,k_surface-1)
879  if ( hl + hr > 0.) then
880  heff(k_surface-1) = 2. * hl * hr / ( hl + hr )
881  else
882  heff(k_surface-1) = 0.
883  endif
884  endif
885 
886  enddo neutral_surfaces
887 
Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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 660 of file MOM_neutral_diffusion.F90.

Referenced by plm_diff(), and test_fv_diff().

660  real, intent(in) :: hkm1 !< Left cell width
661  real, intent(in) :: hk !< Center cell width
662  real, intent(in) :: hkp1 !< Right cell width
663  real, intent(in) :: skm1 !< Left cell average value
664  real, intent(in) :: sk !< Center cell average value
665  real, intent(in) :: skp1 !< Right cell average value
666 
667  ! Local variables
668  real :: h_sum, hp, hm
669 
670  h_sum = ( hkm1 + hkp1 ) + hk
671  if (h_sum /= 0.) h_sum = 1./ h_sum
672  hm = hkm1 + hk
673  if (hm /= 0.) hm = 1./ hm
674  hp = hkp1 + hk
675  if (hp /= 0.) hp = 1./ hp
676  fv_diff = ( hk * h_sum ) * &
677  ( ( 2. * hkm1 + hk ) * hp * ( skp1 - sk ) &
678  + ( 2. * hkp1 + hk ) * hm * ( sk - skm1 ) )
Here is the caller graph for this function:

◆ 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 686 of file MOM_neutral_diffusion.F90.

Referenced by plm_diff(), and test_fvlsq_slope().

686  real, intent(in) :: hkm1 !< Left cell width
687  real, intent(in) :: hk !< Center cell width
688  real, intent(in) :: hkp1 !< Right cell width
689  real, intent(in) :: skm1 !< Left cell average value
690  real, intent(in) :: sk !< Center cell average value
691  real, intent(in) :: skp1 !< Right cell average value
692 
693  ! Local variables
694  real :: xkm1, xkp1
695  real :: h_sum, hx_sum, hxsq_sum, hxy_sum, hy_sum, det
696 
697  xkm1 = -0.5 * ( hk + hkm1 )
698  xkp1 = 0.5 * ( hk + hkp1 )
699  h_sum = ( hkm1 + hkp1 ) + hk
700  hx_sum = hkm1*xkm1 + hkp1*xkp1
701  hxsq_sum = hkm1*(xkm1**2) + hkp1*(xkp1**2)
702  hxy_sum = hkm1*xkm1*skm1 + hkp1*xkp1*skp1
703  hy_sum = ( hkm1*skm1 + hkp1*skp1 ) + hk*sk
704  det = h_sum * hxsq_sum - hx_sum**2
705  if (det /= 0.) then
706  !a = ( hxsq_sum * hy_sum - hx_sum*hxy_sum ) / det ! a would be mean of straight line fit
707  fvlsq_slope = ( h_sum * hxy_sum - hx_sum*hy_sum ) / det ! Gradient of straight line fit
708  else
709  fvlsq_slope = 0. ! Adcroft's reciprocal rule
710  endif
Here is the caller graph for this function:

◆ 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 
)
private

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

Parameters
[in]nkNumber of levels
[in]hLayer thickness (H units)
[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

Definition at line 477 of file MOM_neutral_diffusion.F90.

References plm_diff(), and ppm_edge().

Referenced by neutral_diffusion_calc_coeffs(), neutral_diffusion_unit_tests(), and neutral_surface_flux().

477  integer, intent(in) :: nk !< Number of levels
478  real, dimension(nk), intent(in) :: h !< Layer thickness (H units)
479  real, dimension(nk), intent(in) :: s !< Layer scalar (conc, e.g. ppt)
480  real, dimension(nk+1), intent(inout) :: si !< Interface scalar (conc, e.g. ppt)
481  integer, intent(in) :: i_method !< =1 use average of PLM edges
482  !! =2 use continuous PPM edge interpolation
483  ! Local variables
484  integer :: k, km2, kp1
485  real, dimension(nk) :: diff
486  real :: sb, sa
487 
488  call plm_diff(nk, h, s, 2, 1, diff)
489  si(1) = s(1) - 0.5 * diff(1)
490  if (i_method==1) then
491  do k = 2, nk
492  ! Average of the two edge values (will be bounded and,
493  ! when slopes are unlimited, notionally second-order accurate)
494  sa = s(k-1) + 0.5 * diff(k-1) ! Lower edge value of a PLM reconstruction for layer above
495  sb = s(k) - 0.5 * diff(k) ! Upper edge value of a PLM reconstruction for layer below
496  si(k) = 0.5 * ( sa + sb )
497  enddo
498  elseif (i_method==2) then
499  do k = 2, nk
500  ! PPM quasi-fourth order interpolation for edge values following
501  ! equation 1.6 in Colella & Woodward, 1984: JCP 54, 174-201.
502  km2 = max(1, k-2)
503  kp1 = min(nk, k+1)
504  si(k) = ppm_edge(h(km2), h(k-1), h(k), h(kp1), s(k-1), s(k), diff(k-1), diff(k))
505  enddo
506  endif
507  si(nk+1) = s(nk) + 0.5 * diff(nk)
508 
Here is the call graph for this function:
Here is the caller graph for this function:

◆ interpolate_for_nondim_position()

real function mom_neutral_diffusion::interpolate_for_nondim_position ( real, intent(in)  dRhoNeg,
real, intent(in)  Pneg,
real, intent(in)  dRhoPos,
real, intent(in)  Ppos 
)
private

Returns the non-dimensional position between Pneg and Ppos where the interpolated density difference equals zero. The result is always bounded to be between 0 and 1.

Parameters
[in]drhonegNegative density difference
[in]pnegPosition of negative density difference
[in]drhoposPositive density difference
[in]pposPosition of positive density difference

Definition at line 928 of file MOM_neutral_diffusion.F90.

Referenced by find_neutral_surface_positions(), and test_ifndp().

928  real, intent(in) :: drhoneg !< Negative density difference
929  real, intent(in) :: pneg !< Position of negative density difference
930  real, intent(in) :: drhopos !< Positive density difference
931  real, intent(in) :: ppos !< Position of positive density difference
932 
933  if (ppos<pneg) stop 'interpolate_for_nondim_position: Houston, we have a problem! Ppos<Pneg'
934  if (drhoneg>drhopos) write(0,*) 'dRhoNeg, Pneg, dRhoPos, Ppos=',drhoneg, pneg, drhopos, ppos
935  if (drhoneg>drhopos) stop 'interpolate_for_nondim_position: Houston, we have a problem! dRhoNeg>dRhoPos'
936  if (ppos<=pneg) then ! Handle vanished or inverted layers
937  interpolate_for_nondim_position = 0.5
938  elseif ( drhopos - drhoneg > 0. ) then
939  interpolate_for_nondim_position = min( 1., max( 0., -drhoneg / ( drhopos - drhoneg ) ) )
940  elseif ( drhopos - drhoneg == 0) then
941  if (drhoneg>0.) then
942  interpolate_for_nondim_position = 0.
943  elseif (drhoneg<0.) then
944  interpolate_for_nondim_position = 1.
945  else ! dRhoPos = dRhoNeg = 0
946  interpolate_for_nondim_position = 0.5
947  endif
948  else ! dRhoPos - dRhoNeg < 0
949  interpolate_for_nondim_position = 0.5
950  endif
951  if ( interpolate_for_nondim_position < 0. ) stop 'interpolate_for_nondim_position: Houston, we have a problem! Pint < Pneg'
952  if ( interpolate_for_nondim_position > 1. ) stop 'interpolate_for_nondim_position: Houston, we have a problem! Pint > Ppos'
Here is the caller graph for this function:

◆ 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, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(inout)  Tracer,
integer, intent(in)  m,
real, intent(in)  dt,
character(len=32), intent(in)  name,
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 units)
[in]coef_xdt * Kh * dy / dx at u-points (m^2)
[in]coef_ydt * Kh * dx / dy at u-points (m^2)
[in,out]tracerTracer concentration
[in]mTracer number
[in]dtTracer time step * I_numitts (I_numitts in tracer_hordiff)
[in]nameTracer name
csNeutral diffusion control structure

Definition at line 316 of file MOM_neutral_diffusion.F90.

References neutral_surface_flux().

Referenced by mom_tracer_hor_diff::tracer_hordiff().

316  type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
317  type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
318  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness (H units)
319  real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: coef_x !< dt * Kh * dy / dx at u-points (m^2)
320  real, dimension(SZI_(G),SZJB_(G)), intent(in) :: coef_y !< dt * Kh * dx / dy at u-points (m^2)
321  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: tracer !< Tracer concentration
322  integer, intent(in) :: m !< Tracer number
323  real, intent(in) :: dt !< Tracer time step * I_numitts (I_numitts in tracer_hordiff)
324  character(len=32), intent(in) :: name !< Tracer name
325  type(neutral_diffusion_cs), pointer :: cs !< Neutral diffusion control structure
326 
327  ! Local variables
328  real, dimension(SZIB_(G),SZJ_(G),2*G%ke+1) :: uflx ! Zonal flux of tracer (concentration * H)
329  real, dimension(SZI_(G),SZJB_(G),2*G%ke+1) :: vflx ! Meridional flux of tracer (concentration * H)
330  real, dimension(SZI_(G),SZJ_(G),G%ke) :: tendency ! tendency array for diagn
331  real, dimension(SZI_(G),SZJ_(G)) :: tendency_2d ! depth integrated content tendency for diagn
332  real, dimension(SZIB_(G),SZJ_(G)) :: trans_x_2d ! depth integrated diffusive tracer x-transport diagn
333  real, dimension(SZI_(G),SZJB_(G)) :: trans_y_2d ! depth integrated diffusive tracer y-transport diagn
334  real, dimension(G%ke) :: dtracer ! change in tracer concentration due to ndiffusion
335  integer :: i, j, k, ks, nk
336  real :: ppt2mks, idt, convert
337 
338  nk = gv%ke
339 
340  ! for diagnostics
341  if(cs%id_neutral_diff_tracer_conc_tend(m) > 0 .or. &
342  cs%id_neutral_diff_tracer_cont_tend(m) > 0 .or. &
343  cs%id_neutral_diff_tracer_cont_tend_2d(m) > 0 .or. &
344  cs%id_neutral_diff_tracer_trans_x_2d(m) > 0 .or. &
345  cs%id_neutral_diff_tracer_trans_y_2d(m) > 0) then
346  ppt2mks = 0.001
347  idt = 1.0/dt
348  tendency(:,:,:) = 0.0
349  tendency_2d(:,:) = 0.0
350  trans_x_2d(:,:) = 0.0
351  trans_y_2d(:,:) = 0.0
352  convert = 1.0
353  if(trim(name) == 'T') convert = cs%C_p * gv%H_to_kg_m2
354  if(trim(name) == 'S') convert = ppt2mks * gv%H_to_kg_m2
355  endif
356 
357 
358  ! x-flux
359  do j = g%jsc,g%jec ; do i = g%isc-1,g%iec
360  if (g%mask2dCu(i,j)>0.) then
361  call neutral_surface_flux(nk, h(i,j,:), h(i+1,j,:), &
362  tracer(i,j,:), tracer(i+1,j,:), &
363  cs%uPoL(i,j,:), cs%uPoR(i,j,:), &
364  cs%uKoL(i,j,:), cs%uKoR(i,j,:), &
365  cs%uhEff(i,j,:), uflx(i,j,:))
366  else
367  uflx(i,j,:) = 0.
368  endif
369  enddo ; enddo
370 
371  ! y-flux
372  do j = g%jsc-1,g%jec ; do i = g%isc,g%iec
373  if (g%mask2dCv(i,j)>0.) then
374  call neutral_surface_flux(nk, h(i,j,:), h(i,j+1,:), &
375  tracer(i,j,:), tracer(i,j+1,:), &
376  cs%vPoL(i,j,:), cs%vPoR(i,j,:), &
377  cs%vKoL(i,j,:), cs%vKoR(i,j,:), &
378  cs%vhEff(i,j,:), vflx(i,j,:))
379  else
380  vflx(i,j,:) = 0.
381  endif
382  enddo ; enddo
383 
384  ! Update the tracer concentration from divergence of neutral diffusive flux components
385  do j = g%jsc,g%jec ; do i = g%isc,g%iec
386  if (g%mask2dT(i,j)>0.) then
387 
388  dtracer(:) = 0.
389  do ks = 1,2*nk+1 ;
390  k = cs%uKoL(i,j,ks)
391  dtracer(k) = dtracer(k) + coef_x(i,j) * uflx(i,j,ks)
392  k = cs%uKoR(i-1,j,ks)
393  dtracer(k) = dtracer(k) - coef_x(i-1,j) * uflx(i-1,j,ks)
394  k = cs%vKoL(i,j,ks)
395  dtracer(k) = dtracer(k) + coef_y(i,j) * vflx(i,j,ks)
396  k = cs%vKoR(i,j-1,ks)
397  dtracer(k) = dtracer(k) - coef_y(i,j-1) * vflx(i,j-1,ks)
398  enddo
399  do k = 1, gv%ke
400  tracer(i,j,k) = tracer(i,j,k) + dtracer(k) * &
401  ( g%IareaT(i,j) / ( h(i,j,k) + gv%H_subroundoff ) )
402  enddo
403 
404  if(cs%id_neutral_diff_tracer_conc_tend(m) > 0 .or. &
405  cs%id_neutral_diff_tracer_cont_tend(m) > 0 .or. &
406  cs%id_neutral_diff_tracer_cont_tend_2d(m) > 0 ) then
407  do k = 1, gv%ke
408  tendency(i,j,k) = dtracer(k) * g%IareaT(i,j) * idt
409  enddo
410  endif
411 
412  endif
413  enddo ; enddo
414 
415 
416  ! Diagnose vertically summed zonal flux, giving zonal tracer transport from ndiff.
417  ! Note sign corresponds to downgradient flux convention.
418  if(cs%id_neutral_diff_tracer_trans_x_2d(m) > 0) then
419  do j = g%jsc,g%jec ; do i = g%isc-1,g%iec
420  trans_x_2d(i,j) = 0.
421  if (g%mask2dCu(i,j)>0.) then
422  do ks = 1,2*nk+1 ;
423  trans_x_2d(i,j) = trans_x_2d(i,j) - coef_x(i,j) * uflx(i,j,ks)
424  enddo
425  trans_x_2d(i,j) = trans_x_2d(i,j) * idt * convert
426  endif
427  enddo ; enddo
428  call post_data(cs%id_neutral_diff_tracer_trans_x_2d(m), trans_x_2d(:,:), cs%diag)
429  endif
430 
431  ! Diagnose vertically summed merid flux, giving meridional tracer transport from ndiff.
432  ! Note sign corresponds to downgradient flux convention.
433  if(cs%id_neutral_diff_tracer_trans_y_2d(m) > 0) then
434  do j = g%jsc-1,g%jec ; do i = g%isc,g%iec
435  trans_y_2d(i,j) = 0.
436  if (g%mask2dCv(i,j)>0.) then
437  do ks = 1,2*nk+1 ;
438  trans_y_2d(i,j) = trans_y_2d(i,j) - coef_y(i,j) * vflx(i,j,ks)
439  enddo
440  trans_y_2d(i,j) = trans_y_2d(i,j) * idt * convert
441  endif
442  enddo ; enddo
443  call post_data(cs%id_neutral_diff_tracer_trans_y_2d(m), trans_y_2d(:,:), cs%diag)
444  endif
445 
446  ! post tendency of tracer content
447  if(cs%id_neutral_diff_tracer_cont_tend(m) > 0) then
448  call post_data(cs%id_neutral_diff_tracer_cont_tend(m), tendency(:,:,:)*convert, cs%diag)
449  endif
450 
451  ! post depth summed tendency for tracer content
452  if(cs%id_neutral_diff_tracer_cont_tend_2d(m) > 0) then
453  do j = g%jsc,g%jec ; do i = g%isc,g%iec
454  do k = 1, gv%ke
455  tendency_2d(i,j) = tendency_2d(i,j) + tendency(i,j,k)
456  enddo
457  enddo ; enddo
458  call post_data(cs%id_neutral_diff_tracer_cont_tend_2d(m), tendency_2d(:,:)*convert, cs%diag)
459  endif
460 
461  ! post tendency of tracer concentration; this step must be
462  ! done after posting tracer content tendency, since we alter
463  ! the tendency array.
464  if(cs%id_neutral_diff_tracer_conc_tend(m) > 0) then
465  do k = 1, gv%ke ; do j = g%jsc,g%jec ; do i = g%isc,g%iec
466  tendency(i,j,k) = tendency(i,j,k) / ( h(i,j,k) + gv%H_subroundoff )
467  enddo ; enddo ; enddo
468  call post_data(cs%id_neutral_diff_tracer_conc_tend(m), tendency, cs%diag)
469  endif
470 
471 
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:19
Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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(eos_type), pointer  EOS,
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 units)
[in]tPotential temperature (degC)
[in]sSalinity (ppt)
eosEquation of state structure
csNeutral diffusion control structure

Definition at line 256 of file MOM_neutral_diffusion.F90.

References mom_eos::calculate_density_derivs(), find_neutral_surface_positions(), and interface_scalar().

Referenced by mom_tracer_hor_diff::tracer_hordiff().

256  type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
257  type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
258  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness (H units)
259  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: t !< Potential temperature (degC)
260  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: s !< Salinity (ppt)
261  type(eos_type), pointer :: eos !< Equation of state structure
262  type(neutral_diffusion_cs), pointer :: cs !< Neutral diffusion control structure
263 
264  ! Local variables
265  integer :: i, j, k
266  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: tint ! Interface T (degC)
267  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: sint ! Interface S (ppt)
268  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: pint ! Interface pressure (Pa)
269  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: drdt ! Interface thermal expansion coefficient (kg/m3/degC)
270  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: drds ! Interface haline expansion coefficient (kg/m3/ppt)
271 
272  do j = g%jsc-1, g%jec+1
273  ! Interpolate state to interface
274  do i = g%isc-1, g%iec+1
275  call interface_scalar(g%ke, h(i,j,:), t(i,j,:), tint(i,j,:), 2)
276  call interface_scalar(g%ke, h(i,j,:), s(i,j,:), sint(i,j,:), 2)
277  enddo
278 
279  ! Calculate interface properties
280  pint(:,j,1) = 0. ! Assume P=0 (Pa) at surface - needs correcting for atmospheric and ice loading - AJA
281  do k = 1, g%ke+1
282  call calculate_density_derivs(tint(:,j,k), sint(:,j,k), pint(:,j,k), &
283  drdt(:,j,k), drds(:,j,k), g%isc-1, g%iec-g%isc+3, eos)
284  if (k<=g%ke) pint(:,j,k+1) = pint(:,j,k) + h(:,j,k) * gv%H_to_Pa ! Pressure at next interface, k+1 (Pa)
285  enddo
286  enddo
287 
288  ! Neutral surface factors at U points
289  do j = g%jsc, g%jec
290  do i = g%isc-1, g%iec
291  call find_neutral_surface_positions(g%ke, &
292  pint(i,j,:), tint(i,j,:), sint(i,j,:), drdt(i,j,:), drds(i,j,:), &
293  pint(i+1,j,:), tint(i+1,j,:), sint(i+1,j,:), drdt(i+1,j,:), drds(i+1,j,:),&
294  cs%uPoL(i,j,:), cs%uPoR(i,j,:), cs%uKoL(i,j,:), cs%uKoR(i,j,:), cs%uhEff(i,j,:) )
295  enddo
296  enddo
297 
298  ! Neutral surface factors at V points
299  do j = g%jsc-1, g%jec
300  do i = g%isc, g%iec
301  call find_neutral_surface_positions(g%ke, &
302  pint(i,j,:), tint(i,j,:), sint(i,j,:), drdt(i,j,:), drds(i,j,:), &
303  pint(i,j+1,:), tint(i,j+1,:), sint(i,j+1,:), drdt(i,j+1,:), drds(i,j+1,:),&
304  cs%vPoL(i,j,:), cs%vPoR(i,j,:), cs%vKoL(i,j,:), cs%vKoR(i,j,:), cs%vhEff(i,j,:) )
305  enddo
306  enddo
307 
308  cs%uhEff(:,:,:) = cs%uhEff(:,:,:) / gv%H_to_pa
309  cs%vhEff(:,:,:) = cs%vhEff(:,:,:) / gv%H_to_pa
310 
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:19
Here is the call graph for this function:
Here is the caller graph for this function:

◆ neutral_diffusion_diag_init()

subroutine, public mom_neutral_diffusion::neutral_diffusion_diag_init ( type(time_type), intent(in), target  Time,
type(ocean_grid_type), intent(in)  G,
type(diag_ctrl), intent(in)  diag,
real, intent(in)  C_p,
type(tracer_registry_type), intent(in)  Reg,
type(neutral_diffusion_cs), pointer  CS 
)

Diagnostic handles for neutral diffusion tendencies.

Parameters
[in]timeTime structure
[in]gGrid structure
[in]diagDiagnostics control structure
[in]regTracer structure
[in]c_pSeawater heat capacity
csNeutral diffusion control structure

Definition at line 119 of file MOM_neutral_diffusion.F90.

References mom_diag_mediator::register_diag_field().

Referenced by mom::initialize_mom().

119  type(time_type),target, intent(in) :: time !< Time structure
120  type(ocean_grid_type), intent(in) :: g !< Grid structure
121  type(diag_ctrl), intent(in) :: diag !< Diagnostics control structure
122  type(tracer_registry_type), intent(in) :: reg !< Tracer structure
123  real, intent(in) :: c_p !< Seawater heat capacity
124  type(neutral_diffusion_cs), pointer :: cs !< Neutral diffusion control structure
125 
126  ! local
127  integer :: n,ntr
128 
129  if(.not. associated(cs)) return
130 
131  ntr = reg%ntr
132  cs%C_p = c_p
133 
134  allocate(cs%id_neutral_diff_tracer_conc_tend(ntr))
135  allocate(cs%id_neutral_diff_tracer_cont_tend(ntr))
136  allocate(cs%id_neutral_diff_tracer_cont_tend_2d(ntr))
137  allocate(cs%id_neutral_diff_tracer_trans_x_2d(ntr))
138  allocate(cs%id_neutral_diff_tracer_trans_y_2d(ntr))
139  cs%id_neutral_diff_tracer_conc_tend(:) = -1
140  cs%id_neutral_diff_tracer_cont_tend(:) = -1
141  cs%id_neutral_diff_tracer_cont_tend_2d(:) = -1
142  cs%id_neutral_diff_tracer_trans_x_2d(:) = -1
143  cs%id_neutral_diff_tracer_trans_y_2d(:) = -1
144 
145  do n=1,ntr
146 
147  if(trim(reg%Tr(n)%name) == 'T') then
148 
149  cs%id_neutral_diff_tracer_conc_tend(n) = register_diag_field('ocean_model', &
150  'ndiff_tracer_conc_tendency_'//trim(reg%Tr(n)%name), diag%axesTL, time, &
151  'Neutral diffusion tracer concentration tendency for '//trim(reg%Tr(n)%name),&
152  'Degree C per second')
153 
154  cs%id_neutral_diff_tracer_cont_tend(n) = register_diag_field('ocean_model', &
155  'ndiff_tracer_cont_tendency_'//trim(reg%Tr(n)%name), diag%axesTL, time, &
156  'Neutral diffusion tracer content tendency for '//trim(reg%Tr(n)%name), &
157  'Watts/m2',cmor_field_name='opottemppmdiff', cmor_units='W m-2', &
158  cmor_standard_name= &
159  'tendency_of_sea_water_potential_temperature_expressed_as_heat_content_due_to_parameterized_mesocale_diffusion', &
160  cmor_long_name = &
161  'Tendency of sea water potential temperature expressed as heat content due to parameterized mesocale diffusion', &
162  v_extensive=.true.)
163 
164  cs%id_neutral_diff_tracer_cont_tend_2d(n) = register_diag_field('ocean_model', &
165  'ndiff_tracer_cont_tendency_2d_'//trim(reg%Tr(n)%name), diag%axesT1, time, &
166  'Depth integrated neutral diffusion tracer content tendency for '//trim(reg%Tr(n)%name), &
167  'Watts/m2',cmor_field_name='opottemppmdiff_2d', cmor_units='W m-2', &
168  cmor_standard_name= &
169  'tendency_of_sea_water_potential_temperature_expressed_as_heat_content_due_to_parameterized_mesocale_diffusion_depth_integrated',&
170  cmor_long_name = &
171  'Tendency of sea water potential temperature expressed as heat content due to parameterized mesocale diffusion depth integrated')
172 
173  cs%id_neutral_diff_tracer_trans_x_2d(n) = register_diag_field('ocean_model', &
174  'ndiff_tracer_trans_x_2d_'//trim(reg%Tr(n)%name), diag%axesCu1, time, &
175  'Depth integrated neutral diffusion zonal tracer transport for '//trim(reg%Tr(n)%name),&
176  'Watts')
177 
178  cs%id_neutral_diff_tracer_trans_y_2d(n) = register_diag_field('ocean_model', &
179  'ndiff_tracer_trans_y_2d_'//trim(reg%Tr(n)%name), diag%axesCv1, time, &
180  'Depth integrated neutral diffusion merid tracer transport for '//trim(reg%Tr(n)%name),&
181  'Watts')
182 
183  elseif(trim(reg%Tr(n)%name) == 'S') then
184 
185  cs%id_neutral_diff_tracer_conc_tend(n) = register_diag_field('ocean_model', &
186  'ndiff_tracer_conc_tendency_'//trim(reg%Tr(n)%name), diag%axesTL, time, &
187  'Neutral diffusion tracer concentration tendency for '//trim(reg%Tr(n)%name),&
188  'tracer concentration units per second')
189 
190  cs%id_neutral_diff_tracer_cont_tend(n) = register_diag_field('ocean_model', &
191  'ndiff_tracer_cont_tendency_'//trim(reg%Tr(n)%name), diag%axesTL, time, &
192  'Neutral diffusion tracer content tendency for '//trim(reg%Tr(n)%name), &
193  'kg m-2 s-1',cmor_field_name='osaltpmdiff', cmor_units='kg m-2 s-1', &
194  cmor_standard_name= &
195  'tendency_of_sea_water_salinity_expressed_as_salt_content_due_to_parameterized_mesocale_diffusion', &
196  cmor_long_name = &
197  'Tendency of sea water salinity expressed as salt content due to parameterized mesocale diffusion', &
198  v_extensive=.true.)
199 
200  cs%id_neutral_diff_tracer_cont_tend_2d(n) = register_diag_field('ocean_model', &
201  'ndiff_tracer_cont_tendency_2d_'//trim(reg%Tr(n)%name), diag%axesT1, time, &
202  'Depth integrated neutral diffusion tracer content tendency for '//trim(reg%Tr(n)%name), &
203  'kg m-2 s-1',cmor_field_name='osaltpmdiff_2d', cmor_units='kg m-2 s-1', &
204  cmor_standard_name= &
205  'tendency_of_sea_water_salinity_expressed_as_salt_content_due_to_parameterized_mesocale_diffusion_depth_integrated',&
206  cmor_long_name = &
207  'Tendency of sea water salinity expressed as salt content due to parameterized mesocale diffusion depth integrated')
208 
209  cs%id_neutral_diff_tracer_trans_x_2d(n) = register_diag_field('ocean_model', &
210  'ndiff_tracer_trans_x_2d_'//trim(reg%Tr(n)%name), diag%axesCu1, time, &
211  'Depth integrated neutral diffusion zonal tracer transport for '//trim(reg%Tr(n)%name),&
212  'kg/s')
213 
214  cs%id_neutral_diff_tracer_trans_y_2d(n) = register_diag_field('ocean_model', &
215  'ndiff_tracer_trans_y_2d_'//trim(reg%Tr(n)%name), diag%axesCv1, time, &
216  'Depth integrated neutral diffusion merid tracer transport for '//trim(reg%Tr(n)%name),&
217  'kg/s')
218 
219  else
220 
221  cs%id_neutral_diff_tracer_conc_tend(n) = register_diag_field('ocean_model', &
222  'ndiff_tracer_conc_tendency_'//trim(reg%Tr(n)%name), diag%axesTL, time, &
223  'Neutral diffusion tracer concentration tendency for '//trim(reg%Tr(n)%name),&
224  'tracer concentration * m-2 s-1')
225 
226  cs%id_neutral_diff_tracer_cont_tend(n) = register_diag_field('ocean_model',&
227  'ndiff_tracer_cont_tendency_'//trim(reg%Tr(n)%name), diag%axesTL, time, &
228  'Neutral diffusion tracer content tendency for '//trim(reg%Tr(n)%name), &
229  'tracer content * m-2 s-1', v_extensive=.true.)
230 
231  cs%id_neutral_diff_tracer_cont_tend_2d(n) = register_diag_field('ocean_model', &
232  'ndiff_tracer_cont_tendency_2d_'//trim(reg%Tr(n)%name), diag%axesTL, time, &
233  'Depth integrated neutral diffusion tracer content tendency for '//trim(reg%Tr(n)%name),&
234  'tracer content * m-2 s-1')
235 
236  cs%id_neutral_diff_tracer_trans_x_2d(n) = register_diag_field('ocean_model', &
237  'ndiff_tracer_trans_x_2d_'//trim(reg%Tr(n)%name), diag%axesCu1, time, &
238  'Depth integrated neutral diffusion zonal tracer transport for '//trim(reg%Tr(n)%name),&
239  'kg/s')
240 
241  cs%id_neutral_diff_tracer_trans_y_2d(n) = register_diag_field('ocean_model', &
242  'ndiff_tracer_trans_y_2d_'//trim(reg%Tr(n)%name), diag%axesCv1, time, &
243  'Depth integrated neutral diffusion merid tracer transport for '//trim(reg%Tr(n)%name),&
244  'kg/s')
245 
246  endif
247 
248  enddo
249 
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:19
Here is the call graph for this function:
Here is the caller graph for this function:

◆ neutral_diffusion_end()

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

Deallocates neutral_diffusion control structure.

Definition at line 1513 of file MOM_neutral_diffusion.F90.

1513  type(neutral_diffusion_cs), pointer :: cs
1514 
1515  if (associated(cs)) deallocate(cs)
1516 

◆ 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(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
csNeutral diffusion control structure

Definition at line 65 of file MOM_neutral_diffusion.F90.

References mdl, and mom_error_handler::mom_error().

65  type(time_type), target, intent(in) :: time !< Time structure
66  type(ocean_grid_type), intent(in) :: g !< Grid structure
67  type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure
68  type(param_file_type), intent(in) :: param_file !< Parameter file structure
69  type(neutral_diffusion_cs), pointer :: cs !< Neutral diffusion control structure
70 
71  ! Local variables
72  character(len=256) :: mesg ! Message for error messages.
73 
74  if (associated(cs)) then
75  call mom_error(fatal, "neutral_diffusion_init called with associated control structure.")
76  return
77  endif
78 
79  ! Log this module and master switch for turning it on/off
80  call log_version(param_file, mdl, version, &
81  "This module implements neutral diffusion of tracers")
82  call get_param(param_file, mdl, "USE_NEUTRAL_DIFFUSION", neutral_diffusion_init, &
83  "If true, enables the neutral diffusion module.", &
84  default=.false.)
85 
86  if (.not.neutral_diffusion_init) then
87  return
88  endif
89 
90  allocate(cs)
91  cs%diag => diag
92 
93 
94  ! Read all relevant parameters and write them to the model log.
95 ! call openParameterBlock(param_file,'NEUTRAL_DIFF')
96 ! call get_param(param_file, mdl, "KHTR", CS%KhTr, &
97 ! "The background along-isopycnal tracer diffusivity.", &
98 ! units="m2 s-1", default=0.0)
99 ! call closeParameterBlock(param_file)
100 
101  ! U-points
102  allocate(cs%uPoL(g%isd:g%ied,g%jsd:g%jed,2*g%ke+2)); cs%uPoL(g%isc-1:g%iec,g%jsc:g%jec,:) = 0.
103  allocate(cs%uPoR(g%isd:g%ied,g%jsd:g%jed,2*g%ke+2)); cs%uPoR(g%isc-1:g%iec,g%jsc:g%jec,:) = 0.
104  allocate(cs%uKoL(g%isd:g%ied,g%jsd:g%jed,2*g%ke+2)); cs%uKoL(g%isc-1:g%iec,g%jsc:g%jec,:) = 0
105  allocate(cs%uKoR(g%isd:g%ied,g%jsd:g%jed,2*g%ke+2)); cs%uKoR(g%isc-1:g%iec,g%jsc:g%jec,:) = 0
106  allocate(cs%uHeff(g%isd:g%ied,g%jsd:g%jed,2*g%ke+1)); cs%uHeff(g%isc-1:g%iec,g%jsc:g%jec,:) = 0
107  ! V-points
108  allocate(cs%vPoL(g%isd:g%ied,g%jsd:g%jed,2*g%ke+2)); cs%vPoL(g%isc:g%iec,g%jsc-1:g%jec,:) = 0.
109  allocate(cs%vPoR(g%isd:g%ied,g%jsd:g%jed,2*g%ke+2)); cs%vPoR(g%isc:g%iec,g%jsc-1:g%jec,:) = 0.
110  allocate(cs%vKoL(g%isd:g%ied,g%jsd:g%jed,2*g%ke+2)); cs%vKoL(g%isc:g%iec,g%jsc-1:g%jec,:) = 0
111  allocate(cs%vKoR(g%isd:g%ied,g%jsd:g%jed,2*g%ke+2)); cs%vKoR(g%isc:g%iec,g%jsc-1:g%jec,:) = 0
112  allocate(cs%vHeff(g%isd:g%ied,g%jsd:g%jed,2*g%ke+1)); cs%vHeff(g%isc:g%iec,g%jsc-1:g%jec,:) = 0
113 
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:19
Here is the call graph for this function:

◆ 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]verboseIt true, write results to stdout

Definition at line 1049 of file MOM_neutral_diffusion.F90.

References absolute_positions(), find_neutral_surface_positions(), interface_scalar(), neutral_surface_flux(), test_data1d(), test_fv_diff(), test_fvlsq_slope(), test_ifndp(), and test_nsp().

Referenced by mom_unit_tests::unit_tests().

1049  logical, intent(in) :: verbose !< It true, write results to stdout
1050  ! Local variables
1051  integer, parameter :: nk = 4
1052  real, dimension(nk+1) :: til, tir1, tir2, tir4, tio ! Test interface temperatures
1053  real, dimension(nk) :: tl ! Test layer temperatures
1054  real, dimension(nk+1) :: sil ! Test interface salinities
1055  real, dimension(nk+1) :: pil, pir4 ! Test interface positions
1056  real, dimension(2*nk+2) :: pilro, pirlo ! Test positions
1057  integer, dimension(2*nk+2) :: kol, kor ! Test indexes
1058  real, dimension(2*nk+1) :: heff ! Test positions
1059  real, dimension(2*nk+1) :: flx ! Test flux
1060  integer :: k
1061  logical :: v
1062 
1063  v = verbose
1064 
1065  neutral_diffusion_unit_tests = .false. ! Normally return false
1066  write(*,*) '==== MOM_neutral_diffusion: neutral_diffusion_unit_tests ='
1067 
1068  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fv_diff(v,1.,1.,1., 0.,1.,2., 1., 'FV: Straight line on uniform grid')
1069  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fv_diff(v,1.,1.,0., 0.,4.,8., 7., 'FV: Vanished right cell')
1070  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fv_diff(v,0.,1.,1., 0.,4.,8., 7., 'FV: Vanished left cell')
1071  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fv_diff(v,1.,2.,4., 0.,3.,9., 4., 'FV: Stretched grid')
1072  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fv_diff(v,2.,0.,2., 0.,1.,2., 0., 'FV: Vanished middle cell')
1073  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fv_diff(v,0.,1.,0., 0.,1.,2., 2., 'FV: Vanished on both sides')
1074  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fv_diff(v,1.,0.,0., 0.,1.,2., 0., 'FV: Two vanished cell sides')
1075  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fv_diff(v,0.,0.,0., 0.,1.,2., 0., 'FV: All vanished cells')
1076 
1077  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fvlsq_slope(v,1.,1.,1., 0.,1.,2., 1., 'LSQ: Straight line on uniform grid')
1078  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fvlsq_slope(v,1.,1.,0., 0.,1.,2., 1., 'LSQ: Vanished right cell')
1079  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fvlsq_slope(v,0.,1.,1., 0.,1.,2., 1., 'LSQ: Vanished left cell')
1080  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fvlsq_slope(v,1.,2.,4., 0.,3.,9., 2., 'LSQ: Stretched grid')
1081  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fvlsq_slope(v,1.,0.,1., 0.,1.,2., 2., 'LSQ: Vanished middle cell')
1082  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fvlsq_slope(v,0.,1.,0., 0.,1.,2., 0., 'LSQ: Vanished on both sides')
1083  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fvlsq_slope(v,1.,0.,0., 0.,1.,2., 0., 'LSQ: Two vanished cell sides')
1084  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_fvlsq_slope(v,0.,0.,0., 0.,1.,2., 0., 'LSQ: All vanished cells')
1085 
1086  call interface_scalar(4, (/10.,10.,10.,10./), (/24.,18.,12.,6./), tio, 1)
1087  !neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_data1d(5, Tio, (/27.,21.,15.,9.,3./), 'Linear profile, interface temperatures')
1088  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_data1d(v,5, tio, (/24.,22.5,15.,7.5,6./), 'Linear profile, linear interface temperatures')
1089  call interface_scalar(4, (/10.,10.,10.,10./), (/24.,18.,12.,6./), tio, 2)
1090  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_data1d(v,5, tio, (/24.,22.,15.,8.,6./), 'Linear profile, PPM interface temperatures')
1091 
1092  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(v,-1.0, 0., 1.0, 1.0, 0.5, 'Check mid-point')
1093  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(v, 0.0, 0., 1.0, 1.0, 0.0, 'Check bottom')
1094  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(v, 0.1, 0., 1.1, 1.0, 0.0, 'Check below')
1095  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(v,-1.0, 0., 0.0, 1.0, 1.0, 'Check top')
1096  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(v,-1.0, 0., -0.1, 1.0, 1.0, 'Check above')
1097  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(v,-1.0, 0., 3.0, 1.0, 0.25, 'Check 1/4')
1098  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(v,-3.0, 0., 1.0, 1.0, 0.75, 'Check 3/4')
1099  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(v, 1.0, 0., 1.0, 1.0, 0.0, 'Check dRho=0 below')
1100  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(v,-1.0, 0., -1.0, 1.0, 1.0, 'Check dRho=0 above')
1101  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(v, 0.0, 0., 0.0, 1.0, 0.5, 'Check dRho=0 mid')
1102  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_ifndp(v,-2.0, .5, 5.0, 0.5, 0.5, 'Check dP=0')
1103 
1104  ! Identical columns
1105  call find_neutral_surface_positions(3, &
1106  (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Left positions, T and S
1107  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
1108  (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Right positions, T and S
1109  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
1110  pilro, pirlo, kol, kor, heff)
1111  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(v, 3, kol, kor, pilro, pirlo, heff, &
1112  (/1,1,2,2,3,3,3,3/), & ! KoL
1113  (/1,1,2,2,3,3,3,3/), & ! KoR
1114  (/0.,0.,0.,0.,0.,0.,1.,1./), & ! pL
1115  (/0.,0.,0.,0.,0.,0.,1.,1./), & ! pR
1116  (/0.,10.,0.,10.,0.,10.,0./), & ! hEff
1117  'Identical columns')
1118  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_data1d(v, 8, &
1119  absolute_positions(3, (/0.,10.,20.,30./), kol, pilro), &
1120  (/0.,0.,10.,10.,20.,20.,30.,30./), '... left positions')
1121  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_data1d(v, 8, &
1122  absolute_positions(3, (/0.,10.,20.,30./), kor, pirlo), &
1123  (/0.,0.,10.,10.,20.,20.,30.,30./), '... right positions')
1124  call neutral_surface_flux(3, (/10.,10.,10./), (/10.,10.,10./), & ! nk, hL, hR
1125  (/20.,16.,12./), (/20.,16.,12./), & ! Tl, Tr
1126  pilro, pirlo, kol, kor, heff, flx)
1127  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_data1d(v, 7, flx, &
1128  (/0.,0.,0.,0.,0.,0.,0./), 'Identical columns, rho flux (=0)')
1129  call neutral_surface_flux(3, (/10.,10.,10./), (/10.,10.,10./), & ! nk, hL, hR
1130  (/-1.,-1.,-1./), (/1.,1.,1./), & ! Sl, Sr
1131  pilro, pirlo, kol, kor, heff, flx)
1132  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_data1d(v, 7, flx, &
1133  (/0.,20.,0.,20.,0.,20.,0./), 'Identical columns, S flux')
1134 
1135  ! Right column slightly cooler than left
1136  call find_neutral_surface_positions(3, &
1137  (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Left positions, T and S
1138  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
1139  (/0.,10.,20.,30./), (/20.,16.,12.,8./), (/0.,0.,0.,0./), & ! Right positions, T and S
1140  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
1141  pilro, pirlo, kol, kor, heff)
1142  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(v, 3, kol, kor, pilro, pirlo, heff, &
1143  (/1,1,2,2,3,3,3,3/), & ! kL
1144  (/1,1,1,2,2,3,3,3/), & ! kR
1145  (/0.,0.5,0.,0.5,0.,0.5,1.,1./), & ! pL
1146  (/0.,0.,0.5,0.,0.5,0.,0.5,1./), & ! pR
1147  (/0.,5.,5.,5.,5.,5.,0./), & ! hEff
1148  'Right column slightly cooler')
1149  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_data1d(v, 8, &
1150  absolute_positions(3, (/0.,10.,20.,30./), kol, pilro), &
1151  (/0.,5.,10.,15.,20.,25.,30.,30./), '... left positions')
1152  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_data1d(v, 8, &
1153  absolute_positions(3, (/0.,10.,20.,30./), kor, pirlo), &
1154  (/0.,0.,5.,10.,15.,20.,25.,30./), '... right positions')
1155 
1156  ! Right column slightly warmer than left
1157  call find_neutral_surface_positions(3, &
1158  (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Left positions, T and S
1159  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
1160  (/0.,10.,20.,30./), (/24.,20.,16.,12./), (/0.,0.,0.,0./), & ! Right positions, T and S
1161  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
1162  pilro, pirlo, kol, kor, heff)
1163  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(v, 3, kol, kor, pilro, pirlo, heff, &
1164  (/1,1,1,2,2,3,3,3/), & ! kL
1165  (/1,1,2,2,3,3,3,3/), & ! kR
1166  (/0.,0.,0.5,0.,0.5,0.,0.5,1./), & ! pL
1167  (/0.,0.5,0.,0.5,0.,0.5,1.,1./), & ! pR
1168  (/0.,5.,5.,5.,5.,5.,0./), & ! hEff
1169  'Right column slightly warmer')
1170 
1171  ! Right column somewhat cooler than left
1172  call find_neutral_surface_positions(3, &
1173  (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Left positions, T and S
1174  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
1175  (/0.,10.,20.,30./), (/16.,12.,8.,4./), (/0.,0.,0.,0./), & ! Right positions, T and S
1176  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
1177  pilro, pirlo, kol, kor, heff)
1178  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(v, 3, kol, kor, pilro, pirlo, heff, &
1179  (/1,2,2,3,3,3,3,3/), & ! kL
1180  (/1,1,1,1,2,2,3,3/), & ! kR
1181  (/0.,0.,0.5,0.,0.5,1.,1.,1./), & ! pL
1182  (/0.,0.,0.,0.5,0.,0.5,0.,1./), & ! pR
1183  (/0.,0.,5.,5.,5.,0.,0./), & ! hEff
1184  'Right column somewhat cooler')
1185 
1186  ! Right column much colder than left with no overlap
1187  call find_neutral_surface_positions(3, &
1188  (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Left positions, T and S
1189  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
1190  (/0.,10.,20.,30./), (/9.,7.,5.,3./), (/0.,0.,0.,0./), & ! Right positions, T and S
1191  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
1192  pilro, pirlo, kol, kor, heff)
1193  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(v, 3, kol, kor, pilro, pirlo, heff, &
1194  (/1,2,3,3,3,3,3,3/), & ! kL
1195  (/1,1,1,1,1,2,3,3/), & ! kR
1196  (/0.,0.,0.,1.,1.,1.,1.,1./), & ! pL
1197  (/0.,0.,0.,0.,0.,0.,0.,1./), & ! pR
1198  (/0.,0.,0.,0.,0.,0.,0./), & ! hEff
1199  'Right column much cooler')
1200 
1201  ! Right column with mixed layer
1202  call find_neutral_surface_positions(3, &
1203  (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Left positions, T and S
1204  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
1205  (/0.,10.,20.,30./), (/14.,14.,10.,2./), (/0.,0.,0.,0./), & ! Right positions, T and S
1206  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
1207  pilro, pirlo, kol, kor, heff)
1208  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(v, 3, kol, kor, pilro, pirlo, heff, &
1209  (/1,2,3,3,3,3,3,3/), & ! kL
1210  (/1,1,1,1,2,3,3,3/), & ! kR
1211  (/0.,0.,0.,0.,0.,1.,1.,1./), & ! pL
1212  (/0.,0.,0.,0.,0.,0.,0.,1./), & ! pR
1213  (/0.,0.,0.,0.,10.,0.,0./), & ! hEff
1214  'Right column with mixed layer')
1215 
1216  ! Identical columns with mixed layer
1217  call find_neutral_surface_positions(3, &
1218  (/0.,10.,20.,30./), (/14.,14.,10.,2./), (/0.,0.,0.,0./), & ! Left positions, T and S
1219  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
1220  (/0.,10.,20.,30./), (/14.,14.,10.,2./), (/0.,0.,0.,0./), & ! Right positions, T and S
1221  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
1222  pilro, pirlo, kol, kor, heff)
1223  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(v, 3, kol, kor, pilro, pirlo, heff, &
1224  (/1,1,2,2,3,3,3,3/), & ! kL
1225  (/1,1,2,2,3,3,3,3/), & ! kR
1226  (/0.,0.,0.,0.,0.,0.,1.,1./), & ! pL
1227  (/0.,0.,0.,0.,0.,0.,1.,1./), & ! pR
1228  (/0.,10.,0.,10.,0.,10.,0./), & ! hEff
1229  'Indentical columns with mixed layer')
1230 
1231  ! Right column with unstable mixed layer
1232  call find_neutral_surface_positions(3, &
1233  (/0.,10.,20.,30./), (/14.,14.,10.,2./), (/0.,0.,0.,0./), & ! Left positions, T and S
1234  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
1235  (/0.,10.,20.,30./), (/10.,14.,12.,4./), (/0.,0.,0.,0./), & ! Right positions, T and S
1236  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
1237  pilro, pirlo, kol, kor, heff)
1238  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(v, 3, kol, kor, pilro, pirlo, heff, &
1239  (/1,2,3,3,3,3,3,3/), & ! kL
1240  (/1,1,1,2,3,3,3,3/), & ! kR
1241  (/0.,0.,0.,0.,0.,0.,.75,1./), & ! pL
1242  (/0.,0.,0.,0.,0.,0.25,1.,1./), & ! pR
1243  (/0.,0.,0.,0.,0.,7.5,0./), & ! hEff
1244  'Right column with unstable mixed layer')
1245 
1246  ! Left column with unstable mixed layer
1247  call find_neutral_surface_positions(3, &
1248  (/0.,10.,20.,30./), (/10.,14.,12.,4./), (/0.,0.,0.,0./), & ! Left positions, T and S
1249  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
1250  (/0.,10.,20.,30./), (/14.,14.,10.,2./), (/0.,0.,0.,0./), & ! Right positions, T and S
1251  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
1252  pilro, pirlo, kol, kor, heff)
1253  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(v, 3, kol, kor, pilro, pirlo, heff, &
1254  (/1,1,1,2,3,3,3,3/), & ! kL
1255  (/1,2,3,3,3,3,3,3/), & ! kR
1256  (/0.,0.,0.,0.,0.,0.25,1.,1./), & ! pL
1257  (/0.,0.,0.,0.,0.,0.,.75,1./), & ! pR
1258  (/0.,0.,0.,0.,0.,7.5,0./), & ! hEff
1259  'Left column with unstable mixed layer')
1260 
1261  ! Two unstable mixed layers
1262  call find_neutral_surface_positions(3, &
1263  (/0.,10.,20.,30./), (/8.,12.,10.,2./), (/0.,0.,0.,0./), & ! Left positions, T and S
1264  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
1265  (/0.,10.,20.,30./), (/10.,14.,12.,4./), (/0.,0.,0.,0./), & ! Right positions, T and S
1266  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
1267  pilro, pirlo, kol, kor, heff)
1268  neutral_diffusion_unit_tests = neutral_diffusion_unit_tests .or. test_nsp(v, 3, kol, kor, pilro, pirlo, heff, &
1269  (/1,1,1,1,2,3,3,3/), & ! kL
1270  (/1,2,3,3,3,3,3,3/), & ! kR
1271  (/0.,0.,0.,0.,0.,0.,0.75,1./), & ! pL
1272  (/0.,0.,0.,0.5,0.5,0.5,1.,1./), & ! pR
1273  (/0.,0.,0.,0.,0.,6.,0./), & ! hEff
1274  'Two unstable mixed layers')
1275 
1276  if (.not. neutral_diffusion_unit_tests) write(*,*) 'Pass'
1277 
Here is the call graph for this function:
Here is the caller graph for this function:

◆ neutral_surface_flux()

subroutine mom_neutral_diffusion::neutral_surface_flux ( integer, intent(in)  nk,
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(2*nk+2), intent(in)  PiL,
real, dimension(2*nk+2), intent(in)  PiR,
integer, dimension(2*nk+2), intent(in)  KoL,
integer, dimension(2*nk+2), intent(in)  KoR,
real, dimension(2*nk+1), intent(in)  hEff,
real, dimension(2*nk+1), intent(inout)  Flx 
)
private

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

Parameters
[in]nkNumber of levels
[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)

Definition at line 958 of file MOM_neutral_diffusion.F90.

References interface_scalar(), ppm_ave(), and signum().

Referenced by neutral_diffusion(), and neutral_diffusion_unit_tests().

958  integer, intent(in) :: nk !< Number of levels
959  real, dimension(nk), intent(in) :: hl !< Left-column layer thickness (Pa)
960  real, dimension(nk), intent(in) :: hr !< Right-column layer thickness (Pa)
961  real, dimension(nk), intent(in) :: tl !< Left-column layer tracer (conc, e.g. degC)
962  real, dimension(nk), intent(in) :: tr !< Right-column layer tracer (conc, e.g. degC)
963  real, dimension(2*nk+2), intent(in) :: pil !< Fractional position of neutral surface within layer KoL of left column
964  real, dimension(2*nk+2), intent(in) :: pir !< Fractional position of neutral surface within layer KoR of right column
965  integer, dimension(2*nk+2), intent(in) :: kol !< Index of first left interface above neutral surface
966  integer, dimension(2*nk+2), intent(in) :: kor !< Index of first right interface above neutral surface
967  real, dimension(2*nk+1), intent(in) :: heff !< Effective thickness between two neutral surfaces (Pa)
968  real, dimension(2*nk+1), intent(inout) :: flx !< Flux of tracer between pairs of neutral layers (conc H)
969 
970  ! Local variables
971  integer :: k_sublayer, klb, klt, krb, krt, k
972  real :: t_right_top, t_right_bottom, t_right_layer
973  real :: t_left_top, t_left_bottom, t_left_layer
974  real :: dt_top, dt_bottom, dt_layer, dt_ave
975  real, dimension(nk+1) :: til !< Left-column interface tracer (conc, e.g. degC)
976  real, dimension(nk+1) :: tir !< Right-column interface tracer (conc, e.g. degC)
977  real, dimension(nk) :: al_l !< Left-column left edge value of tracer (conc, e.g. degC)
978  real, dimension(nk) :: ar_l !< Left-column right edge value of tracer (conc, e.g. degC)
979  real, dimension(nk) :: al_r !< Right-column left edge value of tracer (conc, e.g. degC)
980  real, dimension(nk) :: ar_r !< Right-column right edge value of tracer (conc, e.g. degC)
981 
982  call interface_scalar(nk, hl, tl, til, 2)
983  call interface_scalar(nk, hr, tr, tir, 2)
984 
985  ! Setup reconstruction edge values
986  do k = 1, nk
987  al_l(k) = til(k)
988  ar_l(k) = til(k+1)
989  if ( signum(1., ar_l(k) - tl(k))*signum(1., tl(k) - al_l(k)) <= 0.0 ) then
990  al_l(k) = tl(k)
991  ar_l(k) = tl(k)
992  elseif ( sign(3., ar_l(k) - al_l(k)) * ( (tl(k) - al_l(k)) + (tl(k) - ar_l(k))) > abs(ar_l(k) - al_l(k)) ) then
993  al_l(k) = tl(k) + 2.0 * ( tl(k) - ar_l(k) )
994  elseif ( sign(3., ar_l(k) - al_l(k)) * ( (tl(k) - al_l(k)) + (tl(k) - ar_l(k))) < -abs(ar_l(k) - al_l(k)) ) then
995  ar_l(k) = tl(k) + 2.0 * ( tl(k) - al_l(k) )
996  endif
997  al_r(k) = tir(k)
998  ar_r(k) = tir(k+1)
999  if ( signum(1., ar_r(k) - tr(k))*signum(1., tr(k) - al_r(k)) <= 0.0 ) then
1000  al_r(k) = tr(k)
1001  ar_r(k) = tr(k)
1002  elseif ( sign(3., ar_r(k) - al_r(k)) * ( (tr(k) - al_r(k)) + (tr(k) - ar_r(k))) > abs(ar_r(k) - al_r(k)) ) then
1003  al_r(k) = tr(k) + 2.0 * ( tr(k) - ar_r(k) )
1004  elseif ( sign(3., ar_r(k) - al_r(k)) * ( (tr(k) - al_r(k)) + (tr(k) - ar_r(k))) < -abs(ar_r(k) - al_r(k)) ) then
1005  ar_r(k) = tr(k) + 2.0 * ( tr(k) - al_r(k) )
1006  endif
1007  enddo
1008 
1009  do k_sublayer = 1, 2*nk+1
1010  if (heff(k_sublayer) == 0.) then
1011  flx(k_sublayer) = 0.
1012  else
1013 
1014  klb = kol(k_sublayer+1)
1015  t_left_bottom = ( 1. - pil(k_sublayer+1) ) * til(klb) + pil(k_sublayer+1) * til(klb+1)
1016 
1017  klt = kol(k_sublayer)
1018  t_left_top = ( 1. - pil(k_sublayer) ) * til(klt) + pil(k_sublayer) * til(klt+1)
1019 
1020  t_left_layer = ppm_ave(pil(k_sublayer), pil(k_sublayer+1) + real(klb-klt), &
1021  al_l(klt), ar_l(klt), tl(klt))
1022 
1023  krb = kor(k_sublayer+1)
1024  t_right_bottom = ( 1. - pir(k_sublayer+1) ) * tir(krb) + pir(k_sublayer+1) * tir(krb+1)
1025 
1026  krt = kor(k_sublayer)
1027  t_right_top = ( 1. - pir(k_sublayer) ) * tir(krt) + pir(k_sublayer) * tir(krt+1)
1028 
1029  t_right_layer = ppm_ave(pir(k_sublayer), pir(k_sublayer+1) + real(krb-krt), &
1030  al_r(krt), ar_r(krt), tr(krt))
1031 
1032  dt_top = t_right_top - t_left_top
1033  dt_bottom = t_right_bottom - t_left_bottom
1034  dt_ave = 0.5 * ( dt_top + dt_bottom )
1035  dt_layer = t_right_layer - t_left_layer
1036  if (signum(1.,dt_top) * signum(1.,dt_bottom) <= 0. .or. signum(1.,dt_ave) * signum(1.,dt_layer) <= 0. ) then
1037  dt_ave = 0.
1038  else
1039  dt_ave = dt_layer
1040  endif
1041  flx(k_sublayer) = dt_ave * heff(k_sublayer)
1042  endif
1043  enddo
1044 
Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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 units)
[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 593 of file MOM_neutral_diffusion.F90.

References fv_diff(), fvlsq_slope(), and signum().

Referenced by interface_scalar().

593  integer, intent(in) :: nk !< Number of levels
594  real, dimension(nk), intent(in) :: h !< Layer thickness (H units)
595  real, dimension(nk), intent(in) :: s !< Layer salinity (conc, e.g. ppt)
596  integer, intent(in) :: c_method !< Method to use for the centered difference
597  integer, intent(in) :: b_method !< =1, use PCM in first/last cell, =2 uses linear extrapolation
598  real, dimension(nk), intent(inout) :: diff !< Scalar difference across layer (conc, e.g. ppt)
599  !! determined by the following values for c_method:
600  !! 1. Second order finite difference (not recommended)
601  !! 2. Second order finite volume (used in original PPM)
602  !! 3. Finite-volume weighted least squares linear fit
603  !! \todo The use of c_method to choose a scheme is inefficient
604  !! and should eventually be moved up the call tree.
605 
606  ! Local variables
607  integer :: k
608  real :: hkm1, hk, hkp1, skm1, sk, skp1, diff_l, diff_r, diff_c
609 
610  do k = 2, nk-1
611  hkm1 = h(k-1)
612  hk = h(k)
613  hkp1 = h(k+1)
614 
615  if ( ( hkp1 + hk ) * ( hkm1 + hk ) > 0.) then
616  skm1 = s(k-1)
617  sk = s(k)
618  skp1 = s(k+1)
619  if (c_method==1) then
620  ! Simple centered diff (from White)
621  if ( hk + 0.5 * (hkm1 + hkp1) /= 0. ) then
622  diff_c = ( skp1 - skm1 ) * ( hk / ( hk + 0.5 * (hkm1 + hkp1) ) )
623  else
624  diff_c = 0.
625  endif
626  elseif (c_method==2) then
627  ! Second order accurate centered FV slope (from Colella and Woodward, JCP 1984)
628  diff_c = fv_diff(hkm1, hk, hkp1, skm1, sk, skp1)
629  elseif (c_method==3) then
630  ! Second order accurate finite-volume least squares slope
631  diff_c = hk * fvlsq_slope(hkm1, hk, hkp1, skm1, sk, skp1)
632  endif
633  ! Limit centered slope by twice the side differenced slopes
634  diff_l = 2. * ( sk - skm1 )
635  diff_r = 2. * ( skp1 - sk )
636  if ( signum(1., diff_l) * signum(1., diff_r) <= 0. ) then
637  diff(k) = 0. ! PCM for local extrema
638  else
639  diff(k) = sign( min( abs(diff_l), abs(diff_c), abs(diff_r) ), diff_c )
640  endif
641  else
642  diff(k) = 0. ! PCM next to vanished layers
643  endif
644  enddo
645  if (b_method==1) then ! PCM for top and bottom layer
646  diff(1) = 0.
647  diff(nk) = 0.
648  elseif (b_method==2) then ! Linear extrapolation for top and bottom interfaces
649  diff(1) = ( s(2) - s(1) ) * 2. * ( h(1) / ( h(1) + h(2) ) )
650  diff(nk) = s(nk) - s(nk-1) * 2. * ( h(nk) / ( h(nk-1) + h(nk) ) )
651  endif
652 
Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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 554 of file MOM_neutral_diffusion.F90.

Referenced by neutral_surface_flux().

554  real, intent(in) :: xl !< Fraction position of left bound (0,1)
555  real, intent(in) :: xr !< Fraction position of right bound (0,1)
556  real, intent(in) :: al !< Left edge scalar value, at x=0
557  real, intent(in) :: ar !< Right edge scalar value, at x=1
558  real, intent(in) :: amean !< Average scalar value of cell
559 
560  ! Local variables
561  real :: dx, xave, a6, a6o3
562 
563  dx = xr - xl
564  xave = 0.5 * ( xr + xl )
565  a6o3 = 2. * amean - ( al + ar ) ! a6 / 3.
566  a6 = 3. * a6o3
567 
568  if (dx<0.) then
569  stop 'ppm_ave: dx<0 should not happend!'
570  elseif (dx>1.) then
571  stop 'ppm_ave: dx>1 should not happend!'
572  elseif (dx==0.) then
573  ppm_ave = al + ( ar - al ) * xr + a6 * xr * ( 1. - xr )
574  else
575  ppm_ave = ( al + xave * ( ( ar - al ) + a6 ) ) - a6o3 * ( xr**2 + xr * xl + xl**2 )
576  endif
577 
Here is the caller graph for this function:

◆ 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 
)
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

Definition at line 514 of file MOM_neutral_diffusion.F90.

Referenced by interface_scalar().

514  real, intent(in) :: hkm1 !< Width of cell k-1
515  real, intent(in) :: hk !< Width of cell k
516  real, intent(in) :: hkp1 !< Width of cell k+1
517  real, intent(in) :: hkp2 !< Width of cell k+2
518  real, intent(in) :: ak !< Average scalar value of cell k
519  real, intent(in) :: akp1 !< Average scalar value of cell k+1
520  real, intent(in) :: pk !< PLM slope for cell k
521  real, intent(in) :: pkp1 !< PLM slope for cell k+1
522 
523  ! Local variables
524  real :: r_hk_hkp1, r_2hk_hkp1, r_hk_2hkp1, f1, f2, f3, f4
525  real, parameter :: h_neglect = 1.e-30
526 
527  r_hk_hkp1 = hk + hkp1
528  if (r_hk_hkp1 <= 0.) then
529  ppm_edge = 0.5 * ( ak + akp1 )
530  return
531  endif
532  r_hk_hkp1 = 1. / r_hk_hkp1
533  if (hk<hkp1) then
534  ppm_edge = ak + ( hk * r_hk_hkp1 ) * ( akp1 - ak )
535  else
536  ppm_edge = akp1 + ( hkp1 * r_hk_hkp1 ) * ( ak - akp1 )
537  endif
538 
539  r_2hk_hkp1 = 1. / ( ( 2. * hk + hkp1 ) + h_neglect )
540  r_hk_2hkp1 = 1. / ( ( hk + 2. * hkp1 ) + h_neglect )
541  f1 = 1./ ( ( hk + hkp1) + ( hkm1 + hkp2 ) )
542  f2 = 2. * ( hkp1 * hk ) * r_hk_hkp1 * &
543  ( ( hkm1 + hk ) * r_2hk_hkp1 - ( hkp2 + hkp1 ) * r_hk_2hkp1 )
544  f3 = hk * ( hkm1 + hk ) * r_2hk_hkp1
545  f4 = hkp1 * ( hkp1 + hkp2 ) * r_hk_2hkp1
546 
547  ppm_edge = ppm_edge + f1 * ( f2 * ( akp1 - ak ) - ( f3 * pkp1 - f4 * pk ) )
548 
Here is the caller graph for this function:

◆ 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 582 of file MOM_neutral_diffusion.F90.

Referenced by neutral_surface_flux(), and plm_diff().

582  real, intent(in) :: a !< The magnitude argument
583  real, intent(in) :: x !< The sign (or zero) argument
584 
585  signum = sign(a,x)
586  if (x==0.) signum = 0.
587 
Here is the caller graph for this function:

◆ 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 1376 of file MOM_neutral_diffusion.F90.

Referenced by neutral_diffusion_unit_tests().

1376  logical, intent(in) :: verbose !< If true, write results to stdout
1377  integer, intent(in) :: nk !< Number of layers
1378  real, dimension(nk), intent(in) :: po !< Calculated answer
1379  real, dimension(nk), intent(in) :: ptrue !< True answer
1380  character(len=*), intent(in) :: title !< Title for messages
1381 
1382  ! Local variables
1383  integer :: k, stdunit
1384 
1385  test_data1d = .false.
1386  do k = 1,nk
1387  if (po(k) /= ptrue(k)) test_data1d = .true.
1388  enddo
1389 
1390  if (test_data1d .or. verbose) then
1391  stdunit = 6
1392  if (test_data1d) stdunit = 0 ! In case of wrong results, write to error stream
1393  write(stdunit,'(a)') title
1394  do k = 1,nk
1395  if (po(k) /= ptrue(k)) then
1396  test_data1d = .true.
1397  write(stdunit,'(a,i2,2(x,a,f20.16),x,a,1pe22.15,x,a)') 'k=',k,'Po=',po(k),'Ptrue=',ptrue(k),'err=',po(k)-ptrue(k),'WRONG!'
1398  else
1399  if (verbose) &
1400  write(stdunit,'(a,i2,2(x,a,f20.16),x,a,1pe22.15)') 'k=',k,'Po=',po(k),'Ptrue=',ptrue(k),'err=',po(k)-ptrue(k)
1401  endif
1402  enddo
1403  endif
1404 
Here is the caller graph for this function:

◆ 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 1409 of file MOM_neutral_diffusion.F90.

1409  logical, intent(in) :: verbose !< If true, write results to stdout
1410  integer, intent(in) :: nk !< Number of layers
1411  integer, dimension(nk), intent(in) :: po !< Calculated answer
1412  integer, dimension(nk), intent(in) :: ptrue !< True answer
1413  character(len=*), intent(in) :: title !< Title for messages
1414 
1415  ! Local variables
1416  integer :: k, stdunit
1417 
1418  test_data1di = .false.
1419  do k = 1,nk
1420  if (po(k) /= ptrue(k)) test_data1di = .true.
1421  enddo
1422 
1423  if (test_data1di .or. verbose) then
1424  stdunit = 6
1425  if (test_data1di) stdunit = 0 ! In case of wrong results, write to error stream
1426  write(stdunit,'(a)') title
1427  do k = 1,nk
1428  if (po(k) /= ptrue(k)) then
1429  test_data1di = .true.
1430  write(stdunit,'(a,i2,2(x,a,i5),x,a)') 'k=',k,'Io=',po(k),'Itrue=',ptrue(k),'WRONG!'
1431  else
1432  if (verbose) &
1433  write(stdunit,'(a,i2,2(x,a,i5))') 'k=',k,'Io=',po(k),'Itrue=',ptrue(k)
1434  endif
1435  enddo
1436  endif
1437 

◆ 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 1282 of file MOM_neutral_diffusion.F90.

References fv_diff().

Referenced by neutral_diffusion_unit_tests().

1282  logical, intent(in) :: verbose !< If true, write results to stdout
1283  real, intent(in) :: hkm1 !< Left cell width
1284  real, intent(in) :: hk !< Center cell width
1285  real, intent(in) :: hkp1 !< Right cell width
1286  real, intent(in) :: skm1 !< Left cell average value
1287  real, intent(in) :: sk !< Center cell average value
1288  real, intent(in) :: skp1 !< Right cell average value
1289  real, intent(in) :: ptrue !< True answer (Pa)
1290  character(len=*), intent(in) :: title !< Title for messages
1291 
1292  ! Local variables
1293  integer :: stdunit
1294  real :: pret
1295 
1296  pret = fv_diff(hkm1, hk, hkp1, skm1, sk, skp1)
1297  test_fv_diff = (pret /= ptrue)
1298 
1299  if (test_fv_diff .or. verbose) then
1300  stdunit = 6
1301  if (test_fv_diff) stdunit = 0 ! In case of wrong results, write to error stream
1302  write(stdunit,'(a)') title
1303  if (test_fv_diff) then
1304  write(stdunit,'(2(x,a,f20.16),x,a)') 'pRet=',pret,'pTrue=',ptrue,'WRONG!'
1305  else
1306  write(stdunit,'(2(x,a,f20.16))') 'pRet=',pret,'pTrue=',ptrue
1307  endif
1308  endif
1309 
Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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 1314 of file MOM_neutral_diffusion.F90.

References fvlsq_slope().

Referenced by neutral_diffusion_unit_tests().

1314  logical, intent(in) :: verbose !< If true, write results to stdout
1315  real, intent(in) :: hkm1 !< Left cell width
1316  real, intent(in) :: hk !< Center cell width
1317  real, intent(in) :: hkp1 !< Right cell width
1318  real, intent(in) :: skm1 !< Left cell average value
1319  real, intent(in) :: sk !< Center cell average value
1320  real, intent(in) :: skp1 !< Right cell average value
1321  real, intent(in) :: ptrue !< True answer (Pa)
1322  character(len=*), intent(in) :: title !< Title for messages
1323 
1324  ! Local variables
1325  integer :: stdunit
1326  real :: pret
1327 
1328  pret = fvlsq_slope(hkm1, hk, hkp1, skm1, sk, skp1)
1329  test_fvlsq_slope = (pret /= ptrue)
1330 
1331  if (test_fvlsq_slope .or. verbose) then
1332  stdunit = 6
1333  if (test_fvlsq_slope) stdunit = 0 ! In case of wrong results, write to error stream
1334  write(stdunit,'(a)') title
1335  if (test_fvlsq_slope) then
1336  write(stdunit,'(2(x,a,f20.16),x,a)') 'pRet=',pret,'pTrue=',ptrue,'WRONG!'
1337  else
1338  write(stdunit,'(2(x,a,f20.16))') 'pRet=',pret,'pTrue=',ptrue
1339  endif
1340  endif
1341 
Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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/m3)
[in]pnegInterface position of lighter density (Pa)
[in]rhoposHeavier density (kg/m3)
[in]pposInterface position of heavier density (Pa)
[in]ptrueTrue answer (Pa)
[in]titleTitle for messages

Definition at line 1346 of file MOM_neutral_diffusion.F90.

References interpolate_for_nondim_position().

Referenced by neutral_diffusion_unit_tests().

1346  logical, intent(in) :: verbose !< If true, write results to stdout
1347  real, intent(in) :: rhoneg !< Lighter density (kg/m3)
1348  real, intent(in) :: pneg !< Interface position of lighter density (Pa)
1349  real, intent(in) :: rhopos !< Heavier density (kg/m3)
1350  real, intent(in) :: ppos !< Interface position of heavier density (Pa)
1351  real, intent(in) :: ptrue !< True answer (Pa)
1352  character(len=*), intent(in) :: title !< Title for messages
1353 
1354  ! Local variables
1355  integer :: stdunit
1356  real :: pret
1357 
1358  pret = interpolate_for_nondim_position(rhoneg, pneg, rhopos, ppos)
1359  test_ifndp = (pret /= ptrue)
1360 
1361  if (test_ifndp .or. verbose) then
1362  stdunit = 6
1363  if (test_ifndp) stdunit = 0 ! In case of wrong results, write to error stream
1364  write(stdunit,'(a)') title
1365  if (test_ifndp) then
1366  write(stdunit,'(4(x,a,f20.16),2(x,a,1pe22.15),x,a)') 'r1=',rhoneg,'p1=',pneg,'r2=',rhopos,'p2=',ppos,'pRet=',pret,'pTrue=',ptrue,'WRONG!'
1367  else
1368  write(stdunit,'(4(x,a,f20.16),2(x,a,1pe22.15))') 'r1=',rhoneg,'p1=',pneg,'r2=',rhopos,'p2=',ppos,'pRet=',pret,'pTrue=',ptrue
1369  endif
1370  endif
1371 
Here is the call graph for this function:
Here is the caller graph for this function:

◆ test_nsp()

logical function mom_neutral_diffusion::test_nsp ( logical, intent(in)  verbose,
integer, intent(in)  nk,
integer, dimension(2*nk+2), intent(in)  KoL,
integer, dimension(2*nk+2), intent(in)  KoR,
real, dimension(2*nk+2), intent(in)  pL,
real, dimension(2*nk+2), intent(in)  pR,
real, dimension(2*nk+1), intent(in)  hEff,
integer, dimension(2*nk+2), intent(in)  KoL0,
integer, dimension(2*nk+2), intent(in)  KoR0,
real, dimension(2*nk+2), intent(in)  pL0,
real, dimension(2*nk+2), intent(in)  pR0,
real, dimension(2*nk+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]nkNumber of layers
[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 1442 of file MOM_neutral_diffusion.F90.

References compare_nsp_row().

Referenced by neutral_diffusion_unit_tests().

1442  logical, intent(in) :: verbose !< If true, write results to stdout
1443  integer, intent(in) :: nk !< Number of layers
1444  integer, dimension(2*nk+2), intent(in) :: kol !< Index of first left interface above neutral surface
1445  integer, dimension(2*nk+2), intent(in) :: kor !< Index of first right interface above neutral surface
1446  real, dimension(2*nk+2), intent(in) :: pl !< Fractional position of neutral surface within layer KoL of left column
1447  real, dimension(2*nk+2), intent(in) :: pr !< Fractional position of neutral surface within layer KoR of right column
1448  real, dimension(2*nk+1), intent(in) :: heff !< Effective thickness between two neutral surfaces (Pa)
1449  integer, dimension(2*nk+2), intent(in) :: kol0 !< Correct value for KoL
1450  integer, dimension(2*nk+2), intent(in) :: kor0 !< Correct value for KoR
1451  real, dimension(2*nk+2), intent(in) :: pl0 !< Correct value for pL
1452  real, dimension(2*nk+2), intent(in) :: pr0 !< Correct value for pR
1453  real, dimension(2*nk+1), intent(in) :: heff0 !< Correct value for hEff
1454  character(len=*), intent(in) :: title !< Title for messages
1455 
1456  ! Local variables
1457  integer :: k, stdunit
1458  logical :: this_row_failed
1459 
1460  test_nsp = .false.
1461  do k = 1,2*nk+2
1462  test_nsp = test_nsp .or. compare_nsp_row(kol(k), kor(k), pl(k), pr(k), kol0(k), kor0(k), pl0(k), pr0(k))
1463  if (k < 2*nk+2) then
1464  if (heff(k) /= heff0(k)) test_nsp = .true.
1465  endif
1466  enddo
1467 
1468  if (test_nsp .or. verbose) then
1469  stdunit = 6
1470  if (test_nsp) stdunit = 0 ! In case of wrong results, write to error stream
1471  write(stdunit,'(a)') title
1472  do k = 1,2*nk+2
1473  this_row_failed = compare_nsp_row(kol(k), kor(k), pl(k), pr(k), kol0(k), kor0(k), pl0(k), pr0(k))
1474  if (this_row_failed) then
1475  write(stdunit,10) k,kol(k),pl(k),kor(k),pr(k),' <-- WRONG!'
1476  write(stdunit,10) k,kol0(k),pl0(k),kor0(k),pr0(k),' <-- should be this'
1477  else
1478  write(stdunit,10) k,kol(k),pl(k),kor(k),pr(k)
1479  endif
1480  if (k < 2*nk+2) then
1481  if (heff(k) /= heff0(k)) then
1482  write(stdunit,'(i3,8x,"layer hEff =",2(f20.16,a))') k,heff(k)," .neq. ",heff0(k),' <-- WRONG!'
1483  else
1484  write(stdunit,'(i3,8x,"layer hEff =",f20.16)') k,heff(k)
1485  endif
1486  endif
1487  enddo
1488  endif
1489 
1490 10 format("ks=",i3," kL=",i3," pL=",f20.16," kR=",i3," pR=",f20.16,a)
Here is the call graph for this function:
Here is the caller graph for this function:

Variable Documentation

◆ debug_this_module

logical, parameter mom_neutral_diffusion::debug_this_module = .false.
private

Definition at line 59 of file MOM_neutral_diffusion.F90.

59 logical, parameter :: debug_this_module = .false.

◆ mdl

character(len=40) mom_neutral_diffusion::mdl = "MOM_neutral_diffusion"

Definition at line 57 of file MOM_neutral_diffusion.F90.

Referenced by neutral_diffusion_init().

57 character(len=40) :: mdl = "MOM_neutral_diffusion" ! module name