MOM6
MOM_neutral_diffusion.F90
Go to the documentation of this file.
1 !> A column-wise toolbox for implementing neutral diffusion
2 module mom_neutral_diffusion
3 
4 ! This file is part of MOM6. See LICENSE.md for the license.
5 
6 use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
7 use mom_cpu_clock, only : clock_module, clock_routine
8 use mom_diag_mediator, only : diag_ctrl, time_type
9 use mom_diag_mediator, only : post_data, register_diag_field
10 use mom_eos, only : eos_type, calculate_compress, calculate_density_derivs
11 use mom_error_handler, only : mom_error, fatal, warning, mom_mesg, is_root_pe
12 use mom_file_parser, only : get_param, log_version, param_file_type
13 use mom_file_parser, only : openparameterblock, closeparameterblock
14 use mom_grid, only : ocean_grid_type
15 use mom_tracer_registry,only : tracer_registry_type
16 use mom_verticalgrid, only : verticalgrid_type
17 
18 implicit none ; private
19 
20 #include <MOM_memory.h>
21 
22 public neutral_diffusion
23 public neutral_diffusion_init
24 public neutral_diffusion_diag_init
25 public neutral_diffusion_end
26 public neutral_diffusion_calc_coeffs
27 public neutral_diffusion_unit_tests
28 
29 type, public :: neutral_diffusion_cs ; private
30  integer :: nkp1 ! Number of interfaces for a column = nk + 1
31  integer :: nkp1x2 ! Number of intersecting interfaces between columns = 2 * nkp1
32 
33  real, allocatable, dimension(:,:,:) :: upol ! Non-dimensional position with left layer uKoL-1, u-point
34  real, allocatable, dimension(:,:,:) :: upor ! Non-dimensional position with right layer uKoR-1, u-point
35  integer, allocatable, dimension(:,:,:) :: ukol ! Index of left interface corresponding to neutral surface, u-point
36  integer, allocatable, dimension(:,:,:) :: ukor ! Index of right interface corresponding to neutral surface, u-point
37  real, allocatable, dimension(:,:,:) :: uheff ! Effective thickness at u-point (H units)
38  real, allocatable, dimension(:,:,:) :: vpol ! Non-dimensional position with left layer uKoL-1, v-point
39  real, allocatable, dimension(:,:,:) :: vpor ! Non-dimensional position with right layer uKoR-1, v-point
40  integer, allocatable, dimension(:,:,:) :: vkol ! Index of left interface corresponding to neutral surface, v-point
41  integer, allocatable, dimension(:,:,:) :: vkor ! Index of right interface corresponding to neutral surface, v-point
42  real, allocatable, dimension(:,:,:) :: vheff ! Effective thickness at v-point (H units)
43 
44  type(diag_ctrl), pointer :: diag ! structure to regulate output
45  integer, allocatable, dimension(:) :: id_neutral_diff_tracer_conc_tend ! tracer concentration tendency
46  integer, allocatable, dimension(:) :: id_neutral_diff_tracer_cont_tend ! tracer content tendency
47  integer, allocatable, dimension(:) :: id_neutral_diff_tracer_cont_tend_2d ! k-summed tracer content tendency
48  integer, allocatable, dimension(:) :: id_neutral_diff_tracer_trans_x_2d ! k-summed ndiff zonal tracer transport
49  integer, allocatable, dimension(:) :: id_neutral_diff_tracer_trans_y_2d ! k-summed ndiff merid tracer transport
50 
51  real :: c_p ! heat capacity of seawater (J kg-1 K-1)
52 
53 end type neutral_diffusion_cs
54 
55 ! This include declares and sets the variable "version".
56 #include "version_variable.h"
57 character(len=40) :: mdl = "MOM_neutral_diffusion" ! module name
58 
59 logical, parameter :: debug_this_module = .false.
60 
61 contains
62 
63 !> Read parameters and allocate control structure for neutral_diffusion module.
64 logical function neutral_diffusion_init(Time, G, param_file, diag, CS)
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 
114 end function neutral_diffusion_init
115 
116 
117 !> Diagnostic handles for neutral diffusion tendencies.
118 subroutine neutral_diffusion_diag_init(Time, G, diag, C_p, Reg, CS)
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 
250 end subroutine neutral_diffusion_diag_init
251 
252 
253 !> Calculate remapping factors for u/v columns used to map adjoining columns to
254 !! a shared coordinate space.
255 subroutine neutral_diffusion_calc_coeffs(G, GV, h, T, S, EOS, CS)
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 
311 end subroutine neutral_diffusion_calc_coeffs
312 
313 
314 !> Update tracer concentration due to neutral diffusion; layer thickness unchanged by this update.
315 subroutine neutral_diffusion(G, GV, h, Coef_x, Coef_y, Tracer, m, dt, name, CS)
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 
472 end subroutine neutral_diffusion
473 
474 
475 !> Returns interface scalar, Si, for a column of layer values, S.
476 subroutine interface_scalar(nk, h, S, Si, i_method)
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 
509 end subroutine interface_scalar
510 
511 !> Returns the PPM quasi-fourth order edge value at k+1/2 following
512 !! equation 1.6 in Colella & Woodward, 1984: JCP 54, 174-201.
513 real function ppm_edge(hkm1, hk, hkp1, hkp2, Ak, Akp1, Pk, Pkp1)
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 
549 end function ppm_edge
550 
551 !> Returns the average of a PPM reconstruction between two
552 !! fractional positions.
553 real function ppm_ave(xL, xR, aL, aR, aMean)
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 
578 end function ppm_ave
579 
580 !> A true signum function that returns either -abs(a), when x<0; or abs(a) when x>0; or 0 when x=0.
581 real function signum(a,x)
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 
588 end function signum
589 
590 !> Returns PLM slopes for a column where the slopes are the difference in value across each cell.
591 !! The limiting follows equation 1.8 in Colella & Woodward, 1984: JCP 54, 174-201.
592 subroutine plm_diff(nk, h, S, c_method, b_method, diff)
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 
653 end subroutine plm_diff
654 
655 !> Returns the cell-centered second-order finite volume (unlimited PLM) slope
656 !! using three consecutive cell widths and average values. Slope is returned
657 !! as a difference across the central cell (i.e. units of scalar S).
658 !! Discretization follows equation 1.7 in Colella & Woodward, 1984: JCP 54, 174-201.
659 real function fv_diff(hkm1, hk, hkp1, Skm1, Sk, Skp1)
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 ) )
679 end function fv_diff
680 
681 
682 !> Returns the cell-centered second-order weighted least squares slope
683 !! using three consecutive cell widths and average values. Slope is returned
684 !! as a gradient (i.e. units of scalar S over width units).
685 real function fvlsq_slope(hkm1, hk, hkp1, Skm1, Sk, Skp1)
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
711 end function fvlsq_slope
712 
713 
714 !> Returns positions within left/right columns of combined interfaces
715 subroutine find_neutral_surface_positions(nk, Pl, Tl, Sl, dRdTl, dRdSl, Pr, Tr, Sr, dRdTr, dRdSr, PoL, PoR, KoL, KoR, hEff)
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 
888 end subroutine find_neutral_surface_positions
889 
890 !> Converts non-dimensional position within a layer to absolute position (for debugging)
891 real function absolute_position(n,Pint,Karr,NParr,k_surface)
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 
904 end function absolute_position
905 
906 !> Converts non-dimensional positions within layers to absolute positions (for debugging)
907 function absolute_positions(n,Pint,Karr,NParr)
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 
922 end function absolute_positions
923 
924 !> Returns the non-dimensional position between Pneg and Ppos where the
925 !! interpolated density difference equals zero.
926 !! The result is always bounded to be between 0 and 1.
927 real function interpolate_for_nondim_position(dRhoNeg, Pneg, dRhoPos, Ppos)
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'
953 end function interpolate_for_nondim_position
954 
955 
956 !> Returns a single column of neutral diffusion fluxes of a tracer.
957 subroutine neutral_surface_flux(nk, hl, hr, Tl, Tr, PiL, PiR, KoL, KoR, hEff, Flx)
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 
1045 end subroutine neutral_surface_flux
1046 
1047 !> Returns true if unit tests of neutral_diffusion functions fail. Otherwise returns false.
1048 logical function neutral_diffusion_unit_tests(verbose)
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 
1278 end function neutral_diffusion_unit_tests
1279 
1280 !> Returns true if a test of fv_diff() fails, and conditionally writes results to stream
1281 logical function test_fv_diff(verbose, hkm1, hk, hkp1, Skm1, Sk, Skp1, Ptrue, title)
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 
1310 end function test_fv_diff
1311 
1312 !> Returns true if a test of fvlsq_slope() fails, and conditionally writes results to stream
1313 logical function test_fvlsq_slope(verbose, hkm1, hk, hkp1, Skm1, Sk, Skp1, Ptrue, title)
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 
1342 end function test_fvlsq_slope
1343 
1344 !> Returns true if a test of interpolate_for_nondim_position() fails, and conditionally writes results to stream
1345 logical function test_ifndp(verbose, rhoNeg, Pneg, rhoPos, Ppos, Ptrue, title)
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 
1372 end function test_ifndp
1373 
1374 !> Returns true if comparison of Po and Ptrue fails, and conditionally writes results to stream
1375 logical function test_data1d(verbose, nk, Po, Ptrue, title)
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 
1405 end function test_data1d
1406 
1407 !> Returns true if comparison of Po and Ptrue fails, and conditionally writes results to stream
1408 logical function test_data1di(verbose, nk, Po, Ptrue, title)
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 
1438 end function test_data1di
1439 
1440 !> Returns true if output of find_neutral_surface_positions() does not match correct values, and conditionally writes results to stream
1441 logical function test_nsp(verbose, nk, KoL, KoR, pL, pR, hEff, KoL0, KoR0, pL0, pR0, hEff0, title)
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)
1491 end function test_nsp
1492 
1493 !> Compares a single row, k, of output from find_neutral_surface_positions()
1494 logical function compare_nsp_row(KoL, KoR, pL, pR, KoL0, KoR0, pL0, pR0)
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.
1509 end function compare_nsp_row
1510 
1511 !> Deallocates neutral_diffusion control structure
1512 subroutine neutral_diffusion_end(CS)
1513  type(neutral_diffusion_cs), pointer :: CS
1514 
1515  if (associated(cs)) deallocate(cs)
1516 
1517 end subroutine neutral_diffusion_end
1518 
1519 end module mom_neutral_diffusion
mom_neutral_diffusion::test_data1d
logical function test_data1d(verbose, nk, Po, Ptrue, title)
Returns true if comparison of Po and Ptrue fails, and conditionally writes results to stream...
Definition: MOM_neutral_diffusion.F90:1376
mom_neutral_diffusion::absolute_positions
real function, dimension(2 *n+2) absolute_positions(n, Pint, Karr, NParr)
Converts non-dimensional positions within layers to absolute positions (for debugging) ...
Definition: MOM_neutral_diffusion.F90:908
mom_neutral_diffusion::fv_diff
real function fv_diff(hkm1, hk, hkp1, Skm1, Sk, Skp1)
Returns the cell-centered second-order finite volume (unlimited PLM) slope using three consecutive ce...
Definition: MOM_neutral_diffusion.F90:660
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:19
mom_neutral_diffusion::test_fv_diff
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.
Definition: MOM_neutral_diffusion.F90:1282
mom_diag_mediator::diag_ctrl
The following data type a list of diagnostic fields an their variants, as well as variables that cont...
Definition: MOM_diag_mediator.F90:139
mom_neutral_diffusion::test_data1di
logical function test_data1di(verbose, nk, Po, Ptrue, title)
Returns true if comparison of Po and Ptrue fails, and conditionally writes results to stream...
Definition: MOM_neutral_diffusion.F90:1409
mom_neutral_diffusion::neutral_surface_flux
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.
Definition: MOM_neutral_diffusion.F90:958
mom_file_parser::param_file_type
Definition: MOM_file_parser.F90:91
mom_grid
Provides the ocean grid type.
Definition: MOM_grid.F90:2
mom_neutral_diffusion::fvlsq_slope
real function fvlsq_slope(hkm1, hk, hkp1, Skm1, Sk, Skp1)
Returns the cell-centered second-order weighted least squares slope using three consecutive cell widt...
Definition: MOM_neutral_diffusion.F90:686
mom_neutral_diffusion::neutral_diffusion_unit_tests
logical function, public neutral_diffusion_unit_tests(verbose)
Returns true if unit tests of neutral_diffusion functions fail. Otherwise returns false...
Definition: MOM_neutral_diffusion.F90:1049
mom_neutral_diffusion::ppm_ave
real function ppm_ave(xL, xR, aL, aR, aMean)
Returns the average of a PPM reconstruction between two fractional positions.
Definition: MOM_neutral_diffusion.F90:554
mom_cpu_clock
Definition: MOM_cpu_clock.F90:1
mom_neutral_diffusion::neutral_diffusion_init
logical function, public neutral_diffusion_init(Time, G, param_file, diag, CS)
Read parameters and allocate control structure for neutral_diffusion module.
Definition: MOM_neutral_diffusion.F90:65
mom_neutral_diffusion::mdl
character(len=40) mdl
Definition: MOM_neutral_diffusion.F90:57
mom_file_parser
Definition: MOM_file_parser.F90:1
mom_neutral_diffusion::interface_scalar
subroutine interface_scalar(nk, h, S, Si, i_method)
Returns interface scalar, Si, for a column of layer values, S.
Definition: MOM_neutral_diffusion.F90:477
mom_neutral_diffusion::ppm_edge
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...
Definition: MOM_neutral_diffusion.F90:514
mom_neutral_diffusion::test_fvlsq_slope
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...
Definition: MOM_neutral_diffusion.F90:1314
mom_neutral_diffusion::test_ifndp
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 t...
Definition: MOM_neutral_diffusion.F90:1346
mom_neutral_diffusion::neutral_diffusion
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...
Definition: MOM_neutral_diffusion.F90:316
mom_eos::calculate_compress
subroutine, public calculate_compress(T, S, pressure, rho, drho_dp, start, npts, EOS)
Calls the appropriate subroutine to calculate the density and compressibility for 1-D array inputs...
Definition: MOM_EOS.F90:293
mom_diag_mediator::post_data
Definition: MOM_diag_mediator.F90:80
mom_neutral_diffusion::absolute_position
real function absolute_position(n, Pint, Karr, NParr, k_surface)
Converts non-dimensional position within a layer to absolute position (for debugging) ...
Definition: MOM_neutral_diffusion.F90:892
mom_tracer_registry
This module contains the tracer_registry_type and the subroutines that handle registration of tracers...
Definition: MOM_tracer_registry.F90:5
mom_neutral_diffusion::plm_diff
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...
Definition: MOM_neutral_diffusion.F90:593
mom_eos::calculate_density_derivs
subroutine, public calculate_density_derivs(T, S, pressure, drho_dT, drho_dS, start, npts, EOS)
Calls the appropriate subroutine to calculate density derivatives for 1-D array inputs.
Definition: MOM_EOS.F90:214
mom_file_parser::closeparameterblock
subroutine, public closeparameterblock(CS)
Definition: MOM_file_parser.F90:1826
mom_neutral_diffusion::neutral_diffusion_calc_coeffs
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 spac...
Definition: MOM_neutral_diffusion.F90:256
mom_file_parser::openparameterblock
subroutine, public openparameterblock(CS, blockName, desc)
Definition: MOM_file_parser.F90:1810
mom_tracer_registry::tracer_registry_type
Type to carry basic tracer information.
Definition: MOM_tracer_registry.F90:65
mom_error_handler::is_root_pe
logical function, public is_root_pe()
Definition: MOM_error_handler.F90:49
mom_diag_mediator
Definition: MOM_diag_mediator.F90:1
mom_error_handler
Definition: MOM_error_handler.F90:1
mom_neutral_diffusion::find_neutral_surface_positions
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.
Definition: MOM_neutral_diffusion.F90:716
mom_error_handler::mom_mesg
subroutine, public mom_mesg(message, verb, all_print)
Definition: MOM_error_handler.F90:57
mom_eos
Provides subroutines for quantities specific to the equation of state.
Definition: MOM_EOS.F90:2
mom_file_parser::log_version
Definition: MOM_file_parser.F90:140
mom_verticalgrid::verticalgrid_type
Definition: MOM_verticalGrid.F90:34
mom_neutral_diffusion::interpolate_for_nondim_position
real function interpolate_for_nondim_position(dRhoNeg, Pneg, dRhoPos, Ppos)
Returns the non-dimensional position between Pneg and Ppos where the interpolated density difference ...
Definition: MOM_neutral_diffusion.F90:928
mom_neutral_diffusion::neutral_diffusion_end
subroutine, public neutral_diffusion_end(CS)
Deallocates neutral_diffusion control structure.
Definition: MOM_neutral_diffusion.F90:1513
mom_neutral_diffusion::debug_this_module
logical, parameter debug_this_module
Definition: MOM_neutral_diffusion.F90:59
mom_neutral_diffusion
A column-wise toolbox for implementing neutral diffusion.
Definition: MOM_neutral_diffusion.F90:2
mom_neutral_diffusion::signum
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...
Definition: MOM_neutral_diffusion.F90:582
mom_neutral_diffusion::test_nsp
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...
Definition: MOM_neutral_diffusion.F90:1442
mom_neutral_diffusion::neutral_diffusion_cs
Definition: MOM_neutral_diffusion.F90:29
mom_verticalgrid
Definition: MOM_verticalGrid.F90:1
mom_diag_mediator::register_diag_field
integer function, public register_diag_field(module_name, field_name, axes, init_time, long_name, units, missing_value, range, mask_variant, standard_name, verbose, do_not_log, err_msg, interp_method, tile_count, cmor_field_name, cmor_long_name, cmor_units, cmor_standard_name, cell_methods, x_cell_method, y_cell_method, v_cell_method, conversion, v_extensive)
Returns the "diag_mediator" handle for a group (native, CMOR, z-coord, ...) of diagnostics derived fr...
Definition: MOM_diag_mediator.F90:1048
mom_error_handler::mom_error
subroutine, public mom_error(level, message, all_print)
Definition: MOM_error_handler.F90:73
mom_neutral_diffusion::neutral_diffusion_diag_init
subroutine, public neutral_diffusion_diag_init(Time, G, diag, C_p, Reg, CS)
Diagnostic handles for neutral diffusion tendencies.
Definition: MOM_neutral_diffusion.F90:119
mom_neutral_diffusion::compare_nsp_row
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()
Definition: MOM_neutral_diffusion.F90:1495
mom_file_parser::get_param
Definition: MOM_file_parser.F90:135
mom_eos::eos_type
A control structure for the equation of state.
Definition: MOM_EOS.F90:55