MOM6
Data Types | Functions/Subroutines
mom_neutral_diffusion_aux Module Reference

Detailed Description

A column-wise toolbox for implementing neutral diffusion.

Data Types

type  ndiff_aux_cs_type
 The control structure for this module. More...
 

Functions/Subroutines

subroutine, public set_ndiff_aux_params (CS, deg, max_iter, drho_tol, xtol, ref_pres, force_brent, EOS, debug)
 Initialize the parameters used to iteratively find the neutral direction. More...
 
subroutine, public mark_unstable_cells (nk, dRdT, dRdS, T, S, stable_cell, ns)
 Given the reconsturcitons of dRdT, dRdS, T, S mark the cells which are stably stratified parts of the water column For an layer to be unstable the top interface must be denser than the bottom or the bottom interface of the layer. More...
 
subroutine, public increment_interface (nk, kl, ki, stable, reached_bottom, searching_this_column, searching_other_column)
 Increments the interface which was just connected and also set flags if the bottom is reached. More...
 
real function, public calc_drho (T1, S1, dRdT1, dRdS1, T2, S2, dRdT2, dRdS2)
 Calculates difference in density at two points (rho1-rho2) with known density derivatives, T, and S. More...
 
subroutine, public drho_at_pos (CS, T_ref, S_ref, alpha_ref, beta_ref, P_top, P_bot, ppoly_T, ppoly_S, x0, delta_rho, P_out, T_out, S_out, alpha_avg_out, beta_avg_out, delta_T_out, delta_S_out)
 Calculate the difference in neutral density between a reference T, S, alpha, and beta at a point on the polynomial reconstructions of T, S. More...
 
subroutine, public search_other_column (dRhoTop, dRhoBot, Ptop, Pbot, lastP, lastK, kl, kl_0, ki, top_connected, bot_connected, out_P, out_K, search_layer)
 Searches the "other" (searched) column for the position of the neutral surface. More...
 
real function, public 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...
 
real function, public refine_nondim_position (CS, T_ref, S_ref, alpha_ref, beta_ref, P_top, P_bot, ppoly_T, ppoly_S, drho_top, drho_bot, min_bound)
 Use root-finding methods to find where dRho = 0, based on the equation of state and the polynomial reconstructions of temperature, salinity. Initial guess is based on the zero crossing of based on linear profiles of dRho, T, and S, between the top and bottom interface. If second derivatives of the EOS are available, it starts with a Newton's method. However, Newton's method is not guaranteed to be bracketed, a check is performed to see if it it diverges outside the interval. In that case (or in the case that second derivatives are not available), Brent's method is used following the implementation found at https://people.sc.fsu.edu/~jburkardt/f_src/brent/brent.f90. More...
 
subroutine, public check_neutral_positions (CS, Ptop_l, Pbot_l, Ptop_r, Pbot_r, PoL, PoR, Tl_coeffs, Tr_coeffs, Sl_coeffs, Sr_coeffs)
 
subroutine, public kahan_sum (sum, summand, c)
 Do a compensated sum to account for roundoff level. More...
 

Function/Subroutine Documentation

◆ calc_drho()

real function, public mom_neutral_diffusion_aux::calc_drho ( real, intent(in)  T1,
real, intent(in)  S1,
real, intent(in)  dRdT1,
real, intent(in)  dRdS1,
real, intent(in)  T2,
real, intent(in)  S2,
real, intent(in)  dRdT2,
real, intent(in)  dRdS2 
)

Calculates difference in density at two points (rho1-rho2) with known density derivatives, T, and S.

Parameters
[in]t1Temperature at point 1
[in]s1Salinity at point 1
[in]drdt1dRhodT at point 1
[in]drds1dRhodS at point 1
[in]t2Temperature at point 2
[in]s2Salinity at point 2
[in]drdt2dRhodT at point 2
[in]drds2dRhodS at point

Definition at line 151 of file MOM_neutral_diffusion_aux.F90.

151  real, intent(in ) :: T1 !< Temperature at point 1
152  real, intent(in ) :: S1 !< Salinity at point 1
153  real, intent(in ) :: dRdT1 !< dRhodT at point 1
154  real, intent(in ) :: dRdS1 !< dRhodS at point 1
155  real, intent(in ) :: T2 !< Temperature at point 2
156  real, intent(in ) :: S2 !< Salinity at point 2
157  real, intent(in ) :: dRdT2 !< dRhodT at point 2
158  real, intent(in ) :: dRdS2 !< dRhodS at point
159 
160  calc_drho = 0.5*( (drdt1+drdt2)*(t1-t2) + (drds1+drds2)*(s1-s2) )

◆ check_neutral_positions()

subroutine, public mom_neutral_diffusion_aux::check_neutral_positions ( type(ndiff_aux_cs_type), intent(in)  CS,
real, intent(in)  Ptop_l,
real, intent(in)  Pbot_l,
real, intent(in)  Ptop_r,
real, intent(in)  Pbot_r,
real, intent(in)  PoL,
real, intent(in)  PoR,
real, dimension(cs%nterm), intent(in)  Tl_coeffs,
real, dimension(cs%nterm), intent(in)  Tr_coeffs,
real, dimension(cs%nterm), intent(in)  Sl_coeffs,
real, dimension(cs%nterm), intent(in)  Sr_coeffs 
)
Parameters
[in]csControl structure with parameters for this module
[in]ptop_lPressure at top interface of left cell
[in]pbot_lPressure at bottom interface of left cell
[in]ptop_rPressure at top interface of right cell
[in]pbot_rPressure at bottom interface of right cell
[in]polNondim position in left cell
[in]porNondim position in right cell
[in]tl_coeffsT polynomial coefficients of left cell
[in]tr_coeffsT polynomial coefficients of right cell
[in]sl_coeffsS polynomial coefficients of left cell
[in]sr_coeffsS polynomial coefficients of right cell

Definition at line 620 of file MOM_neutral_diffusion_aux.F90.

620  type(ndiff_aux_CS_type), intent(in) :: CS !< Control structure with parameters for this module
621  real, intent(in) :: Ptop_l !< Pressure at top interface of left cell
622  real, intent(in) :: Pbot_l !< Pressure at bottom interface of left cell
623  real, intent(in) :: Ptop_r !< Pressure at top interface of right cell
624  real, intent(in) :: Pbot_r !< Pressure at bottom interface of right cell
625  real, intent(in) :: PoL !< Nondim position in left cell
626  real, intent(in) :: PoR !< Nondim position in right cell
627  real, dimension(CS%nterm), intent(in) :: Tl_coeffs !< T polynomial coefficients of left cell
628  real, dimension(CS%nterm), intent(in) :: Tr_coeffs !< T polynomial coefficients of right cell
629  real, dimension(CS%nterm), intent(in) :: Sl_coeffs !< S polynomial coefficients of left cell
630  real, dimension(CS%nterm), intent(in) :: Sr_coeffs !< S polynomial coefficients of right cell
631  ! Local variables
632  real :: Pl, Pr, Tl, Tr, Sl, Sr
633  real :: alpha_l, alpha_r, beta_l, beta_r
634  real :: delta_rho
635 
636  pl = (1. - pol)*ptop_l + pol*pbot_l
637  pr = (1. - por)*ptop_r + por*pbot_r
638  tl = evaluation_polynomial( tl_coeffs, cs%nterm, pol )
639  sl = evaluation_polynomial( sl_coeffs, cs%nterm, pol )
640  tr = evaluation_polynomial( tr_coeffs, cs%nterm, por )
641  sr = evaluation_polynomial( sr_coeffs, cs%nterm, por )
642 
643  if (cs%ref_pres<0.) then
644  call calculate_density_derivs( tl, sl, pl, alpha_l, beta_l, cs%EOS )
645  call calculate_density_derivs( tr, sr, pr, alpha_r, beta_r, cs%EOS )
646  else
647  call calculate_density_derivs( tl, sl, cs%ref_pres, alpha_l, beta_l, cs%EOS )
648  call calculate_density_derivs( tr, sr, cs%ref_pres, alpha_r, beta_r, cs%EOS )
649  endif
650 
651  delta_rho = 0.5*( (alpha_l + alpha_r)*(tl - tr) + (beta_l + beta_r)*(sl - sr) )
652  write(*,'(A,ES23.15)') "Delta-rho: ", delta_rho
653 

◆ drho_at_pos()

subroutine, public mom_neutral_diffusion_aux::drho_at_pos ( type(ndiff_aux_cs_type), intent(in)  CS,
real, intent(in)  T_ref,
real, intent(in)  S_ref,
real, intent(in)  alpha_ref,
real, intent(in)  beta_ref,
real, intent(in)  P_top,
real, intent(in)  P_bot,
real, dimension(cs%nterm), intent(in)  ppoly_T,
real, dimension(cs%nterm), intent(in)  ppoly_S,
real, intent(in)  x0,
real, intent(out)  delta_rho,
real, intent(out), optional  P_out,
real, intent(out), optional  T_out,
real, intent(out), optional  S_out,
real, intent(out), optional  alpha_avg_out,
real, intent(out), optional  beta_avg_out,
real, intent(out), optional  delta_T_out,
real, intent(out), optional  delta_S_out 
)

Calculate the difference in neutral density between a reference T, S, alpha, and beta at a point on the polynomial reconstructions of T, S.

Parameters
[in]csControl structure with parameters for this module
[in]t_refTemperature at reference surface
[in]s_refSalinity at reference surface
[in]alpha_refdRho/dT at reference surface
[in]beta_refdRho/dS at reference surface
[in]p_topPressure (Pa) at top interface of layer to be searched
[in]p_botPressure (Pa) at bottom interface
[in]ppoly_tCoefficients of T reconstruction
[in]ppoly_sCoefficients of S reconstruciton
[in]x0Nondimensional position to evaluate
[out]delta_rhoThe density difference from a reference value
[out]p_outPressure at point x0
[out]t_outTemperature at point x0
[out]s_outSalinity at point x0
[out]alpha_avg_outAverage of alpha between reference and x0
[out]beta_avg_outAverage of beta between reference and x0
[out]delta_t_outDifference in temperature between reference and x0
[out]delta_s_outDifference in salinity between reference and x0

Definition at line 167 of file MOM_neutral_diffusion_aux.F90.

167  type(ndiff_aux_CS_type), intent(in) :: CS !< Control structure with parameters for this module
168  real, intent(in) :: T_ref !< Temperature at reference surface
169  real, intent(in) :: S_ref !< Salinity at reference surface
170  real, intent(in) :: alpha_ref !< dRho/dT at reference surface
171  real, intent(in) :: beta_ref !< dRho/dS at reference surface
172  real, intent(in) :: P_top !< Pressure (Pa) at top interface of layer to be searched
173  real, intent(in) :: P_bot !< Pressure (Pa) at bottom interface
174  real, dimension(CS%nterm), intent(in) :: ppoly_T !< Coefficients of T reconstruction
175  real, dimension(CS%nterm), intent(in) :: ppoly_S !< Coefficients of S reconstruciton
176  real, intent(in) :: x0 !< Nondimensional position to evaluate
177  real, intent(out) :: delta_rho !< The density difference from a reference value
178  real, optional, intent(out) :: P_out !< Pressure at point x0
179  real, optional, intent(out) :: T_out !< Temperature at point x0
180  real, optional, intent(out) :: S_out !< Salinity at point x0
181  real, optional, intent(out) :: alpha_avg_out !< Average of alpha between reference and x0
182  real, optional, intent(out) :: beta_avg_out !< Average of beta between reference and x0
183  real, optional, intent(out) :: delta_T_out !< Difference in temperature between reference and x0
184  real, optional, intent(out) :: delta_S_out !< Difference in salinity between reference and x0
185 
186  real :: alpha, beta, alpha_avg, beta_avg, P_int, T, S, delta_T, delta_S
187 
188  p_int = (1. - x0)*p_top + x0*p_bot
189  t = evaluation_polynomial( ppoly_t, cs%nterm, x0 )
190  s = evaluation_polynomial( ppoly_s, cs%nterm, x0 )
191  ! Interpolated pressure if using locally referenced neutral density
192  if (cs%ref_pres<0.) then
193  call calculate_density_derivs( t, s, p_int, alpha, beta, cs%EOS )
194  else
195  ! Constant reference pressure (isopycnal)
196  call calculate_density_derivs( t, s, cs%ref_pres, alpha, beta, cs%EOS )
197  endif
198 
199  ! Calculate the f(P) term for Newton's method
200  alpha_avg = 0.5*( alpha + alpha_ref )
201  beta_avg = 0.5*( beta + beta_ref )
202  delta_t = t - t_ref
203  delta_s = s - s_ref
204  delta_rho = alpha_avg*delta_t + beta_avg*delta_s
205 
206  ! If doing a Newton step, these quantities are needed, otherwise they can just be optional
207  if (present(p_out)) p_out = p_int
208  if (present(t_out)) t_out = t
209  if (present(s_out)) s_out = s
210  if (present(alpha_avg_out)) alpha_avg_out = alpha_avg
211  if (present(beta_avg_out)) beta_avg_out = beta_avg
212  if (present(delta_t_out)) delta_t_out = delta_t
213  if (present(delta_s_out)) delta_s_out = delta_s
214 

◆ increment_interface()

subroutine, public mom_neutral_diffusion_aux::increment_interface ( integer, intent(in)  nk,
integer, intent(inout)  kl,
integer, intent(inout)  ki,
logical, dimension(nk), intent(in)  stable,
logical, intent(inout)  reached_bottom,
logical, intent(inout)  searching_this_column,
logical, intent(inout)  searching_other_column 
)

Increments the interface which was just connected and also set flags if the bottom is reached.

Parameters
[in]nkNumber of vertical levels
[in,out]klCurrent layer (potentially updated)
[in,out]kiCurrent interface
[in]stableTrue if the cell is stably stratified
[in,out]reached_bottomUpdated when kl == nk and ki == 2
[in,out]searching_this_columnUpdated when kl == nk and ki == 2
[in,out]searching_other_columnUpdated when kl == nk and ki == 2

Definition at line 116 of file MOM_neutral_diffusion_aux.F90.

116  integer, intent(in ) :: nk !< Number of vertical levels
117  integer, intent(inout) :: kl !< Current layer (potentially updated)
118  integer, intent(inout) :: ki !< Current interface
119  logical, dimension(nk), intent(in ) :: stable !< True if the cell is stably stratified
120  logical, intent(inout) :: reached_bottom !< Updated when kl == nk and ki == 2
121  logical, intent(inout) :: searching_this_column !< Updated when kl == nk and ki == 2
122  logical, intent(inout) :: searching_other_column !< Updated when kl == nk and ki == 2
123  integer :: k
124 
125  if (ki == 1) then
126  ki = 2
127  elseif ((ki == 2) .and. (kl < nk) ) then
128  do k = kl+1,nk
129  if (stable(kl)) then
130  kl = k
131  ki = 1
132  exit
133  endif
134  ! If we did not find another stable cell, then the current cell is essentially the bottom
135  ki = 2
136  reached_bottom = .true.
137  searching_this_column = .true.
138  searching_other_column = .false.
139  enddo
140  elseif ((kl == nk) .and. (ki==2)) then
141  reached_bottom = .true.
142  searching_this_column = .true.
143  searching_other_column = .false.
144  else
145  call mom_error(fatal,"Unanticipated eventuality in increment_interface")
146  endif

◆ interpolate_for_nondim_position()

real function, public mom_neutral_diffusion_aux::interpolate_for_nondim_position ( real, intent(in)  dRhoNeg,
real, intent(in)  Pneg,
real, intent(in)  dRhoPos,
real, intent(in)  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.

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

Definition at line 304 of file MOM_neutral_diffusion_aux.F90.

304  real, intent(in) :: dRhoNeg !< Negative density difference
305  real, intent(in) :: Pneg !< Position of negative density difference
306  real, intent(in) :: dRhoPos !< Positive density difference
307  real, intent(in) :: Ppos !< Position of positive density difference
308 
309  if (ppos<pneg) then
310  stop 'interpolate_for_nondim_position: Houston, we have a problem! Ppos<Pneg'
311  elseif (drhoneg>drhopos) then
312  write(0,*) 'dRhoNeg, Pneg, dRhoPos, Ppos=',drhoneg, pneg, drhopos, ppos
313  elseif (drhoneg>drhopos) then
314  stop 'interpolate_for_nondim_position: Houston, we have a problem! dRhoNeg>dRhoPos'
315  endif
316  if (ppos<=pneg) then ! Handle vanished or inverted layers
317  interpolate_for_nondim_position = 0.5
318  elseif ( drhopos - drhoneg > 0. ) then
319  interpolate_for_nondim_position = min( 1., max( 0., -drhoneg / ( drhopos - drhoneg ) ) )
320  elseif ( drhopos - drhoneg == 0) then
321  if (drhoneg>0.) then
322  interpolate_for_nondim_position = 0.
323  elseif (drhoneg<0.) then
324  interpolate_for_nondim_position = 1.
325  else ! dRhoPos = dRhoNeg = 0
326  interpolate_for_nondim_position = 0.5
327  endif
328  else ! dRhoPos - dRhoNeg < 0
329  interpolate_for_nondim_position = 0.5
330  endif
331  if ( interpolate_for_nondim_position < 0. ) &
332  stop 'interpolate_for_nondim_position: Houston, we have a problem! Pint < Pneg'
333  if ( interpolate_for_nondim_position > 1. ) &
334  stop 'interpolate_for_nondim_position: Houston, we have a problem! Pint > Ppos'

◆ kahan_sum()

subroutine, public mom_neutral_diffusion_aux::kahan_sum ( real, intent(inout)  sum,
real, intent(in)  summand,
real, intent(inout)  c 
)

Do a compensated sum to account for roundoff level.

Parameters
[in,out]sumRunning sum
[in]summandTerm to be added
[in,out]cKeep track of roundoff

Definition at line 658 of file MOM_neutral_diffusion_aux.F90.

658  real, intent(inout) :: sum !< Running sum
659  real, intent(in ) :: summand !< Term to be added
660  real ,intent(inout) :: c !< Keep track of roundoff
661  real :: y, t
662  y = summand - c
663  t = sum + y
664  c = (t-sum) - y
665  sum = t
666 

◆ mark_unstable_cells()

subroutine, public mom_neutral_diffusion_aux::mark_unstable_cells ( integer, intent(in)  nk,
real, dimension(nk,2), intent(in)  dRdT,
real, dimension(nk,2), intent(in)  dRdS,
real, dimension(nk,2), intent(in)  T,
real, dimension(nk,2), intent(in)  S,
logical, dimension(nk), intent(out)  stable_cell,
integer, intent(out)  ns 
)

Given the reconsturcitons of dRdT, dRdS, T, S mark the cells which are stably stratified parts of the water column For an layer to be unstable the top interface must be denser than the bottom or the bottom interface of the layer.

Parameters
[in]nkNumber of levels in a column
[in]drdtdrho/dT [kg m-3 degC-1] at interfaces
[in]drdsdrho/dS [kg m-3 ppt-1] at interfaces
[in]tTemperature [degC] at interfaces
[in]sSalinity [ppt] at interfaces
[out]stable_cellTrue if this cell is unstably stratified
[out]nsNumber of neutral surfaces in unmasked part of the column

Definition at line 64 of file MOM_neutral_diffusion_aux.F90.

64  integer, intent(in) :: nk !< Number of levels in a column
65  real, dimension(nk,2), intent(in) :: dRdT !< drho/dT [kg m-3 degC-1] at interfaces
66  real, dimension(nk,2), intent(in) :: dRdS !< drho/dS [kg m-3 ppt-1] at interfaces
67  real, dimension(nk,2), intent(in) :: T !< Temperature [degC] at interfaces
68  real, dimension(nk,2), intent(in) :: S !< Salinity [ppt] at interfaces
69  logical, dimension(nk), intent( out) :: stable_cell !< True if this cell is unstably stratified
70  integer, intent( out) :: ns !< Number of neutral surfaces in unmasked part of the column
71 
72  integer :: k, first_stable, prev_stable
73  real :: delta_rho
74 
75  ! First check to make sure that density profile between the two interfaces of the cell are stable
76  ! Note that we neglect a factor of 0.5 because we only care about the sign of delta_rho not magnitude
77  do k = 1,nk
78  ! Compare density of bottom interface to top interface, should be positive (or zero) if stable
79  delta_rho = (drdt(k,2) + drdt(k,1))*(t(k,2) - t(k,1)) + (drds(k,2) + drds(k,1))*(s(k,2) - s(k,1))
80  stable_cell(k) = delta_rho >= 0.
81  enddo
82 
83  first_stable = 1
84  ! Check to see that bottom interface of upper cell is lighter than the upper interface of the lower cell
85  do k=1,nk
86  if (stable_cell(k)) then
87  first_stable = k
88  exit
89  endif
90  enddo
91  prev_stable = first_stable
92 
93  ! Start either with the first stable cell or the layer just below the surface
94  do k = prev_stable+1, nk
95  ! Don't do anything if the cell has already been marked as unstable
96  if (.not. stable_cell(k)) cycle
97  ! Otherwise, we need to check to see if this cell's upper interface is denser than the previous stable_cell
98  ! Compare top interface of lower cell to bottom interface of upper cell, positive or zero if bottom cell is stable
99  delta_rho = (drdt(k,1) + drdt(prev_stable,2))*(t(k,1) - t(prev_stable,2)) + &
100  (drds(k,1) + drds(prev_stable,2))*(s(k,1) - s(prev_stable,2))
101  stable_cell(k) = delta_rho >= 0.
102  ! If the lower cell is marked as stable, then it should be the next reference cell
103  if (stable_cell(k)) prev_stable = k
104  enddo
105 
106  ! Number of interfaces is the 2 times number of stable cells in the water column
107  ns = 0
108  do k = 1,nk
109  if (stable_cell(k)) ns = ns + 2
110  enddo
111 

◆ refine_nondim_position()

real function, public mom_neutral_diffusion_aux::refine_nondim_position ( type(ndiff_aux_cs_type), intent(in)  CS,
real, intent(in)  T_ref,
real, intent(in)  S_ref,
real, intent(in)  alpha_ref,
real, intent(in)  beta_ref,
real, intent(in)  P_top,
real, intent(in)  P_bot,
real, dimension(:), intent(in)  ppoly_T,
real, dimension(:), intent(in)  ppoly_S,
real, intent(in)  drho_top,
real, intent(in)  drho_bot,
real, intent(in)  min_bound 
)

Use root-finding methods to find where dRho = 0, based on the equation of state and the polynomial reconstructions of temperature, salinity. Initial guess is based on the zero crossing of based on linear profiles of dRho, T, and S, between the top and bottom interface. If second derivatives of the EOS are available, it starts with a Newton's method. However, Newton's method is not guaranteed to be bracketed, a check is performed to see if it it diverges outside the interval. In that case (or in the case that second derivatives are not available), Brent's method is used following the implementation found at https://people.sc.fsu.edu/~jburkardt/f_src/brent/brent.f90.

Parameters
[in]csControl structure with parameters for this module
[in]t_refTemperature of the neutral surface at the searched from interface
[in]s_refSalinity of the neutral surface at the searched from interface
[in]alpha_refdRho/dT of the neutral surface at the searched from interface
[in]beta_refdRho/dS of the neutral surface at the searched from interface
[in]p_topPressure at the top interface in the layer to be searched
[in]p_botPressure at the bottom interface in the layer to be searched
[in]ppoly_tCoefficients of the order N polynomial reconstruction of T within the layer to be searched.
[in]ppoly_sCoefficients of the order N polynomial reconstruction of T within the layer to be searched.
[in]drho_topDelta rho at top interface (or previous position in cell
[in]drho_botDelta rho at bottom interface
[in]min_boundLower bound of position, also serves as the initial guess

Definition at line 346 of file MOM_neutral_diffusion_aux.F90.

346  type(ndiff_aux_CS_type), intent(in) :: CS !< Control structure with parameters for this module
347  real, intent(in) :: T_ref !< Temperature of the neutral surface at the searched from interface
348  real, intent(in) :: S_ref !< Salinity of the neutral surface at the searched from interface
349  real, intent(in) :: alpha_ref !< dRho/dT of the neutral surface at the searched from interface
350  real, intent(in) :: beta_ref !< dRho/dS of the neutral surface at the searched from interface
351  real, intent(in) :: P_top !< Pressure at the top interface in the layer to be searched
352  real, intent(in) :: P_bot !< Pressure at the bottom interface in the layer to be searched
353  real, dimension(:), intent(in) :: ppoly_T !< Coefficients of the order N polynomial reconstruction of T within
354  !! the layer to be searched.
355  real, dimension(:), intent(in) :: ppoly_S !< Coefficients of the order N polynomial reconstruction of T within
356  !! the layer to be searched.
357  real, intent(in) :: drho_top !< Delta rho at top interface (or previous position in cell
358  real, intent(in) :: drho_bot !< Delta rho at bottom interface
359  real, intent(in) :: min_bound !< Lower bound of position, also serves as the initial guess
360 
361  ! Local variables
362  integer :: form_of_EOS
363  integer :: iter
364  logical :: do_newton, do_brent
365 
366  real :: delta_rho, d_delta_rho_dP ! Terms for the Newton iteration
367  real :: P_int, P_min, P_ref ! Interpolated pressure
368  real :: delta_rho_init, delta_rho_final
369  real :: neg_x, neg_fun
370  real :: T, S, alpha, beta, alpha_avg, beta_avg
371  ! Newton's Method with variables
372  real :: dT_dP, dS_dP, delta_T, delta_S, delta_P
373  real :: dbeta_dS, dbeta_dT, dalpha_dT, dalpha_dS, dbeta_dP, dalpha_dP
374  real :: a, b, c, b_last
375  ! Extra Brent's Method variables
376  real :: d, e, f, fa, fb, fc, m, p, q, r, s0, sa, sb, tol, machep
377 
378  real :: P_last
379 
380  machep = epsilon(machep)
381  if (cs%ref_pres>=0.) p_ref = cs%ref_pres
382  delta_p = p_bot-p_top
383  refine_nondim_position = min_bound
384 
385  call extract_member_eos(cs%EOS, form_of_eos = form_of_eos)
386  do_newton = (form_of_eos == eos_linear) .or. (form_of_eos == eos_teos10) .or. (form_of_eos == eos_wright)
387  do_brent = .not. do_newton
388  if (cs%force_brent) then
389  do_newton = .not. cs%force_brent
390  do_brent = cs%force_brent
391  endif
392 
393  ! Calculate the initial values
394  call drho_at_pos(cs, t_ref, s_ref, alpha_ref, beta_ref, p_top, p_bot, ppoly_t, ppoly_s, min_bound, &
395  delta_rho, p_int, t, s, alpha_avg, beta_avg, delta_t, delta_s)
396  delta_rho_init = delta_rho
397  if ( abs(delta_rho_init) <= cs%drho_tol ) then
398  refine_nondim_position = min_bound
399  return
400  endif
401  if (abs(drho_bot) <= cs%drho_tol) then
402  refine_nondim_position = 1.
403  return
404  endif
405 
406  ! Set the initial values to ensure that the algorithm returns a 'negative' value
407  neg_fun = delta_rho
408  neg_x = min_bound
409 
410  if (cs%debug) then
411  write (*,*) "------"
412  write (*,*) "Starting x0, delta_rho: ", min_bound, delta_rho
413  endif
414 
415  ! For now only linear, Wright, and TEOS-10 equations of state have functions providing second derivatives and
416  ! thus can use Newton's method for the equation of state
417  if (do_newton) then
418  refine_nondim_position = min_bound
419  ! Set lower bound of pressure
420  p_min = p_top*(1.-min_bound) + p_bot*(min_bound)
421  fa = delta_rho_init ; a = min_bound
422  fb = delta_rho_init ; b = min_bound
423  fc = drho_bot ; c = 1.
424  ! Iterate over Newton's method for the function: x0 = x0 - delta_rho/d_delta_rho_dP
425  do iter = 1, cs%max_iter
426  p_int = p_top*(1. - b) + p_bot*b
427  ! Evaluate total derivative of delta_rho
428  if (cs%ref_pres<0.) p_ref = p_int
429  call calculate_density_second_derivs( t, s, p_ref, dbeta_ds, dbeta_dt, dalpha_dt, dbeta_dp, dalpha_dp, cs%EOS )
430  ! In the case of a constant reference pressure, no dependence on neutral direction with pressure
431  if (cs%ref_pres>=0.) then
432  dalpha_dp = 0. ; dbeta_dp = 0.
433  endif
434  dalpha_ds = dbeta_dt ! Cross derivatives are identicial
435  ! By chain rule dT_dP= (dT_dz)*(dz/dP) = dT_dz / (Pbot-Ptop)
436  dt_dp = first_derivative_polynomial( ppoly_t, cs%nterm, b ) / delta_p
437  ds_dp = first_derivative_polynomial( ppoly_s, cs%nterm, b ) / delta_p
438  ! Total derivative of d_delta_rho wrt P
439  d_delta_rho_dp = 0.5*( delta_s*(ds_dp*dbeta_ds + dt_dp*dbeta_dt + dbeta_dp) + &
440  ( delta_t*(ds_dp*dalpha_ds + dt_dp*dalpha_dt + dalpha_dp))) + &
441  ds_dp*beta_avg + dt_dp*alpha_avg
442  ! This probably won't happen, but if it does take a bisection
443  if (d_delta_rho_dp == 0.) then
444  b = 0.5*(a+c)
445  else
446  ! Newton step update
447  p_int = p_int - (fb / d_delta_rho_dp)
448  ! This line is equivalent to the next
449  ! refine_nondim_position = (P_top-P_int)/(P_top-P_bot)
450  b_last = b
451  b = (p_int-p_top)/delta_p
452  ! Test to see if it fell out of the bracketing interval. If so, take a bisection step
453  if (b < a .or. b > c) then
454  b = 0.5*(a + c)
455  endif
456  endif
457  call drho_at_pos(cs, t_ref, s_ref, alpha_ref, beta_ref, p_top, p_bot, ppoly_t, ppoly_s, &
458  b, fb, p_int, t, s, alpha_avg, beta_avg, delta_t, delta_s)
459  if (cs%debug) write(*,'(A,I3.3,X,ES23.15,X,ES23.15)') "Iteration, b, fb: ", iter, b, fb
460 
461  if (fb < 0. .and. fb > neg_fun) then
462  neg_fun = fb
463  neg_x = b
464  endif
465 
466  ! For the logic to find neutral surfaces to work properly, the function needs to converge to zero
467  ! or a small negative value
468  if ((fb <= 0.) .and. (fb >= -cs%drho_tol)) then
469  refine_nondim_position = b
470  exit
471  endif
472  ! Exit if method has stalled out
473  if ( abs(b-b_last)<=cs%xtol ) then
474  refine_nondim_position = b
475  exit
476  endif
477 
478  ! Update the bracket
479  if (sign(1.,fa)*sign(1.,fb)<0.) then
480  c = b
481  fc = delta_rho
482  else
483  a = b
484  fa = delta_rho
485  endif
486  enddo
487  refine_nondim_position = b
488  delta_rho = fb
489  endif
490  if (delta_rho > 0.) then
491  refine_nondim_position = neg_x
492  delta_rho = neg_fun
493  endif
494  ! Do Brent if analytic second derivatives don't exist
495  if (do_brent) then
496  sa = max(refine_nondim_position,min_bound) ; fa = delta_rho
497  sb = 1. ; fb = drho_bot
498  c = sa ; fc = fa ; e = sb - sa; d = e
499 
500 
501  ! This is from https://people.sc.fsu.edu/~jburkardt/f_src/brent/brent.f90
502  do iter = 1,cs%max_iter
503  if ( abs( fc ) < abs( fb ) ) then
504  sa = sb
505  sb = c
506  c = sa
507  fa = fb
508  fb = fc
509  fc = fa
510  endif
511  tol = 2. * machep * abs( sb ) + cs%xtol
512  m = 0.5 * ( c - sb )
513  if ( abs( m ) <= tol .or. fb == 0. ) then
514  exit
515  endif
516  if ( abs( e ) < tol .or. abs( fa ) <= abs( fb ) ) then
517  e = m
518  d = e
519  else
520  s0 = fb / fa
521  if ( sa == c ) then
522  p = 2. * m * s0
523  q = 1. - s0
524  else
525  q = fa / fc
526  r = fb / fc
527  p = s0 * ( 2. * m * q * ( q - r ) - ( sb - sa ) * ( r - 1. ) )
528  q = ( q - 1. ) * ( r - 1. ) * ( s0 - 1. )
529  endif
530  if ( 0. < p ) then
531  q = - q
532  else
533  p = - p
534  endif
535  s0 = e
536  e = d
537  if ( 2. * p < 3. * m * q - abs( tol * q ) .and. &
538  p < abs( 0.5 * s0 * q ) ) then
539  d = p / q
540  else
541  e = m
542  d = e
543  endif
544  endif
545  sa = sb
546  fa = fb
547  if ( tol < abs( d ) ) then
548  sb = sb + d
549  elseif ( 0. < m ) then
550  sb = sb + tol
551  else
552  sb = sb - tol
553  endif
554  call drho_at_pos(cs, t_ref, s_ref, alpha_ref, beta_ref, p_top, p_bot, ppoly_t, ppoly_s, &
555  sb, fb)
556  if ( ( 0. < fb .and. 0. < fc ) .or. &
557  ( fb <= 0. .and. fc <= 0. ) ) then
558  c = sa
559  fc = fa
560  e = sb - sa
561  d = e
562  endif
563  enddo
564  ! Modified from original to ensure that the minimum is found
565  fa = abs(fa) ; fb = abs(fb) ; fc = abs(fc)
566  delta_rho = min(fa, fb, fc)
567 
568  if (fb==delta_rho) then
569  refine_nondim_position = max(sb,min_bound)
570  elseif (fa==delta_rho) then
571  refine_nondim_position = max(sa,min_bound)
572  elseif (fc==delta_rho) then
573  refine_nondim_position = max(c, min_bound)
574  endif
575  endif
576 
577  ! Make sure that the result is bounded between 0 and 1
578  if (refine_nondim_position>1.) then
579  if (cs%debug) then
580  write (*,*) "T, T Poly Coeffs: ", t, ppoly_t
581  write (*,*) "S, S Poly Coeffs: ", s, ppoly_s
582  write (*,*) "T_ref, alpha_ref: ", t_ref, alpha_ref
583  write (*,*) "S_ref, beta_ref : ", s_ref, beta_ref
584  write (*,*) "x0: ", min_bound
585  write (*,*) "refine_nondim_position: ", refine_nondim_position
586  endif
587  call mom_error(warning, "refine_nondim_position>1.")
588  refine_nondim_position = 1.
589  endif
590 
591  if (refine_nondim_position<min_bound) then
592  if (cs%debug) then
593  write (*,*) "T, T Poly Coeffs: ", t, ppoly_t
594  write (*,*) "S, S Poly Coeffs: ", s, ppoly_s
595  write (*,*) "T_ref, alpha_ref: ", t_ref, alpha_ref
596  write (*,*) "S_ref, beta_ref : ", s_ref, beta_ref
597  write (*,*) "x0: ", min_bound
598  write (*,*) "refine_nondim_position: ", refine_nondim_position
599  endif
600  call mom_error(warning, "refine_nondim_position<min_bound.")
601  refine_nondim_position = min_bound
602  endif
603 
604  if (cs%debug) then
605  call drho_at_pos(cs, t_ref, s_ref, alpha_ref, beta_ref, p_top, p_bot, ppoly_t, ppoly_s, &
606  refine_nondim_position, delta_rho)
607  write (*,*) "End delta_rho: ", delta_rho
608  write (*,*) "x0, delta_x: ", min_bound, refine_nondim_position-min_bound
609  write (*,*) "refine_nondim_position: ", refine_nondim_position
610  write (*,*) "Iterations: ", iter
611  write (*,*) "T_poly_coeffs: ", ppoly_t
612  write (*,*) "S_poly_coeffs: ", ppoly_s
613  write (*,*) "******"
614  endif
615 

◆ search_other_column()

subroutine, public mom_neutral_diffusion_aux::search_other_column ( real, intent(in)  dRhoTop,
real, intent(in)  dRhoBot,
real, intent(in)  Ptop,
real, intent(in)  Pbot,
real, intent(in)  lastP,
integer, intent(in)  lastK,
integer, intent(in)  kl,
integer, intent(in)  kl_0,
integer, intent(in)  ki,
logical, dimension(:), intent(inout)  top_connected,
logical, dimension(:), intent(inout)  bot_connected,
real, intent(out)  out_P,
integer, intent(out)  out_K,
logical, intent(out)  search_layer 
)

Searches the "other" (searched) column for the position of the neutral surface.

Parameters
[in]drhotopDensity difference across top interface
[in]drhobotDensity difference across top interface
[in]ptopPressure at top interface
[in]pbotPressure at bottom interface
[in]lastpLast position connected in the searched column
[in]lastkLast layer connected in the searched column
[in]klLayer in the searched column
[in]kl_0Layer in the searched column
[in]kiInterface of the searched column
[in,out]top_connectedTrue if the top interface was pointed to
[in,out]bot_connectedTrue if the top interface was pointed to
[out]out_pPosition within searched column
[out]out_kLayer within searched column
[out]search_layerNeutral surface within cell

Definition at line 220 of file MOM_neutral_diffusion_aux.F90.

220  real, intent(in ) :: dRhoTop !< Density difference across top interface
221  real, intent(in ) :: dRhoBot !< Density difference across top interface
222  real, intent(in ) :: Ptop !< Pressure at top interface
223  real, intent(in ) :: Pbot !< Pressure at bottom interface
224  real, intent(in ) :: lastP !< Last position connected in the searched column
225  integer, intent(in ) :: lastK !< Last layer connected in the searched column
226  integer, intent(in ) :: kl !< Layer in the searched column
227  integer, intent(in ) :: kl_0 !< Layer in the searched column
228  integer, intent(in ) :: ki !< Interface of the searched column
229  logical, dimension(:), intent(inout) :: top_connected !< True if the top interface was pointed to
230  logical, dimension(:), intent(inout) :: bot_connected !< True if the top interface was pointed to
231  real, intent( out) :: out_P !< Position within searched column
232  integer, intent( out) :: out_K !< Layer within searched column
233  logical, intent( out) :: search_layer !< Neutral surface within cell
234 
235  search_layer = .false.
236  if (kl > kl_0) then ! Away from top cell
237  if (kl == lastk) then ! Searching in the same layer
238  if (drhotop > 0.) then
239  if (lastk == kl) then
240  out_p = lastp
241  else
242  out_p = 0.
243  endif
244  out_k = kl
245 ! out_P = max(0.,lastP) ; out_K = kl
246  elseif ( drhotop == drhobot ) then
247  if (top_connected(kl)) then
248  out_p = 1. ; out_k = kl
249  else
250  out_p = max(0.,lastp) ; out_k = kl
251  endif
252  elseif (drhotop >= drhobot) then
253  out_p = 1. ; out_k = kl
254  elseif (drhotop < 0. .and. drhobot < 0.)then
255  out_p = 1. ; out_k = kl
256  else
257  out_k = kl
258  out_p = max(interpolate_for_nondim_position( drhotop, ptop, drhobot, pbot ),lastp)
259  search_layer = .true.
260  endif
261  else ! Searching across the interface
262  if (.not. bot_connected(kl-1) ) then
263  out_k = kl-1
264  out_p = 1.
265  else
266  out_k = kl
267  out_p = 0.
268  endif
269  endif
270  else ! At the top cell
271  if (ki == 1) then
272  out_p = 0. ; out_k = kl
273  elseif (drhotop > 0.) then
274  if (lastk == kl) then
275  out_p = lastp
276  else
277  out_p = 0.
278  endif
279  out_k = kl
280 ! out_P = max(0.,lastP) ; out_K = kl
281  elseif ( drhotop == drhobot ) then
282  if (top_connected(kl)) then
283  out_p = 1. ; out_k = kl
284  else
285  out_p = max(0.,lastp) ; out_k = kl
286  endif
287  elseif (drhotop >= drhobot) then
288  out_p = 1. ; out_k = kl
289  elseif (drhotop < 0. .and. drhobot < 0.)then
290  out_p = 1. ; out_k = kl
291  else
292  out_k = kl
293  out_p = max(interpolate_for_nondim_position( drhotop, ptop, drhobot, pbot ),lastp)
294  search_layer = .true.
295  endif
296  endif
297 

◆ set_ndiff_aux_params()

subroutine, public mom_neutral_diffusion_aux::set_ndiff_aux_params ( type(ndiff_aux_cs_type), intent(inout)  CS,
integer, intent(in), optional  deg,
integer, intent(in), optional  max_iter,
real, intent(in), optional  drho_tol,
real, intent(in), optional  xtol,
real, intent(in), optional  ref_pres,
logical, intent(in), optional  force_brent,
type(eos_type), intent(in), optional, target  EOS,
logical, intent(in), optional  debug 
)

Initialize the parameters used to iteratively find the neutral direction.

Parameters
[in,out]csControl structure for refine_pos
[in]degDegree of polynommial used in reconstruction
[in]max_iterMaximum number of iterations
[in]drho_tolTolerance for function convergence
[in]xtolTolerance for change in position
[in]ref_presReference pressure to use
[in]force_brentForce Brent method for linear, TEOS-10, and WRIGHT
[in]debugIf true, print output use to help debug neutral diffusion
[in]eosEquation of state

Definition at line 40 of file MOM_neutral_diffusion_aux.F90.

40  type(ndiff_aux_CS_type), intent(inout) :: CS !< Control structure for refine_pos
41  integer, optional, intent(in ) :: deg !< Degree of polynommial used in reconstruction
42  integer, optional, intent(in ) :: max_iter !< Maximum number of iterations
43  real, optional, intent(in ) :: drho_tol !< Tolerance for function convergence
44  real, optional, intent(in ) :: xtol !< Tolerance for change in position
45  real, optional, intent(in ) :: ref_pres !< Reference pressure to use
46  logical, optional, intent(in ) :: force_brent !< Force Brent method for linear, TEOS-10, and WRIGHT
47  logical, optional, intent(in ) :: debug !< If true, print output use to help debug neutral diffusion
48  type(EOS_type), target, optional, intent(in ) :: EOS !< Equation of state
49 
50  if (present( deg )) cs%nterm = deg + 1
51  if (present( max_iter )) cs%max_iter = max_iter
52  if (present( drho_tol )) cs%drho_tol = drho_tol
53  if (present( xtol )) cs%xtol = xtol
54  if (present( ref_pres )) cs%ref_pres = ref_pres
55  if (present( force_brent )) cs%force_brent = force_brent
56  if (present( eos )) cs%EOS => eos
57  if (present( debug )) cs%debug = debug
58