6 use mom_cpu_clock,
only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
31 implicit none ;
private
33 #include <MOM_memory.h>
44 logical :: continuous_reconstruction = .true.
45 logical :: debug = .false.
46 logical :: hard_fail_heff
51 logical :: interior_only
54 real,
allocatable,
dimension(:,:,:) :: upol
55 real,
allocatable,
dimension(:,:,:) :: upor
56 integer,
allocatable,
dimension(:,:,:) :: ukol
58 integer,
allocatable,
dimension(:,:,:) :: ukor
60 real,
allocatable,
dimension(:,:,:) :: uheff
61 real,
allocatable,
dimension(:,:,:) :: vpol
62 real,
allocatable,
dimension(:,:,:) :: vpor
63 integer,
allocatable,
dimension(:,:,:) :: vkol
65 integer,
allocatable,
dimension(:,:,:) :: vkor
67 real,
allocatable,
dimension(:,:,:) :: vheff
69 real,
allocatable,
dimension(:,:,:,:) :: ppoly_coeffs_t
70 real,
allocatable,
dimension(:,:,:,:) :: ppoly_coeffs_s
72 real,
allocatable,
dimension(:,:,:) :: drdt
73 real,
allocatable,
dimension(:,:,:) :: drds
74 real,
allocatable,
dimension(:,:,:) :: tint
75 real,
allocatable,
dimension(:,:,:) :: sint
76 real,
allocatable,
dimension(:,:,:) :: pint
78 real,
allocatable,
dimension(:,:,:,:) :: t_i
79 real,
allocatable,
dimension(:,:,:,:) :: s_i
80 real,
allocatable,
dimension(:,:,:,:) :: p_i
81 real,
allocatable,
dimension(:,:,:,:) :: drdt_i
82 real,
allocatable,
dimension(:,:,:,:) :: drds_i
83 integer,
allocatable,
dimension(:,:) :: ns
84 logical,
allocatable,
dimension(:,:,:) :: stable_cell
87 integer :: neutral_pos_method
88 character(len=40) :: delta_rho_form
91 integer :: id_uheff_2d = -1
92 integer :: id_vheff_2d = -1
97 type(
kpp_cs),
pointer :: kpp_csp => null()
102 #include "version_variable.h"
103 character(len=40) ::
mdl =
"MOM_neutral_diffusion"
109 type(time_type),
target,
intent(in) :: time
111 type(
diag_ctrl),
target,
intent(inout) :: diag
113 type(
eos_type),
target,
intent(in) :: eos
118 character(len=256) :: mesg
119 character(len=80) :: string
120 logical :: boundary_extrap
122 if (
associated(cs))
then
123 call mom_error(fatal,
"neutral_diffusion_init called with associated control structure.")
130 "This module implements neutral diffusion of tracers")
132 "If true, enables the neutral diffusion module.", &
145 call get_param(param_file,
mdl,
"NDIFF_CONTINUOUS", cs%continuous_reconstruction, &
146 "If true, uses a continuous reconstruction of T and S when "//&
147 "finding neutral surfaces along which diffusion will happen. "//&
148 "If false, a PPM discontinuous reconstruction of T and S "//&
149 "is done which results in a higher order routine but exacts "//&
150 "a higher computational cost.", default=.true.)
151 call get_param(param_file,
mdl,
"NDIFF_REF_PRES", cs%ref_pres, &
152 "The reference pressure (Pa) used for the derivatives of "//&
153 "the equation of state. If negative (default), local "//&
154 "pressure is used.", &
156 call get_param(param_file,
mdl,
"NDIFF_INTERIOR_ONLY", cs%interior_only, &
157 "If true, only applies neutral diffusion in the ocean interior."//&
158 "That is, the algorithm will exclude the surface and bottom"//&
159 "boundary layers.",default = .false.)
162 if (cs%continuous_reconstruction .eqv. .false.)
then
163 call get_param(param_file,
mdl,
"NDIFF_BOUNDARY_EXTRAP", boundary_extrap, &
164 "Extrapolate at the top and bottommost cells, otherwise \n"// &
165 "assume boundaries are piecewise constant", &
167 call get_param(param_file,
mdl,
"NDIFF_REMAPPING_SCHEME", string, &
168 "This sets the reconstruction scheme used "//&
169 "for vertical remapping for all variables. "//&
170 "It can be one of the following schemes: "//&
172 call initialize_remapping( cs%remap_CS, string, boundary_extrapolation = boundary_extrap )
173 call extract_member_remapping_cs(cs%remap_CS, degree=cs%deg)
174 call get_param(param_file,
mdl,
"NEUTRAL_POS_METHOD", cs%neutral_pos_method, &
175 "Method used to find the neutral position \n"// &
176 "1. Delta_rho varies linearly, find 0 crossing \n"// &
177 "2. Alpha and beta vary linearly from top to bottom, \n"// &
178 " Newton's method for neutral position \n"// &
179 "3. Full nonlinear equation of state, use regula falsi \n"// &
180 " for neutral position", default=3)
181 if (cs%neutral_pos_method > 4 .or. cs%neutral_pos_method < 0)
then
182 call mom_error(fatal,
"Invalid option for NEUTRAL_POS_METHOD")
185 call get_param(param_file,
mdl,
"DELTA_RHO_FORM", cs%delta_rho_form, &
186 "Determine how the difference in density is calculated \n"// &
187 " full : Difference of in-situ densities \n"// &
188 " no_pressure: Calculated from dRdT, dRdS, but no \n"// &
189 " pressure dependence", &
190 default=
"mid_pressure")
191 if (cs%neutral_pos_method > 1)
then
192 call get_param(param_file,
mdl,
"NDIFF_DRHO_TOL", cs%drho_tol, &
193 "Sets the convergence criterion for finding the neutral\n"// &
194 "position within a layer in kg m-3.", &
196 call get_param(param_file,
mdl,
"NDIFF_X_TOL", cs%x_tol, &
197 "Sets the convergence criterion for a change in nondim\n"// &
198 "position within a layer.", &
200 call get_param(param_file,
mdl,
"NDIFF_MAX_ITER", cs%max_iter, &
201 "The maximum number of iterations to be done before \n"// &
202 "exiting the iterative loop to find the neutral surface", &
205 call get_param(param_file,
mdl,
"NDIFF_DEBUG", cs%debug, &
206 "Turns on verbose output for discontinuous neutral "//&
207 "diffusion routines.", &
209 call get_param(param_file,
mdl,
"HARD_FAIL_HEFF", cs%hard_fail_heff, &
210 "Bring down the model if a problem with heff is detected",&
214 if (cs%interior_only)
then
217 if ( .not.
ASSOCIATED(cs%energetic_PBL_CSp) .and. .not.
ASSOCIATED(cs%KPP_CSp) )
then
218 call mom_error(fatal,
"NDIFF_INTERIOR_ONLY is true, but no valid boundary layer scheme was found")
226 if (cs%continuous_reconstruction)
then
228 allocate(cs%dRdT(szi_(g),szj_(g),szk_(g)+1)) ; cs%dRdT(:,:,:) = 0.
229 allocate(cs%dRdS(szi_(g),szj_(g),szk_(g)+1)) ; cs%dRdS(:,:,:) = 0.
232 allocate(cs%T_i(szi_(g),szj_(g),szk_(g),2)) ; cs%T_i(:,:,:,:) = 0.
233 allocate(cs%S_i(szi_(g),szj_(g),szk_(g),2)) ; cs%S_i(:,:,:,:) = 0.
234 allocate(cs%P_i(szi_(g),szj_(g),szk_(g),2)) ; cs%P_i(:,:,:,:) = 0.
235 allocate(cs%dRdT_i(szi_(g),szj_(g),szk_(g),2)) ; cs%dRdT_i(:,:,:,:) = 0.
236 allocate(cs%dRdS_i(szi_(g),szj_(g),szk_(g),2)) ; cs%dRdS_i(:,:,:,:) = 0.
237 allocate(cs%ppoly_coeffs_T(szi_(g),szj_(g),szk_(g),cs%deg+1)) ; cs%ppoly_coeffs_T(:,:,:,:) = 0.
238 allocate(cs%ppoly_coeffs_S(szi_(g),szj_(g),szk_(g),cs%deg+1)) ; cs%ppoly_coeffs_S(:,:,:,:) = 0.
239 allocate(cs%ns(szi_(g),szj_(g))) ; cs%ns(:,:) = 0.
242 allocate(cs%Tint(szi_(g),szj_(g),szk_(g)+1)) ; cs%Tint(:,:,:) = 0.
243 allocate(cs%Sint(szi_(g),szj_(g),szk_(g)+1)) ; cs%Sint(:,:,:) = 0.
244 allocate(cs%Pint(szi_(g),szj_(g),szk_(g)+1)) ; cs%Pint(:,:,:) = 0.
245 allocate(cs%stable_cell(szi_(g),szj_(g),szk_(g))) ; cs%stable_cell(:,:,:) = .true.
247 allocate(cs%uPoL(g%isd:g%ied,g%jsd:g%jed, cs%nsurf)); cs%uPoL(g%isc-1:g%iec,g%jsc:g%jec,:) = 0.
248 allocate(cs%uPoR(g%isd:g%ied,g%jsd:g%jed, cs%nsurf)); cs%uPoR(g%isc-1:g%iec,g%jsc:g%jec,:) = 0.
249 allocate(cs%uKoL(g%isd:g%ied,g%jsd:g%jed, cs%nsurf)); cs%uKoL(g%isc-1:g%iec,g%jsc:g%jec,:) = 0
250 allocate(cs%uKoR(g%isd:g%ied,g%jsd:g%jed, cs%nsurf)); cs%uKoR(g%isc-1:g%iec,g%jsc:g%jec,:) = 0
251 allocate(cs%uHeff(g%isd:g%ied,g%jsd:g%jed,cs%nsurf-1)); cs%uHeff(g%isc-1:g%iec,g%jsc:g%jec,:) = 0
253 allocate(cs%vPoL(g%isd:g%ied,g%jsd:g%jed, cs%nsurf)); cs%vPoL(g%isc:g%iec,g%jsc-1:g%jec,:) = 0.
254 allocate(cs%vPoR(g%isd:g%ied,g%jsd:g%jed, cs%nsurf)); cs%vPoR(g%isc:g%iec,g%jsc-1:g%jec,:) = 0.
255 allocate(cs%vKoL(g%isd:g%ied,g%jsd:g%jed, cs%nsurf)); cs%vKoL(g%isc:g%iec,g%jsc-1:g%jec,:) = 0
256 allocate(cs%vKoR(g%isd:g%ied,g%jsd:g%jed, cs%nsurf)); cs%vKoR(g%isc:g%iec,g%jsc-1:g%jec,:) = 0
257 allocate(cs%vHeff(g%isd:g%ied,g%jsd:g%jed,cs%nsurf-1)); cs%vHeff(g%isc:g%iec,g%jsc-1:g%jec,:) = 0
267 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
268 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: t
269 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: s
275 real,
dimension(SZK_(G),2) :: ppoly_r_s
276 real,
dimension(SZI_(G), SZJ_(G)) :: heff_sum
277 real,
dimension(SZI_(G),SZJ_(G)) :: hbl
279 real,
dimension(SZI_(G)) :: ref_pres
280 real,
dimension(SZI_(G)) :: rho_tmp
281 real :: h_neglect, h_neglect_edge
282 integer,
dimension(SZI_(G), SZJ_(G)) :: k_top
283 real,
dimension(SZI_(G), SZJ_(G)) :: zeta_top
285 integer,
dimension(SZI_(G), SZJ_(G)) :: k_bot
286 real,
dimension(SZI_(G), SZJ_(G)) :: zeta_bot
289 pa_to_h = 1. / gv%H_to_pa
291 k_top(:,:) = 1 ; k_bot(:,:) = 1
292 zeta_top(:,:) = 0. ; zeta_bot(:,:) = 1.
295 if (cs%interior_only)
then
297 if (
ASSOCIATED(cs%KPP_CSp))
call kpp_get_bld(cs%KPP_CSp, hbl, g)
301 do j=g%jsc-1, g%jec+1 ;
do i=g%isc-1,g%iec+1
302 call boundary_k_range(
surface, g%ke, h(i,j,:), hbl(i,j), k_top(i,j), zeta_top(i,j), k_bot(i,j), zeta_bot(i,j))
308 if (gv%Boussinesq)
then
309 h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
311 h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10
315 if (cs%ref_pres>=0.)
then
316 ref_pres(:) = cs%ref_pres
319 if (cs%continuous_reconstruction)
then
325 cs%dRdT_i(:,:,:,:) = 0.
326 cs%dRdS_i(:,:,:,:) = 0.
328 cs%stable_cell(:,:,:) = .true.
333 do k=1,g%ke ;
do j=g%jsc-1, g%jec+1 ;
do i=g%isc-1,g%iec+1
334 cs%Pint(i,j,k+1) = cs%Pint(i,j,k) + h(i,j,k)*gv%H_to_Pa
335 enddo ;
enddo ;
enddo
339 if (.not. cs%continuous_reconstruction)
then
340 do j=g%jsc-1, g%jec+1 ;
do i=g%isc-1,g%iec+1
342 cs%P_i(i,j,1,2) = h(i,j,1)*gv%H_to_Pa
344 do k=2,g%ke ;
do j=g%jsc-1, g%jec+1 ;
do i=g%isc-1,g%iec+1
345 cs%P_i(i,j,k,1) = cs%P_i(i,j,k-1,2)
346 cs%P_i(i,j,k,2) = cs%P_i(i,j,k-1,2) + h(i,j,k)*gv%H_to_Pa
347 enddo ;
enddo ;
enddo
350 do j = g%jsc-1, g%jec+1
352 do i = g%isc-1, g%iec+1
353 if (cs%continuous_reconstruction)
then
354 call interface_scalar(g%ke, h(i,j,:), t(i,j,:), cs%Tint(i,j,:), 2, h_neglect)
355 call interface_scalar(g%ke, h(i,j,:), s(i,j,:), cs%Sint(i,j,:), 2, h_neglect)
357 call build_reconstructions_1d( cs%remap_CS, g%ke, h(i,j,:), t(i,j,:), cs%ppoly_coeffs_T(i,j,:,:), &
358 cs%T_i(i,j,:,:), ppoly_r_s, imethod, h_neglect, h_neglect_edge )
359 call build_reconstructions_1d( cs%remap_CS, g%ke, h(i,j,:), s(i,j,:), cs%ppoly_coeffs_S(i,j,:,:), &
360 cs%S_i(i,j,:,:), ppoly_r_s, imethod, h_neglect, h_neglect_edge )
373 if (cs%continuous_reconstruction)
then
375 if (cs%ref_pres<0) ref_pres(:) = cs%Pint(:,j,k)
377 cs%dRdT(:,j,k), cs%dRdS(:,j,k), g%isc-1, g%iec-g%isc+3, cs%EOS)
381 if (cs%ref_pres<0) ref_pres(:) = cs%Pint(:,j,k)
384 cs%dRdT_i(:,j,k,1), cs%dRdS_i(:,j,k,1), g%isc-1, g%iec-g%isc+3, cs%EOS)
385 if (cs%ref_pres<0)
then
386 ref_pres(:) = cs%Pint(:,j,k+1)
390 cs%dRdT_i(:,j,k,2), cs%dRdS_i(:,j,k,2), g%isc-1, g%iec-g%isc+3, cs%EOS)
395 if (.not. cs%continuous_reconstruction)
then
396 do j = g%jsc-1, g%jec+1 ;
do i = g%isc-1, g%iec+1
397 call mark_unstable_cells( cs, g%ke, cs%T_i(i,j,:,:), cs%S_i(i,j,:,:), cs%P_i(i,j,:,:), cs%stable_cell(i,j,:) )
398 if (cs%interior_only)
then
399 if (.not. cs%stable_cell(i,j,k_bot(i,j))) zeta_bot(i,j) = -1.
401 do k = 1, k_bot(i,j)-1
402 cs%stable_cell(i,j,k) = .false.
420 do j = g%jsc, g%jec ;
do i = g%isc-1, g%iec
421 if (g%mask2dCu(i,j) > 0.)
then
422 if (cs%continuous_reconstruction)
then
424 cs%Pint(i,j,:), cs%Tint(i,j,:), cs%Sint(i,j,:), cs%dRdT(i,j,:), cs%dRdS(i,j,:), &
425 cs%Pint(i+1,j,:), cs%Tint(i+1,j,:), cs%Sint(i+1,j,:), cs%dRdT(i+1,j,:), cs%dRdS(i+1,j,:), &
426 cs%uPoL(i,j,:), cs%uPoR(i,j,:), cs%uKoL(i,j,:), cs%uKoR(i,j,:), cs%uhEff(i,j,:), &
427 k_bot(i,j), k_bot(i+1,j), 1.-zeta_bot(i,j), 1.-zeta_bot(i+1,j))
430 cs%P_i(i,j,:,:), h(i,j,:), cs%T_i(i,j,:,:), cs%S_i(i,j,:,:), cs%ppoly_coeffs_T(i,j,:,:), &
431 cs%ppoly_coeffs_S(i,j,:,:),cs%stable_cell(i,j,:), &
432 cs%P_i(i+1,j,:,:), h(i+1,j,:), cs%T_i(i+1,j,:,:), cs%S_i(i+1,j,:,:), cs%ppoly_coeffs_T(i+1,j,:,:), &
433 cs%ppoly_coeffs_S(i+1,j,:,:), cs%stable_cell(i+1,j,:), &
434 cs%uPoL(i,j,:), cs%uPoR(i,j,:), cs%uKoL(i,j,:), cs%uKoR(i,j,:), cs%uhEff(i,j,:), &
435 hard_fail_heff = cs%hard_fail_heff)
441 do j = g%jsc-1, g%jec ;
do i = g%isc, g%iec
442 if (g%mask2dCv(i,j) > 0.)
then
443 if (cs%continuous_reconstruction)
then
445 cs%Pint(i,j,:), cs%Tint(i,j,:), cs%Sint(i,j,:), cs%dRdT(i,j,:), cs%dRdS(i,j,:), &
446 cs%Pint(i,j+1,:), cs%Tint(i,j+1,:), cs%Sint(i,j+1,:), cs%dRdT(i,j+1,:), cs%dRdS(i,j+1,:), &
447 cs%vPoL(i,j,:), cs%vPoR(i,j,:), cs%vKoL(i,j,:), cs%vKoR(i,j,:), cs%vhEff(i,j,:), &
448 k_bot(i,j), k_bot(i,j+1), 1.-zeta_bot(i,j), 1.-zeta_bot(i,j+1))
451 cs%P_i(i,j,:,:), h(i,j,:), cs%T_i(i,j,:,:), cs%S_i(i,j,:,:), cs%ppoly_coeffs_T(i,j,:,:), &
452 cs%ppoly_coeffs_S(i,j,:,:),cs%stable_cell(i,j,:), &
453 cs%P_i(i,j+1,:,:), h(i,j+1,:), cs%T_i(i,j+1,:,:), cs%S_i(i,j+1,:,:), cs%ppoly_coeffs_T(i,j+1,:,:), &
454 cs%ppoly_coeffs_S(i,j+1,:,:), cs%stable_cell(i,j+1,:), &
455 cs%vPoL(i,j,:), cs%vPoR(i,j,:), cs%vKoL(i,j,:), cs%vKoR(i,j,:), cs%vhEff(i,j,:), &
456 hard_fail_heff = cs%hard_fail_heff)
465 if (cs%continuous_reconstruction)
then
466 do k = 1, cs%nsurf-1 ;
do j = g%jsc, g%jec ;
do i = g%isc-1, g%iec
467 if (g%mask2dCu(i,j) > 0.) cs%uhEff(i,j,k) = cs%uhEff(i,j,k) * pa_to_h
468 enddo ;
enddo ;
enddo
469 do k = 1, cs%nsurf-1 ;
do j = g%jsc-1, g%jec ;
do i = g%isc, g%iec
470 if (g%mask2dCv(i,j) > 0.) cs%vhEff(i,j,k) = cs%vhEff(i,j,k) * pa_to_h
471 enddo ;
enddo ;
enddo
474 if (cs%id_uhEff_2d>0)
then
476 do k = 1,cs%nsurf-1 ;
do j=g%jsc,g%jec ;
do i=g%isc-1,g%iec
477 heff_sum(i,j) = heff_sum(i,j) + cs%uhEff(i,j,k)
478 enddo ;
enddo ;
enddo
479 call post_data(cs%id_uhEff_2d, heff_sum, cs%diag)
481 if (cs%id_vhEff_2d>0)
then
483 do k = 1,cs%nsurf-1 ;
do j=g%jsc-1,g%jec ;
do i=g%isc,g%iec
484 heff_sum(i,j) = heff_sum(i,j) + cs%vhEff(i,j,k)
485 enddo ;
enddo ;
enddo
486 call post_data(cs%id_vhEff_2d, heff_sum, cs%diag)
495 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
496 real,
dimension(SZIB_(G),SZJ_(G)),
intent(in) :: coef_x
497 real,
dimension(SZI_(G),SZJB_(G)),
intent(in) :: coef_y
498 real,
intent(in) :: dt
505 real,
dimension(SZIB_(G),SZJ_(G),CS%nsurf-1) :: uflx
506 real,
dimension(SZI_(G),SZJB_(G),CS%nsurf-1) :: vflx
508 real,
dimension(SZI_(G),SZJ_(G),G%ke) :: tendency
509 real,
dimension(SZI_(G),SZJ_(G)) :: tendency_2d
510 real,
dimension(SZIB_(G),SZJ_(G)) :: trans_x_2d
511 real,
dimension(SZI_(G),SZJB_(G)) :: trans_y_2d
512 real,
dimension(G%ke) :: dtracer
516 integer :: i, j, k, m, ks, nk
518 real :: h_neglect, h_neglect_edge
521 h_neglect_edge = gv%m_to_H*1.0e-10
522 h_neglect = gv%m_to_H*1.0e-30
531 if (tracer%id_dfxy_conc > 0 .or. tracer%id_dfxy_cont > 0 .or. tracer%id_dfxy_cont_2d > 0 .or. &
532 tracer%id_dfx_2d > 0 .or. tracer%id_dfy_2d > 0)
then
534 tendency(:,:,:) = 0.0
541 do j = g%jsc,g%jec ;
do i = g%isc-1,g%iec
542 if (g%mask2dCu(i,j)>0.)
then
544 tracer%t(i,j,:), tracer%t(i+1,j,:), &
545 cs%uPoL(i,j,:), cs%uPoR(i,j,:), &
546 cs%uKoL(i,j,:), cs%uKoR(i,j,:), &
547 cs%uhEff(i,j,:), uflx(i,j,:), &
548 cs%continuous_reconstruction, h_neglect, cs%remap_CS, h_neglect_edge)
553 do j = g%jsc-1,g%jec ;
do i = g%isc,g%iec
554 if (g%mask2dCv(i,j)>0.)
then
556 tracer%t(i,j,:), tracer%t(i,j+1,:), &
557 cs%vPoL(i,j,:), cs%vPoR(i,j,:), &
558 cs%vKoL(i,j,:), cs%vKoR(i,j,:), &
559 cs%vhEff(i,j,:), vflx(i,j,:), &
560 cs%continuous_reconstruction, h_neglect, cs%remap_CS, h_neglect_edge)
565 do j = g%jsc,g%jec ;
do i = g%isc,g%iec
566 if (g%mask2dT(i,j)>0.)
then
571 dtracer(k) = dtracer(k) + coef_x(i,j) * uflx(i,j,ks)
572 k = cs%uKoR(i-1,j,ks)
573 dtracer(k) = dtracer(k) - coef_x(i-1,j) * uflx(i-1,j,ks)
575 dtracer(k) = dtracer(k) + coef_y(i,j) * vflx(i,j,ks)
576 k = cs%vKoR(i,j-1,ks)
577 dtracer(k) = dtracer(k) - coef_y(i,j-1) * vflx(i,j-1,ks)
580 tracer%t(i,j,k) = tracer%t(i,j,k) + dtracer(k) * &
581 ( g%IareaT(i,j) / ( h(i,j,k) + gv%H_subroundoff ) )
584 if (tracer%id_dfxy_conc > 0 .or. tracer%id_dfxy_cont > 0 .or. tracer%id_dfxy_cont_2d > 0 )
then
586 tendency(i,j,k) = dtracer(k) * g%IareaT(i,j) * idt
595 if (tracer%id_dfx_2d > 0)
then
596 do j = g%jsc,g%jec ;
do i = g%isc-1,g%iec
598 if (g%mask2dCu(i,j)>0.)
then
600 trans_x_2d(i,j) = trans_x_2d(i,j) - coef_x(i,j) * uflx(i,j,ks)
602 trans_x_2d(i,j) = trans_x_2d(i,j) * idt
605 call post_data(tracer%id_dfx_2d, trans_x_2d(:,:), cs%diag)
610 if (tracer%id_dfy_2d > 0)
then
611 do j = g%jsc-1,g%jec ;
do i = g%isc,g%iec
613 if (g%mask2dCv(i,j)>0.)
then
615 trans_y_2d(i,j) = trans_y_2d(i,j) - coef_y(i,j) * vflx(i,j,ks)
617 trans_y_2d(i,j) = trans_y_2d(i,j) * idt
620 call post_data(tracer%id_dfy_2d, trans_y_2d(:,:), cs%diag)
624 if (tracer%id_dfxy_cont > 0)
then
625 call post_data(tracer%id_dfxy_cont, tendency(:,:,:), cs%diag)
629 if (tracer%id_dfxy_cont_2d > 0)
then
630 tendency_2d(:,:) = 0.
631 do j = g%jsc,g%jec ;
do i = g%isc,g%iec
633 tendency_2d(i,j) = tendency_2d(i,j) + tendency(i,j,k)
636 call post_data(tracer%id_dfxy_cont_2d, tendency_2d(:,:), cs%diag)
642 if (tracer%id_dfxy_conc > 0)
then
643 do k = 1, gv%ke ;
do j = g%jsc,g%jec ;
do i = g%isc,g%iec
644 tendency(i,j,k) = tendency(i,j,k) / ( h(i,j,k) + gv%H_subroundoff )
645 enddo ;
enddo ;
enddo
646 call post_data(tracer%id_dfxy_conc, tendency, cs%diag)
654 integer,
intent(in) :: nk
655 real,
dimension(nk),
intent(in) :: h
656 real,
dimension(nk),
intent(in) :: S
657 real,
dimension(nk+1),
intent(inout) :: Si
658 integer,
intent(in) :: i_method
660 real,
intent(in) :: h_neglect
662 integer :: k, km2, kp1
663 real,
dimension(nk) :: diff
667 si(1) = s(1) - 0.5 * diff(1)
668 if (i_method==1)
then
672 sa = s(k-1) + 0.5 * diff(k-1)
673 sb = s(k) - 0.5 * diff(k)
674 si(k) = 0.5 * ( sa + sb )
676 elseif (i_method==2)
then
682 si(k) =
ppm_edge(h(km2), h(k-1), h(k), h(kp1), s(k-1), s(k), diff(k-1), diff(k), h_neglect)
685 si(nk+1) = s(nk) + 0.5 * diff(nk)
691 real function ppm_edge(hkm1, hk, hkp1, hkp2, Ak, Akp1, Pk, Pkp1, h_neglect)
692 real,
intent(in) :: hkm1
693 real,
intent(in) :: hk
694 real,
intent(in) :: hkp1
695 real,
intent(in) :: hkp2
696 real,
intent(in) :: ak
697 real,
intent(in) :: akp1
698 real,
intent(in) :: pk
699 real,
intent(in) :: pkp1
700 real,
intent(in) :: h_neglect
703 real :: r_hk_hkp1, r_2hk_hkp1, r_hk_2hkp1, f1, f2, f3, f4
705 r_hk_hkp1 = hk + hkp1
706 if (r_hk_hkp1 <= 0.)
then
710 r_hk_hkp1 = 1. / r_hk_hkp1
712 ppm_edge = ak + ( hk * r_hk_hkp1 ) * ( akp1 - ak )
714 ppm_edge = akp1 + ( hkp1 * r_hk_hkp1 ) * ( ak - akp1 )
717 r_2hk_hkp1 = 1. / ( ( 2. * hk + hkp1 ) + h_neglect )
718 r_hk_2hkp1 = 1. / ( ( hk + 2. * hkp1 ) + h_neglect )
719 f1 = 1./ ( ( hk + hkp1) + ( hkm1 + hkp2 ) )
720 f2 = 2. * ( hkp1 * hk ) * r_hk_hkp1 * &
721 ( ( hkm1 + hk ) * r_2hk_hkp1 - ( hkp2 + hkp1 ) * r_hk_2hkp1 )
722 f3 = hk * ( hkm1 + hk ) * r_2hk_hkp1
723 f4 = hkp1 * ( hkp1 + hkp2 ) * r_hk_2hkp1
731 real function ppm_ave(xL, xR, aL, aR, aMean)
732 real,
intent(in) :: xl
733 real,
intent(in) :: xr
734 real,
intent(in) :: al
735 real,
intent(in) :: ar
736 real,
intent(in) :: amean
739 real :: dx, xave, a6, a6o3
742 xave = 0.5 * ( xr + xl )
743 a6o3 = 2. * amean - ( al + ar )
747 stop
'ppm_ave: dx<0 should not happend!'
749 stop
'ppm_ave: dx>1 should not happend!'
751 ppm_ave = al + ( ar - al ) * xr + a6 * xr * ( 1. - xr )
753 ppm_ave = ( al + xave * ( ( ar - al ) + a6 ) ) - a6o3 * ( xr**2 + xr * xl + xl**2 )
759 real,
intent(in) :: a
760 real,
intent(in) :: x
769 subroutine plm_diff(nk, h, S, c_method, b_method, diff)
770 integer,
intent(in) :: nk
771 real,
dimension(nk),
intent(in) :: h
772 real,
dimension(nk),
intent(in) :: S
773 integer,
intent(in) :: c_method
774 integer,
intent(in) :: b_method
775 real,
dimension(nk),
intent(inout) :: diff
785 real :: hkm1, hk, hkp1, Skm1, Sk, Skp1, diff_l, diff_r, diff_c
792 if ( ( hkp1 + hk ) * ( hkm1 + hk ) > 0.)
then
796 if (c_method==1)
then
798 if ( hk + 0.5 * (hkm1 + hkp1) /= 0. )
then
799 diff_c = ( skp1 - skm1 ) * ( hk / ( hk + 0.5 * (hkm1 + hkp1) ) )
803 elseif (c_method==2)
then
805 diff_c =
fv_diff(hkm1, hk, hkp1, skm1, sk, skp1)
806 elseif (c_method==3)
then
808 diff_c = hk *
fvlsq_slope(hkm1, hk, hkp1, skm1, sk, skp1)
811 diff_l = 2. * ( sk - skm1 )
812 diff_r = 2. * ( skp1 - sk )
813 if (
signum(1., diff_l) *
signum(1., diff_r) <= 0. )
then
816 diff(k) = sign( min( abs(diff_l), abs(diff_c), abs(diff_r) ), diff_c )
822 if (b_method==1)
then
825 elseif (b_method==2)
then
826 diff(1) = ( s(2) - s(1) ) * 2. * ( h(1) / ( h(1) + h(2) ) )
827 diff(nk) = s(nk) - s(nk-1) * 2. * ( h(nk) / ( h(nk-1) + h(nk) ) )
836 real function fv_diff(hkm1, hk, hkp1, Skm1, Sk, Skp1)
837 real,
intent(in) :: hkm1
838 real,
intent(in) :: hk
839 real,
intent(in) :: hkp1
840 real,
intent(in) :: skm1
841 real,
intent(in) :: sk
842 real,
intent(in) :: skp1
845 real :: h_sum, hp, hm
847 h_sum = ( hkm1 + hkp1 ) + hk
848 if (h_sum /= 0.) h_sum = 1./ h_sum
850 if (hm /= 0.) hm = 1./ hm
852 if (hp /= 0.) hp = 1./ hp
854 ( ( 2. * hkm1 + hk ) * hp * ( skp1 - sk ) &
855 + ( 2. * hkp1 + hk ) * hm * ( sk - skm1 ) )
862 real function fvlsq_slope(hkm1, hk, hkp1, Skm1, Sk, Skp1)
863 real,
intent(in) :: hkm1
864 real,
intent(in) :: hk
865 real,
intent(in) :: hkp1
866 real,
intent(in) :: skm1
867 real,
intent(in) :: sk
868 real,
intent(in) :: skp1
872 real :: h_sum, hx_sum, hxsq_sum, hxy_sum, hy_sum, det
874 xkm1 = -0.5 * ( hk + hkm1 )
875 xkp1 = 0.5 * ( hk + hkp1 )
876 h_sum = ( hkm1 + hkp1 ) + hk
877 hx_sum = hkm1*xkm1 + hkp1*xkp1
878 hxsq_sum = hkm1*(xkm1**2) + hkp1*(xkp1**2)
879 hxy_sum = hkm1*xkm1*skm1 + hkp1*xkp1*skp1
880 hy_sum = ( hkm1*skm1 + hkp1*skp1 ) + hk*sk
881 det = h_sum * hxsq_sum - hx_sum**2
884 fvlsq_slope = ( h_sum * hxy_sum - hx_sum*hy_sum ) / det
893 dRdTr, dRdSr, PoL, PoR, KoL, KoR, hEff, bl_kl, bl_kr, bl_zl, bl_zr)
894 integer,
intent(in) :: nk
895 real,
dimension(nk+1),
intent(in) :: Pl
896 real,
dimension(nk+1),
intent(in) :: Tl
897 real,
dimension(nk+1),
intent(in) :: Sl
898 real,
dimension(nk+1),
intent(in) :: dRdTl
899 real,
dimension(nk+1),
intent(in) :: dRdSl
900 real,
dimension(nk+1),
intent(in) :: Pr
901 real,
dimension(nk+1),
intent(in) :: Tr
902 real,
dimension(nk+1),
intent(in) :: Sr
903 real,
dimension(nk+1),
intent(in) :: dRdTr
904 real,
dimension(nk+1),
intent(in) :: dRdSr
905 real,
dimension(2*nk+2),
intent(inout) :: PoL
907 real,
dimension(2*nk+2),
intent(inout) :: PoR
909 integer,
dimension(2*nk+2),
intent(inout) :: KoL
910 integer,
dimension(2*nk+2),
intent(inout) :: KoR
911 real,
dimension(2*nk+1),
intent(inout) :: hEff
912 integer,
optional,
intent(in) :: bl_kl
913 integer,
optional,
intent(in) :: bl_kr
914 real,
optional,
intent(in) :: bl_zl
915 real,
optional,
intent(in) :: bl_zr
923 logical :: searching_left_column
924 logical :: searching_right_column
925 logical :: reached_bottom
926 integer :: krm1, klm1
927 real :: dRho, dRhoTop, dRhoBot, hL, hR
928 integer :: lastK_left, lastK_right
929 real :: lastP_left, lastP_right
930 logical :: interior_limit
941 reached_bottom = .false.
944 interior_limit =
PRESENT(bl_kl) .and.
PRESENT(bl_kr) .and.
PRESENT(bl_zr) .and.
PRESENT(bl_zl)
947 neutral_surfaces:
do k_surface = 1, ns
949 if (klm1>nk) stop
'find_neutral_surface_positions(): klm1 went out of bounds!'
951 if (krm1>nk) stop
'find_neutral_surface_positions(): krm1 went out of bounds!'
954 drho = 0.5 * ( ( drdtr(kr) + drdtl(kl) ) * ( tr(kr) - tl(kl) ) &
955 + ( drdsr(kr) + drdsl(kl) ) * ( sr(kr) - sl(kl) ) )
957 if (.not. reached_bottom)
then
959 searching_left_column = .true.
960 searching_right_column = .false.
961 elseif (drho > 0.)
then
962 searching_right_column = .true.
963 searching_left_column = .false.
965 if (kl + kr == 2)
then
966 searching_left_column = .true.
967 searching_right_column = .false.
969 searching_left_column = .not. searching_left_column
970 searching_right_column = .not. searching_right_column
975 if (searching_left_column)
then
978 drhotop = 0.5 * ( ( drdtl(klm1) + drdtr(kr) ) * ( tl(klm1) - tr(kr) ) &
979 + ( drdsl(klm1) + drdsr(kr) ) * ( sl(klm1) - sr(kr) ) )
981 drhobot = 0.5 * ( ( drdtl(klm1+1) + drdtr(kr) ) * ( tl(klm1+1) - tr(kr) ) &
982 + ( drdsl(klm1+1) + drdsr(kr) ) * ( sl(klm1+1) - sr(kr) ) )
986 if (drhotop > 0. .or. kr+kl==2)
then
988 elseif (drhotop >= drhobot)
then
995 if (pol(k_surface)>=1. .and. klm1<nk)
then
997 pol(k_surface) = pol(k_surface) - 1.
999 if (real(klm1-lastk_left)+(pol(k_surface)-lastp_left)<0.)
then
1000 pol(k_surface) = lastp_left
1003 kol(k_surface) = klm1
1014 reached_bottom = .true.
1015 searching_right_column = .true.
1016 searching_left_column = .false.
1018 elseif (searching_right_column)
then
1021 drhotop = 0.5 * ( ( drdtr(krm1) + drdtl(kl) ) * ( tr(krm1) - tl(kl) ) &
1022 + ( drdsr(krm1) + drdsl(kl) ) * ( sr(krm1) - sl(kl) ) )
1024 drhobot = 0.5 * ( ( drdtr(krm1+1) + drdtl(kl) ) * ( tr(krm1+1) - tl(kl) ) &
1025 + ( drdsr(krm1+1) + drdsl(kl) ) * ( sr(krm1+1) - sl(kl) ) )
1029 if (drhotop >= 0. .or. kr+kl==2)
then
1031 elseif (drhotop >= drhobot)
then
1038 if (por(k_surface)>=1. .and. krm1<nk)
then
1040 por(k_surface) = por(k_surface) - 1.
1042 if (real(krm1-lastk_right)+(por(k_surface)-lastp_right)<0.)
then
1043 por(k_surface) = lastp_right
1046 kor(k_surface) = krm1
1057 reached_bottom = .true.
1058 searching_right_column = .false.
1059 searching_left_column = .true.
1064 if (interior_limit)
then
1065 if (kol(k_surface)<=bl_kl)
then
1066 kol(k_surface) = bl_kl
1067 if (pol(k_surface)<bl_zl)
then
1068 pol(k_surface) = bl_zl
1071 if (kor(k_surface)<=bl_kr)
then
1072 kor(k_surface) = bl_kr
1073 if (por(k_surface)<bl_zr)
then
1074 por(k_surface) = bl_zr
1079 lastk_left = kol(k_surface) ; lastp_left = pol(k_surface)
1080 lastk_right = kor(k_surface) ; lastp_right = por(k_surface)
1084 if (k_surface>1)
then
1085 hl =
absolute_position(nk,ns,pl,kol,pol,k_surface) -
absolute_position(nk,ns,pl,kol,pol,k_surface-1)
1086 hr =
absolute_position(nk,ns,pr,kor,por,k_surface) -
absolute_position(nk,ns,pr,kor,por,k_surface-1)
1087 if ( hl + hr > 0.)
then
1088 heff(k_surface-1) = 2. * hl * hr / ( hl + hr )
1090 heff(k_surface-1) = 0.
1094 enddo neutral_surfaces
1101 real,
intent(in) :: drhoneg
1102 real,
intent(in) :: pneg
1103 real,
intent(in) :: drhopos
1104 real,
intent(in) :: ppos
1107 stop
'interpolate_for_nondim_position: Houston, we have a problem! Ppos<Pneg'
1108 elseif (drhoneg>drhopos)
then
1109 write(0,*)
'dRhoNeg, Pneg, dRhoPos, Ppos=',drhoneg, pneg, drhopos, ppos
1110 elseif (drhoneg>drhopos)
then
1111 stop
'interpolate_for_nondim_position: Houston, we have a problem! dRhoNeg>dRhoPos'
1113 if (ppos<=pneg)
then
1115 elseif ( drhopos - drhoneg > 0. )
then
1117 elseif ( drhopos - drhoneg == 0)
then
1118 if (drhoneg>0.)
then
1120 elseif (drhoneg<0.)
then
1129 stop
'interpolate_for_nondim_position: Houston, we have a problem! Pint < Pneg'
1131 stop
'interpolate_for_nondim_position: Houston, we have a problem! Pint > Ppos'
1138 Pres_r, hcol_r, Tr, Sr, ppoly_T_r, ppoly_S_r, stable_r,&
1139 PoL, PoR, KoL, KoR, hEff, zeta_bot_L, zeta_bot_R, &
1140 k_bot_L, k_bot_R, hard_fail_heff)
1143 integer,
intent(in) :: nk
1144 real,
dimension(nk,2),
intent(in) :: Pres_l
1145 real,
dimension(nk),
intent(in) :: hcol_l
1146 real,
dimension(nk,2),
intent(in) :: Tl
1147 real,
dimension(nk,2),
intent(in) :: Sl
1148 real,
dimension(:,:),
intent(in) :: ppoly_T_l
1149 real,
dimension(:,:),
intent(in) :: ppoly_S_l
1150 logical,
dimension(nk),
intent(in) :: stable_l
1151 real,
dimension(nk,2),
intent(in) :: Pres_r
1152 real,
dimension(nk),
intent(in) :: hcol_r
1153 real,
dimension(nk,2),
intent(in) :: Tr
1154 real,
dimension(nk,2),
intent(in) :: Sr
1155 real,
dimension(:,:),
intent(in) :: ppoly_T_r
1156 real,
dimension(:,:),
intent(in) :: ppoly_S_r
1157 logical,
dimension(nk),
intent(in) :: stable_r
1158 real,
dimension(4*nk),
intent(inout) :: PoL
1160 real,
dimension(4*nk),
intent(inout) :: PoR
1162 integer,
dimension(4*nk),
intent(inout) :: KoL
1163 integer,
dimension(4*nk),
intent(inout) :: KoR
1164 real,
dimension(4*nk-1),
intent(inout) :: hEff
1165 real,
optional,
intent(in) :: zeta_bot_L
1167 real,
optional,
intent(in) :: zeta_bot_R
1170 integer,
optional,
intent(in) :: k_bot_L
1171 integer,
optional,
intent(in) :: k_bot_R
1172 logical,
optional,
intent(in) :: hard_fail_heff
1176 integer :: k_surface
1177 integer :: kl_left, kl_right
1178 integer :: ki_left, ki_right
1179 logical :: searching_left_column
1180 logical :: searching_right_column
1181 logical :: reached_bottom
1182 logical :: search_layer
1183 logical :: fail_heff
1184 real :: dRho, dRhoTop, dRhoBot, hL, hR
1186 real :: dRdT_from_top, dRdS_from_top
1187 real :: dRdT_from_bot, dRdS_from_bot
1188 real :: dRdT_to_top, dRdS_to_top
1189 real :: dRdT_to_bot, dRdS_to_bot
1190 real :: T_ref, S_ref, P_ref, P_top, P_bot
1191 real :: lastP_left, lastP_right
1192 integer :: k_init_L, k_init_R
1193 real :: p_init_L, p_init_R
1202 reached_bottom = .false.
1203 searching_left_column = .false.
1204 searching_right_column = .false.
1207 if (
PRESENT(hard_fail_heff)) fail_heff = hard_fail_heff
1209 if (
PRESENT(k_bot_l) .and.
PRESENT(k_bot_r) .and.
PRESENT(zeta_bot_l) .and.
PRESENT(zeta_bot_r))
then
1210 k_init_l = k_bot_l; k_init_r = k_bot_r
1211 p_init_l = zeta_bot_l; p_init_r = zeta_bot_r
1212 lastp_left = zeta_bot_l; lastp_right = zeta_bot_r
1213 kl_left = k_bot_l; kl_right = k_bot_r
1215 k_init_l = 1 ; k_init_r = 1
1216 p_init_l = 0. ; p_init_r = 0.
1219 neutral_surfaces:
do k_surface = 1, ns
1221 if (k_surface == ns)
then
1227 elseif (.not. stable_l(kl_left))
then
1228 if (k_surface > 1)
then
1229 pol(k_surface) = ki_left - 1
1230 kol(k_surface) = kl_left
1231 por(k_surface) = por(k_surface-1)
1232 kor(k_surface) = kor(k_surface-1)
1234 por(k_surface) = p_init_r
1235 kor(k_surface) = k_init_r
1236 pol(k_surface) = p_init_l
1237 kol(k_surface) = k_init_l
1239 call increment_interface(nk, kl_left, ki_left, reached_bottom, searching_left_column, searching_right_column)
1240 searching_left_column = .true.
1241 searching_right_column = .false.
1242 elseif (.not. stable_r(kl_right))
then
1243 if (k_surface > 1)
then
1244 por(k_surface) = ki_right - 1
1245 kor(k_surface) = kl_right
1246 pol(k_surface) = pol(k_surface-1)
1247 kol(k_surface) = kol(k_surface-1)
1254 call increment_interface(nk, kl_right, ki_right, reached_bottom, searching_right_column, searching_left_column)
1255 searching_left_column = .false.
1256 searching_right_column = .true.
1262 tr(kl_right, ki_right), sr(kl_right, ki_right), pres_r(kl_right,ki_right), &
1263 tl(kl_left, ki_left), sl(kl_left, ki_left) , pres_l(kl_left,ki_left), &
1265 if (cs%debug)
write(*,
'(A,I2,A,E12.4,A,I2,A,I2,A,I2,A,I2)')
"k_surface=",k_surface,
" dRho=",drho, &
1266 "kl_left=",kl_left,
" ki_left=",ki_left,
" kl_right=",kl_right,
" ki_right=",ki_right
1268 if (.not. reached_bottom)
then
1270 searching_left_column = .true.
1271 searching_right_column = .false.
1272 elseif (drho > 0.)
then
1273 searching_left_column = .false.
1274 searching_right_column = .true.
1276 if ( ( kl_left + kl_right == 2 ) .and. (ki_left + ki_right == 2) )
then
1277 searching_left_column = .true.
1278 searching_right_column = .false.
1280 searching_left_column = .not. searching_left_column
1281 searching_right_column = .not. searching_right_column
1285 if (searching_left_column)
then
1287 por(k_surface) = ki_right - 1.
1288 kor(k_surface) = kl_right
1290 tr(kl_right, ki_right), sr(kl_right, ki_right), pres_r(kl_right, ki_right), &
1291 tl(kl_left,1), sl(kl_left,1), pres_l(kl_left,1), &
1292 tl(kl_left,2), sl(kl_left,2), pres_l(kl_left,2), &
1293 ppoly_t_l(kl_left,:), ppoly_s_l(kl_left,:))
1294 kol(k_surface) = kl_left
1297 write(*,
'(A,I2)')
"Searching left layer ", kl_left
1298 write(*,
'(A,I2,X,I2)')
"Searching from right: ", kl_right, ki_right
1299 write(*,*)
"Temp/Salt Reference: ", tr(kl_right,ki_right), sr(kl_right,ki_right)
1300 write(*,*)
"Temp/Salt Top L: ", tl(kl_left,1), sl(kl_left,1)
1301 write(*,*)
"Temp/Salt Bot L: ", tl(kl_left,2), sl(kl_left,2)
1303 call increment_interface(nk, kl_right, ki_right, reached_bottom, searching_right_column, searching_left_column)
1304 lastp_left = pol(k_surface)
1306 if ( kl_right == (kor(k_surface) + 1) ) lastp_right = 0.
1308 elseif (searching_right_column)
then
1310 pol(k_surface) = ki_left - 1.
1311 kol(k_surface) = kl_left
1313 tl(kl_left, ki_left), sl(kl_left, ki_left), pres_l(kl_left, ki_left), &
1314 tr(kl_right,1), sr(kl_right,1), pres_r(kl_right,1), &
1315 tr(kl_right,2), sr(kl_right,2), pres_r(kl_right,2), &
1316 ppoly_t_r(kl_right,:), ppoly_s_r(kl_right,:))
1317 kor(k_surface) = kl_right
1320 write(*,
'(A,I2)')
"Searching right layer ", kl_right
1321 write(*,
'(A,I2,X,I2)')
"Searching from left: ", kl_left, ki_left
1322 write(*,*)
"Temp/Salt Reference: ", tl(kl_left,ki_left), sl(kl_left,ki_left)
1323 write(*,*)
"Temp/Salt Top L: ", tr(kl_right,1), sr(kl_right,1)
1324 write(*,*)
"Temp/Salt Bot L: ", tr(kl_right,2), sr(kl_right,2)
1326 call increment_interface(nk, kl_left, ki_left, reached_bottom, searching_left_column, searching_right_column)
1327 lastp_right = por(k_surface)
1329 if ( kl_left == (kol(k_surface) + 1) ) lastp_left = 0.
1333 if (cs%debug)
write(*,
'(A,I3,A,ES16.6,A,I2,A,ES16.6)')
"KoL:", kol(k_surface),
" PoL:", pol(k_surface), &
1334 " KoR:", kor(k_surface),
" PoR:", por(k_surface)
1337 if (k_surface>1)
then
1338 if ( kol(k_surface) == kol(k_surface-1) .and. kor(k_surface) == kor(k_surface-1) )
then
1339 hl = (pol(k_surface) - pol(k_surface-1))*hcol_l(kol(k_surface))
1340 hr = (por(k_surface) - por(k_surface-1))*hcol_r(kor(k_surface))
1341 if (hl < 0. .or. hr < 0.)
then
1343 call mom_error(fatal,
"Negative thicknesses in neutral diffusion")
1345 if (searching_left_column)
then
1346 pol(k_surface) = pol(k_surface-1)
1347 kol(k_surface) = kol(k_surface-1)
1348 elseif (searching_right_column)
then
1349 por(k_surface) = por(k_surface-1)
1350 kor(k_surface) = kor(k_surface-1)
1353 elseif ( hl + hr == 0. )
then
1354 heff(k_surface-1) = 0.
1356 heff(k_surface-1) = 2. * ( (hl * hr) / ( hl + hr ) )
1357 if ( kol(k_surface) /= kol(k_surface-1) )
then
1358 call mom_error(fatal,
"Neutral sublayer spans multiple layers")
1360 if ( kor(k_surface) /= kor(k_surface-1) )
then
1361 call mom_error(fatal,
"Neutral sublayer spans multiple layers")
1365 heff(k_surface-1) = 0.
1368 enddo neutral_surfaces
1374 integer,
intent(in) :: nk
1375 real,
dimension(nk,2),
intent(in) :: T
1376 real,
dimension(nk,2),
intent(in) :: S
1377 real,
dimension(nk,2),
intent(in) :: P
1378 logical,
dimension(nk),
intent( out) :: stable_cell
1380 integer :: k, first_stable, prev_stable
1385 t(k,1), s(k,1), max(p(k,1),cs%ref_pres), delta_rho )
1386 stable_cell(k) = delta_rho > 0.
1391 real function search_other_column(CS, ksurf, pos_last, T_from, S_from, P_from, T_top, S_top, P_top, &
1392 T_bot, S_bot, P_bot, T_poly, S_poly )
result(pos)
1394 integer,
intent(in ) :: ksurf
1395 real,
intent(in ) :: pos_last
1397 real,
intent(in ) :: t_from
1398 real,
intent(in ) :: s_from
1399 real,
intent(in ) :: p_from
1400 real,
intent(in ) :: t_top
1401 real,
intent(in ) :: s_top
1402 real,
intent(in ) :: p_top
1403 real,
intent(in ) :: t_bot
1404 real,
intent(in ) :: s_bot
1405 real,
intent(in ) :: p_bot
1406 real,
dimension(:),
intent(in ) :: t_poly
1407 real,
dimension(:),
intent(in ) :: s_poly
1409 real :: drhotop, drhobot
1410 real :: drdt_top, drds_top, drdt_bot, drds_bot
1411 real :: drdt_from, drds_from
1415 if (cs%neutral_pos_method == 1 .or. cs%neutral_pos_method == 3)
then
1418 elseif (cs%neutral_pos_method == 2)
then
1420 drdt_top, drds_top, drdt_from, drds_from)
1422 drdt_bot, drds_bot, drdt_from, drds_from)
1426 if ( (drhotop > 0.) .or. (ksurf == 1) )
then
1428 elseif ( drhotop > drhobot )
then
1430 elseif ( drhotop < 0. .and. drhobot < 0.)
then
1432 elseif ( drhotop == 0. .and. drhobot == 0. )
then
1434 elseif ( drhobot == 0. )
then
1436 elseif ( drhotop == 0. )
then
1445 if (cs%neutral_pos_method==1)
then
1449 elseif (cs%neutral_pos_method == 2)
then
1451 p_top, drdt_top, drds_top, &
1452 p_bot, drdt_bot, drds_bot, t_poly, s_poly )
1453 elseif (cs%neutral_pos_method == 3)
then
1460 subroutine increment_interface(nk, kl, ki, reached_bottom, searching_this_column, searching_other_column)
1461 integer,
intent(in ) :: nk
1462 integer,
intent(inout) :: kl
1463 integer,
intent(inout) :: ki
1464 logical,
intent(inout) :: reached_bottom
1465 logical,
intent(inout) :: searching_this_column
1466 logical,
intent(inout) :: searching_other_column
1469 reached_bottom = .false.
1471 if ((ki == 2) .and. (kl < nk) )
then
1474 elseif ((kl == nk) .and. (ki==2))
then
1475 reached_bottom = .true.
1476 searching_this_column = .false.
1477 searching_other_column = .true.
1482 call mom_error(fatal,
"Unanticipated eventuality in increment_interface")
1495 P_top, dRdT_top, dRdS_top, &
1496 P_bot, dRdT_bot, dRdS_bot, ppoly_T, ppoly_S )
result( z )
1498 real,
intent(in) :: z0
1499 real,
intent(in) :: t_ref
1500 real,
intent(in) :: s_ref
1501 real,
intent(in) :: p_ref
1502 real,
intent(in) :: drdt_ref
1503 real,
intent(in) :: drds_ref
1504 real,
intent(in) :: p_top
1505 real,
intent(in) :: drdt_top
1506 real,
intent(in) :: drds_top
1507 real,
intent(in) :: p_bot
1508 real,
intent(in) :: drdt_bot
1509 real,
intent(in) :: drds_bot
1510 real,
dimension(:),
intent(in) :: ppoly_t
1512 real,
dimension(:),
intent(in) :: ppoly_s
1516 real :: drdt_diff, drds_diff
1517 real :: drho, drho_dz, drdt_z, drds_z, t_z, s_z, deltat, deltas, deltap, dt_dz, ds_dz
1518 real :: drho_min, drho_max, ztest, zmin, zmax, drdt_sum, drds_sum, dz, p_z, dp_dz
1522 real :: t_top, t_bot, s_top, s_bot
1524 nterm =
SIZE(ppoly_t)
1527 drdt_diff = drdt_bot - drdt_top
1528 drds_diff = drds_bot - drds_top
1530 dp_dz = p_bot - p_top
1538 drdt_z = a1*drdt_top + a2*drdt_bot
1539 drds_z = a1*drds_top + a2*drds_bot
1540 p_z = a1*p_top + a2*p_bot
1542 t_ref, s_ref, p_ref, drdt_ref, drds_ref)
1547 t_ref, s_ref, p_ref, drdt_ref, drds_ref)
1549 if (drho_min >= 0.)
then
1552 elseif (drho_max == 0.)
then
1556 if ( sign(1.,drho_min) == sign(1.,drho_max) )
then
1557 call mom_error(fatal,
"drho_min is the same sign as dhro_max")
1562 do iter = 1, cs%max_iter
1566 drdt_z = a1*drdt_top + a2*drdt_bot
1567 drds_z = a1*drds_top + a2*drds_bot
1570 p_z = a1*p_top + a2*p_bot
1571 deltat = t_z - t_ref
1572 deltas = s_z - s_ref
1573 deltap = p_z - p_ref
1574 drdt_sum = drdt_ref + drdt_z
1575 drds_sum = drds_ref + drds_z
1577 t_ref, s_ref, p_ref, drdt_ref, drds_ref)
1580 if (abs(drho) <= cs%drho_tol)
exit
1582 if (drho < 0. .and. drho > drho_min)
then
1585 elseif (drho > 0. .and. drho < drho_max)
then
1593 drho_dz = 0.5*( (drdt_diff*deltat + drdt_sum*dt_dz) + (drds_diff*deltas + drds_sum*ds_dz) )
1595 ztest = z - drho/drho_dz
1597 if (ztest < zmin .or. ztest > zmax)
then
1598 if ( drho < 0. )
then
1599 ztest = 0.5*(z + zmax)
1601 ztest = 0.5*(zmin + z)
1606 if ( abs(z-ztest) <= cs%x_tol )
exit
1615 function find_neutral_pos_full( CS, z0, T_ref, S_ref, P_ref, P_top, P_bot, ppoly_T, ppoly_S )
result( z )
1617 real,
intent(in) :: z0
1618 real,
intent(in) :: t_ref
1619 real,
intent(in) :: s_ref
1620 real,
intent(in) :: p_ref
1621 real,
intent(in) :: p_top
1622 real,
intent(in) :: p_bot
1623 real,
dimension(:),
intent(in) :: ppoly_t
1625 real,
dimension(:),
intent(in) :: ppoly_s
1632 real :: drho_a, drho_b, drho_c
1633 real :: a, b, c, ta, tb, tc, sa, sb, sc, pa, pb, pc
1639 nterm =
SIZE(ppoly_t)
1644 pb = p_top*(1.-b) + p_bot*b
1653 if (drho_b >= 0.)
then
1656 elseif (drho_c == 0.)
then
1660 if ( sign(1.,drho_b) == sign(1.,drho_c) )
then
1665 do iter = 1, cs%max_iter
1667 a = (drho_b*c - drho_c*b)/(drho_b-drho_c)
1670 pa = p_top*(1.-a) + p_bot*a
1672 if (abs(drho_a) < cs%drho_tol)
then
1677 if (drho_a*drho_c > 0.)
then
1678 if ( abs(a-c)<cs%x_tol)
then
1682 c = a ; drho_c = drho_a;
1683 if (side == -1) drho_b = 0.5*drho_b
1685 elseif ( drho_b*drho_a > 0 )
then
1686 if ( abs(a-b)<cs%x_tol)
then
1690 b = a ; drho_b = drho_a
1691 if (side == 1) drho_c = 0.5*drho_c
1705 drdt1_out, drds1_out, drdt2_out, drds2_out )
1707 real,
intent(in ) :: T1
1708 real,
intent(in ) :: S1
1709 real,
intent(in ) :: p1_in
1710 real,
intent(in ) :: T2
1711 real,
intent(in ) :: S2
1712 real,
intent(in ) :: p2_in
1713 real,
intent( out) :: drho
1714 real,
optional,
intent( out) :: dRdT1_out
1715 real,
optional,
intent( out) :: dRdS1_out
1716 real,
optional,
intent( out) :: dRdT2_out
1717 real,
optional,
intent( out) :: dRdS2_out
1719 real :: rho1, rho2, p1, p2, pmid
1720 real :: drdt1, drdt2, drds1, drds2, drdp1, drdp2, rho_dummy
1723 if (cs%ref_pres > 0.)
then
1732 if (trim(cs%delta_rho_form) ==
'full')
then
1733 pmid = 0.5 * (p1 + p2)
1738 elseif (trim(cs%delta_rho_form) ==
'mid_pressure')
then
1739 pmid = 0.5 * (p1 + p2)
1740 if (cs%ref_pres>=0) pmid = cs%ref_pres
1744 elseif (trim(cs%delta_rho_form) ==
'local_pressure')
then
1749 call mom_error(fatal,
"delta_rho_form is not recognized")
1752 if (
PRESENT(drdt1_out)) drdt1_out = drdt1
1753 if (
PRESENT(drds1_out)) drds1_out = drds1
1754 if (
PRESENT(drdt2_out)) drdt2_out = drdt2
1755 if (
PRESENT(drds2_out)) drds2_out = drds2
1764 T2, S2, P2, dRdT2, dRdS2 )
result (drho)
1778 drho = 0.5 * ( (drdt1+drdt2)*(t1-t2) + (drds1+drds2)*(s1-s2))
1783 integer,
intent(in) :: n
1784 integer,
intent(in) :: ns
1785 real,
intent(in) :: pint(n+1)
1786 integer,
intent(in) :: karr(ns)
1787 real,
intent(in) :: nparr(ns)
1788 integer,
intent(in) :: k_surface
1793 if (k>n) stop
'absolute_position: k>nk is out of bounds!'
1800 integer,
intent(in) :: n
1801 integer,
intent(in) :: ns
1802 real,
intent(in) :: pint(n+1)
1803 integer,
intent(in) :: karr(ns)
1804 real,
intent(in) :: nparr(ns)
1809 integer :: k_surface, k
1811 do k_surface = 1, ns
1818 subroutine neutral_surface_flux(nk, nsurf, deg, hl, hr, Tl, Tr, PiL, PiR, KoL, KoR, &
1819 hEff, Flx, continuous, h_neglect, remap_CS, h_neglect_edge)
1820 integer,
intent(in) :: nk
1821 integer,
intent(in) :: nsurf
1822 integer,
intent(in) :: deg
1823 real,
dimension(nk),
intent(in) :: hl
1824 real,
dimension(nk),
intent(in) :: hr
1825 real,
dimension(nk),
intent(in) :: Tl
1826 real,
dimension(nk),
intent(in) :: Tr
1827 real,
dimension(nsurf),
intent(in) :: PiL
1829 real,
dimension(nsurf),
intent(in) :: PiR
1831 integer,
dimension(nsurf),
intent(in) :: KoL
1832 integer,
dimension(nsurf),
intent(in) :: KoR
1833 real,
dimension(nsurf-1),
intent(in) :: hEff
1834 real,
dimension(nsurf-1),
intent(inout) :: Flx
1835 logical,
intent(in) :: continuous
1836 real,
intent(in) :: h_neglect
1841 real,
optional,
intent(in) :: h_neglect_edge
1845 integer :: k_sublayer, klb, klt, krb, krt, k
1846 real :: T_right_top, T_right_bottom, T_right_layer, T_right_sub, T_right_top_int, T_right_bot_int
1847 real :: T_left_top, T_left_bottom, T_left_layer, T_left_sub, T_left_top_int, T_left_bot_int
1848 real :: dT_top, dT_bottom, dT_layer, dT_ave, dT_sublayer, dT_top_int, dT_bot_int
1849 real,
dimension(nk+1) :: Til
1850 real,
dimension(nk+1) :: Tir
1851 real,
dimension(nk) :: aL_l
1852 real,
dimension(nk) :: aR_l
1853 real,
dimension(nk) :: aL_r
1854 real,
dimension(nk) :: aR_r
1857 real,
dimension(nk,2) :: Tid_l
1858 real,
dimension(nk,2) :: Tid_r
1859 real,
dimension(nk,deg+1) :: ppoly_r_coeffs_l
1860 real,
dimension(nk,deg+1) :: ppoly_r_coeffs_r
1861 real,
dimension(nk,deg+1) :: ppoly_r_S_l
1862 real,
dimension(nk,deg+1) :: ppoly_r_S_r
1863 logical :: down_flux
1865 if (continuous)
then
1871 ppoly_r_coeffs_l(:,:) = 0.
1872 ppoly_r_coeffs_r(:,:) = 0.
1876 call build_reconstructions_1d( remap_cs, nk, hl, tl, ppoly_r_coeffs_l, tid_l, &
1877 ppoly_r_s_l, imethod, h_neglect, h_neglect_edge )
1878 call build_reconstructions_1d( remap_cs, nk, hr, tr, ppoly_r_coeffs_r, tid_r, &
1879 ppoly_r_s_r, imethod, h_neglect, h_neglect_edge )
1882 do k_sublayer = 1, nsurf-1
1883 if (heff(k_sublayer) == 0.)
then
1884 flx(k_sublayer) = 0.
1886 if (continuous)
then
1887 klb = kol(k_sublayer+1)
1888 t_left_bottom = ( 1. - pil(k_sublayer+1) ) * til(klb) + pil(k_sublayer+1) * til(klb+1)
1889 klt = kol(k_sublayer)
1890 t_left_top = ( 1. - pil(k_sublayer) ) * til(klt) + pil(k_sublayer) * til(klt+1)
1891 t_left_layer =
ppm_ave(pil(k_sublayer), pil(k_sublayer+1) + real(klb-klt), &
1892 al_l(klt), ar_l(klt), tl(klt))
1894 krb = kor(k_sublayer+1)
1895 t_right_bottom = ( 1. - pir(k_sublayer+1) ) * tir(krb) + pir(k_sublayer+1) * tir(krb+1)
1896 krt = kor(k_sublayer)
1897 t_right_top = ( 1. - pir(k_sublayer) ) * tir(krt) + pir(k_sublayer) * tir(krt+1)
1898 t_right_layer =
ppm_ave(pir(k_sublayer), pir(k_sublayer+1) + real(krb-krt), &
1899 al_r(krt), ar_r(krt), tr(krt))
1900 dt_top = t_right_top - t_left_top
1901 dt_bottom = t_right_bottom - t_left_bottom
1902 dt_ave = 0.5 * ( dt_top + dt_bottom )
1903 dt_layer = t_right_layer - t_left_layer
1909 flx(k_sublayer) = dt_ave * heff(k_sublayer)
1913 ppoly_r_coeffs_l, t_left_top, t_left_bottom, t_left_sub, &
1914 t_left_top_int, t_left_bot_int, t_left_layer)
1916 ppoly_r_coeffs_r, t_right_top, t_right_bottom, t_right_sub, &
1917 t_right_top_int, t_right_bot_int, t_right_layer)
1919 dt_top = t_right_top - t_left_top
1920 dt_bottom = t_right_bottom - t_left_bottom
1921 dt_sublayer = t_right_sub - t_left_sub
1922 dt_top_int = t_right_top_int - t_left_top_int
1923 dt_bot_int = t_right_bot_int - t_left_bot_int
1927 down_flux = dt_top <= 0. .and. dt_bottom <= 0. .and. &
1928 dt_sublayer <= 0. .and. dt_top_int <= 0. .and. &
1930 down_flux = down_flux .or. &
1931 (dt_top >= 0. .and. dt_bottom >= 0. .and. &
1932 dt_sublayer >= 0. .and. dt_top_int >= 0. .and. &
1935 flx(k_sublayer) = dt_sublayer * heff(k_sublayer)
1937 flx(k_sublayer) = 0.
1946 subroutine neutral_surface_t_eval(nk, ns, k_sub, Ks, Ps, T_mean, T_int, deg, iMethod, T_poly, &
1947 T_top, T_bot, T_sub, T_top_int, T_bot_int, T_layer)
1948 integer,
intent(in ) :: nk
1949 integer,
intent(in ) :: ns
1950 integer,
intent(in ) :: k_sub
1951 integer,
dimension(ns),
intent(in ) :: Ks
1952 real,
dimension(ns),
intent(in ) :: Ps
1953 real,
dimension(nk),
intent(in ) :: T_mean
1954 real,
dimension(nk,2),
intent(in ) :: T_int
1955 integer,
intent(in ) :: deg
1956 integer,
intent(in ) :: iMethod
1957 real,
dimension(nk,deg+1),
intent(in ) :: T_poly
1958 real,
intent( out) :: T_top
1959 real,
intent( out) :: T_bot
1960 real,
intent( out) :: T_sub
1961 real,
intent( out) :: T_top_int
1962 real,
intent( out) :: T_bot_int
1963 real,
intent( out) :: T_layer
1965 integer :: kl, ks_top, ks_bot
1969 if ( ks(ks_top) /= ks(ks_bot) )
then
1970 call mom_error(fatal,
"Neutral surfaces span more than one layer")
1974 if ( (ps(ks_top) == 0.) .and. (ps(ks_bot) == 1.))
then
1979 if ( (kl > 1) .and. (ps(ks_top) == 0.) )
then
1980 t_top = t_int(kl-1,2)
1985 if ( (kl < nk) .and. (ps(ks_bot) == 1.) )
then
1986 t_bot = t_int(kl+1,1)
1991 t_sub =
average_value_ppoly(nk, t_mean, t_int, t_poly, imethod, kl, ps(ks_top), ps(ks_bot))
1994 t_layer = t_mean(kl)
2000 integer,
intent(in) :: nk
2001 real,
dimension(nk),
intent(in) :: Tl
2002 real,
dimension(nk+1),
intent(in) :: Ti
2003 real,
dimension(nk),
intent(inout) :: aL
2004 real,
dimension(nk),
intent(inout) :: aR
2011 if (
signum(1., ar(k) - tl(k))*
signum(1., tl(k) - al(k)) <= 0.0 )
then
2014 elseif ( sign(3., ar(k) - al(k)) * ( (tl(k) - al(k)) + (tl(k) - ar(k))) > abs(ar(k) - al(k)) )
then
2015 al(k) = tl(k) + 2.0 * ( tl(k) - ar(k) )
2016 elseif ( sign(3., ar(k) - al(k)) * ( (tl(k) - al(k)) + (tl(k) - ar(k))) < -abs(ar(k) - al(k)) )
then
2017 ar(k) = tl(k) + 2.0 * ( tl(k) - al(k) )
2024 logical,
intent(in) :: verbose
2034 logical,
intent(in) :: verbose
2036 integer,
parameter :: nk = 4
2037 real,
dimension(nk+1) :: til, tir1, tir2, tir4, tio
2038 real,
dimension(nk) :: tl
2039 real,
dimension(nk+1) :: sil
2040 real,
dimension(nk+1) :: pil, pir4
2041 real,
dimension(2*nk+2) :: pilro, pirlo
2042 integer,
dimension(2*nk+2) :: kol, kor
2043 real,
dimension(2*nk+1) :: heff
2044 real,
dimension(2*nk+1) :: flx
2047 real :: h_neglect, h_neglect_edge
2049 h_neglect_edge = 1.0e-10 ; h_neglect = 1.0e-30
2054 write(*,*)
'==== MOM_neutral_diffusion: ndiff_unit_tests_continuous ='
2057 test_fv_diff(v,1.,1.,1., 0.,1.,2., 1.,
'FV: Straight line on uniform grid')
2059 test_fv_diff(v,1.,1.,0., 0.,4.,8., 7.,
'FV: Vanished right cell')
2061 test_fv_diff(v,0.,1.,1., 0.,4.,8., 7.,
'FV: Vanished left cell')
2063 test_fv_diff(v,1.,2.,4., 0.,3.,9., 4.,
'FV: Stretched grid')
2065 test_fv_diff(v,2.,0.,2., 0.,1.,2., 0.,
'FV: Vanished middle cell')
2067 test_fv_diff(v,0.,1.,0., 0.,1.,2., 2.,
'FV: Vanished on both sides')
2069 test_fv_diff(v,1.,0.,0., 0.,1.,2., 0.,
'FV: Two vanished cell sides')
2071 test_fv_diff(v,0.,0.,0., 0.,1.,2., 0.,
'FV: All vanished cells')
2074 test_fvlsq_slope(v,1.,1.,1., 0.,1.,2., 1.,
'LSQ: Straight line on uniform grid')
2090 call interface_scalar(4, (/10.,10.,10.,10./), (/24.,18.,12.,6./), tio, 1, h_neglect)
2094 test_data1d(v,5, tio, (/24.,22.5,15.,7.5,6./),
'Linear profile, linear interface temperatures')
2095 call interface_scalar(4, (/10.,10.,10.,10./), (/24.,18.,12.,6./), tio, 2, h_neglect)
2097 test_data1d(v,5, tio, (/24.,22.,15.,8.,6./),
'Linear profile, PPM interface temperatures')
2100 test_ifndp(v,-1.0, 0., 1.0, 1.0, 0.5,
'Check mid-point')
2102 test_ifndp(v, 0.0, 0., 1.0, 1.0, 0.0,
'Check bottom')
2104 test_ifndp(v, 0.1, 0., 1.1, 1.0, 0.0,
'Check below')
2106 test_ifndp(v,-1.0, 0., 0.0, 1.0, 1.0,
'Check top')
2108 test_ifndp(v,-1.0, 0., -0.1, 1.0, 1.0,
'Check above')
2110 test_ifndp(v,-1.0, 0., 3.0, 1.0, 0.25,
'Check 1/4')
2112 test_ifndp(v,-3.0, 0., 1.0, 1.0, 0.75,
'Check 3/4')
2114 test_ifndp(v, 1.0, 0., 1.0, 1.0, 0.0,
'Check dRho=0 below')
2116 test_ifndp(v,-1.0, 0., -1.0, 1.0, 1.0,
'Check dRho=0 above')
2118 test_ifndp(v, 0.0, 0., 0.0, 1.0, 0.5,
'Check dRho=0 mid')
2120 test_ifndp(v,-2.0, .5, 5.0, 0.5, 0.5,
'Check dP=0')
2124 (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), &
2125 (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &
2126 (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), &
2127 (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &
2128 pilro, pirlo, kol, kor, heff)
2130 (/1,1,2,2,3,3,3,3/), &
2131 (/1,1,2,2,3,3,3,3/), &
2132 (/0.,0.,0.,0.,0.,0.,1.,1./), &
2133 (/0.,0.,0.,0.,0.,0.,1.,1./), &
2134 (/0.,10.,0.,10.,0.,10.,0./), &
2135 'Identical columns')
2138 (/0.,0.,10.,10.,20.,20.,30.,30./),
'... left positions')
2141 (/0.,0.,10.,10.,20.,20.,30.,30./),
'... right positions')
2143 (/20.,16.,12./), (/20.,16.,12./), &
2144 pilro, pirlo, kol, kor, heff, flx, .true., &
2145 h_neglect, h_neglect_edge=h_neglect_edge)
2147 (/0.,0.,0.,0.,0.,0.,0./),
'Identical columns, rho flux (=0)')
2149 (/-1.,-1.,-1./), (/1.,1.,1./), &
2150 pilro, pirlo, kol, kor, heff, flx, .true., &
2151 h_neglect, h_neglect_edge=h_neglect_edge)
2153 (/0.,20.,0.,20.,0.,20.,0./),
'Identical columns, S flux')
2157 (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), &
2158 (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &
2159 (/0.,10.,20.,30./), (/20.,16.,12.,8./), (/0.,0.,0.,0./), &
2160 (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &
2161 pilro, pirlo, kol, kor, heff)
2163 (/1,1,2,2,3,3,3,3/), &
2164 (/1,1,1,2,2,3,3,3/), &
2165 (/0.,0.5,0.,0.5,0.,0.5,1.,1./), &
2166 (/0.,0.,0.5,0.,0.5,0.,0.5,1./), &
2167 (/0.,5.,5.,5.,5.,5.,0./), &
2168 'Right column slightly cooler')
2171 (/0.,5.,10.,15.,20.,25.,30.,30./),
'... left positions')
2174 (/0.,0.,5.,10.,15.,20.,25.,30./),
'... right positions')
2178 (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), &
2179 (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &
2180 (/0.,10.,20.,30./), (/24.,20.,16.,12./), (/0.,0.,0.,0./), &
2181 (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &
2182 pilro, pirlo, kol, kor, heff)
2184 (/1,1,1,2,2,3,3,3/), &
2185 (/1,1,2,2,3,3,3,3/), &
2186 (/0.,0.,0.5,0.,0.5,0.,0.5,1./), &
2187 (/0.,0.5,0.,0.5,0.,0.5,1.,1./), &
2188 (/0.,5.,5.,5.,5.,5.,0./), &
2189 'Right column slightly warmer')
2193 (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), &
2194 (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &
2195 (/0.,10.,20.,30./), (/16.,12.,8.,4./), (/0.,0.,0.,0./), &
2196 (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &
2197 pilro, pirlo, kol, kor, heff)
2199 (/1,2,2,3,3,3,3,3/), &
2200 (/1,1,1,1,2,2,3,3/), &
2201 (/0.,0.,0.5,0.,0.5,1.,1.,1./), &
2202 (/0.,0.,0.,0.5,0.,0.5,0.,1./), &
2203 (/0.,0.,5.,5.,5.,0.,0./), &
2204 'Right column somewhat cooler')
2208 (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), &
2209 (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &
2210 (/0.,10.,20.,30./), (/9.,7.,5.,3./), (/0.,0.,0.,0./), &
2211 (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &
2212 pilro, pirlo, kol, kor, heff)
2214 (/1,2,3,3,3,3,3,3/), &
2215 (/1,1,1,1,1,2,3,3/), &
2216 (/0.,0.,0.,1.,1.,1.,1.,1./), &
2217 (/0.,0.,0.,0.,0.,0.,0.,1./), &
2218 (/0.,0.,0.,0.,0.,0.,0./), &
2219 'Right column much cooler')
2223 (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), &
2224 (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &
2225 (/0.,10.,20.,30./), (/14.,14.,10.,2./), (/0.,0.,0.,0./), &
2226 (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &
2227 pilro, pirlo, kol, kor, heff)
2229 (/1,2,3,3,3,3,3,3/), &
2230 (/1,1,1,1,2,3,3,3/), &
2231 (/0.,0.,0.,0.,0.,1.,1.,1./), &
2232 (/0.,0.,0.,0.,0.,0.,0.,1./), &
2233 (/0.,0.,0.,0.,10.,0.,0./), &
2234 'Right column with mixed layer')
2238 (/0.,10.,20.,30./), (/14.,14.,10.,2./), (/0.,0.,0.,0./), &
2239 (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &
2240 (/0.,10.,20.,30./), (/14.,14.,10.,2./), (/0.,0.,0.,0./), &
2241 (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &
2242 pilro, pirlo, kol, kor, heff)
2244 (/1,1,2,2,3,3,3,3/), &
2245 (/1,1,2,2,3,3,3,3/), &
2246 (/0.,0.,0.,0.,0.,0.,1.,1./), &
2247 (/0.,0.,0.,0.,0.,0.,1.,1./), &
2248 (/0.,10.,0.,10.,0.,10.,0./), &
2249 'Identical columns with mixed layer')
2253 (/0.,10.,20.,30./), (/14.,14.,10.,2./), (/0.,0.,0.,0./), &
2254 (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &
2255 (/0.,10.,20.,30./), (/10.,14.,12.,4./), (/0.,0.,0.,0./), &
2256 (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &
2257 pilro, pirlo, kol, kor, heff)
2259 (/1,2,3,3,3,3,3,3/), &
2260 (/1,1,1,2,3,3,3,3/), &
2261 (/0.,0.,0.,0.,0.,0.,.75,1./), &
2262 (/0.,0.,0.,0.,0.,0.25,1.,1./), &
2263 (/0.,0.,0.,0.,0.,7.5,0./), &
2264 'Right column with unstable mixed layer')
2268 (/0.,10.,20.,30./), (/10.,14.,12.,4./), (/0.,0.,0.,0./), &
2269 (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &
2270 (/0.,10.,20.,30./), (/14.,14.,10.,2./), (/0.,0.,0.,0./), &
2271 (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &
2272 pilro, pirlo, kol, kor, heff)
2274 (/1,1,1,2,3,3,3,3/), &
2275 (/1,2,3,3,3,3,3,3/), &
2276 (/0.,0.,0.,0.,0.,0.25,1.,1./), &
2277 (/0.,0.,0.,0.,0.,0.,.75,1./), &
2278 (/0.,0.,0.,0.,0.,7.5,0./), &
2279 'Left column with unstable mixed layer')
2283 (/0.,10.,20.,30./), (/8.,12.,10.,2./), (/0.,0.,0.,0./), &
2284 (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &
2285 (/0.,10.,20.,30./), (/10.,14.,12.,4./), (/0.,0.,0.,0./), &
2286 (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &
2287 pilro, pirlo, kol, kor, heff)
2289 (/1,1,1,1,2,3,3,3/), &
2290 (/1,2,3,3,3,3,3,3/), &
2291 (/0.,0.,0.,0.,0.,0.,0.75,1./), &
2292 (/0.,0.,0.,0.5,0.5,0.5,1.,1./), &
2293 (/0.,0.,0.,0.,0.,6.,0./), &
2294 'Two unstable mixed layers')
2301 logical,
intent(in) :: verbose
2303 integer,
parameter :: nk = 3
2304 integer,
parameter :: ns = nk*4
2305 real,
dimension(nk) :: sl, sr, tl, tr, hl, hr
2306 real,
dimension(nk,2) :: til, sil, tir, sir
2307 real,
dimension(nk,2) :: pres_l, pres_r
2308 integer,
dimension(ns) :: kol, kor
2309 real,
dimension(ns) :: pol, por
2310 real,
dimension(ns-1) :: heff, flx
2314 real,
dimension(nk,2) :: ppoly_t_l, ppoly_t_r
2315 real,
dimension(nk,2) :: ppoly_s_l, ppoly_s_r
2316 real,
dimension(nk,2) :: drdt, drds
2317 logical,
dimension(nk) :: stable_l, stable_r
2319 integer :: ns_l, ns_r
2320 real :: h_neglect, h_neglect_edge
2326 write(*,*)
'==== MOM_neutral_diffusion: ndiff_unit_tests_discontinuous ='
2328 h_neglect = 1.0e-30 ; h_neglect_edge = 1.0e-10
2333 call eos_manual_init(cs%EOS, form_of_eos =
eos_linear, drho_dt = -1., drho_ds = 0.)
2334 sl(:) = 0. ; sr(:) = 0. ; ; sil(:,:) = 0. ; sir(:,:) = 0.
2335 ppoly_t_l(:,:) = 0.; ppoly_t_r(:,:) = 0.
2336 ppoly_s_l(:,:) = 0.; ppoly_s_r(:,:) = 0.
2340 hl = (/10.,10.,10./) ; hr = (/10.,10.,10./)
2341 pres_l(1,1) = 0. ; pres_l(1,2) = hl(1) ; pres_r(1,1) = 0. ; pres_r(1,2) = hr(1)
2343 pres_l(k,1) = pres_l(k-1,2)
2344 pres_l(k,2) = pres_l(k,1) + hl(k)
2345 pres_r(k,1) = pres_r(k-1,2)
2346 pres_r(k,2) = pres_r(k,1) + hr(k)
2348 cs%delta_rho_form =
'mid_pressure'
2349 cs%neutral_pos_method = 1
2351 til(1,:) = (/ 22.00, 18.00 /); til(2,:) = (/ 18.00, 14.00 /); til(3,:) = (/ 14.00, 10.00 /);
2352 tir(1,:) = (/ 22.00, 18.00 /); tir(2,:) = (/ 18.00, 14.00 /); tir(3,:) = (/ 14.00, 10.00 /);
2356 pres_r, hr, tir, sir, ppoly_t_r, ppoly_s_r, stable_r, pol, por, kol, kor, heff)
2358 (/ 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3 /), &
2359 (/ 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3 /), &
2360 (/ 0.00, 0.00, 1.00, 1.00, 0.00, 0.00, 1.00, 1.00, 0.00, 0.00, 1.00, 1.00 /), &
2361 (/ 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 1.00, 1.00 /), &
2362 (/ 0.00, 10.00, 0.00, 0.00, 0.00, 10.00, 0.00, 0.00, 0.00, 10.00, 0.00 /), &
2363 'Identical Columns')
2365 til(1,:) = (/ 22.00, 18.00 /); til(2,:) = (/ 18.00, 14.00 /); til(3,:) = (/ 14.00, 10.00 /);
2366 tir(1,:) = (/ 20.00, 16.00 /); tir(2,:) = (/ 16.00, 12.00 /); tir(3,:) = (/ 12.00, 8.00 /);
2370 pres_r, hr, tir, sir, ppoly_t_r, ppoly_s_r, stable_r, pol, por, kol, kor, heff)
2372 (/ 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3 /), &
2373 (/ 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3 /), &
2374 (/ 0.00, 0.50, 1.00, 0.00, 0.50, 0.50, 1.00, 0.00, 0.50, 0.50, 1.00, 1.00 /), &
2375 (/ 0.00, 0.00, 0.50, 0.50, 1.00, 0.00, 0.50, 0.50, 1.00, 0.00, 0.50, 1.00 /), &
2376 (/ 0.00, 5.00, 0.00, 5.00, 0.00, 5.00, 0.00, 5.00, 0.00, 5.00, 0.00 /), &
2377 'Right slightly cooler')
2379 til(1,:) = (/ 20.00, 16.00 /); til(2,:) = (/ 16.00, 12.00 /); til(3,:) = (/ 12.00, 8.00 /);
2380 tir(1,:) = (/ 22.00, 18.00 /); tir(2,:) = (/ 18.00, 14.00 /); tir(3,:) = (/ 14.00, 10.00 /);
2384 pres_r, hr, tir, sir, ppoly_t_r, ppoly_s_r, stable_r, pol, por, kol, kor, heff)
2386 (/ 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3 /), &
2387 (/ 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3 /), &
2388 (/ 0.00, 0.00, 0.50, 0.50, 1.00, 0.00, 0.50, 0.50, 1.00, 0.00, 0.50, 1.00 /), &
2389 (/ 0.00, 0.50, 1.00, 0.00, 0.50, 0.50, 1.00, 0.00, 0.50, 0.50, 1.00, 1.00 /), &
2390 (/ 0.00, 5.00, 0.00, 5.00, 0.00, 5.00, 0.00, 5.00, 0.00, 5.00, 0.00 /), &
2391 'Left slightly cooler')
2393 til(1,:) = (/ 22.00, 20.00 /); til(2,:) = (/ 18.00, 16.00 /); til(3,:) = (/ 14.00, 12.00 /);
2394 tir(1,:) = (/ 32.00, 24.00 /); tir(2,:) = (/ 22.00, 14.00 /); tir(3,:) = (/ 12.00, 4.00 /);
2398 pres_r, hr, tir, sir, ppoly_t_r, ppoly_s_r, stable_r, pol, por, kol, kor, heff)
2400 (/ 1, 1, 1, 1, 1, 2, 2, 3, 3, 3, 3, 3 /), &
2401 (/ 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3 /), &
2402 (/ 0.00, 0.00, 0.00, 0.00, 1.00, 0.00, 1.00, 0.00, 0.00, 1.00, 1.00, 1.00 /), &
2403 (/ 0.00, 1.00, 0.00, 0.00, 0.25, 0.50, 0.75, 1.00, 0.00, 0.00, 0.00, 1.00 /), &
2404 (/ 0.00, 0.00, 0.00, 4.00, 0.00, 4.00, 0.00, 0.00, 0.00, 0.00, 0.00 /), &
2405 'Right more strongly stratified')
2407 til(1,:) = (/ 22.00, 18.00 /); til(2,:) = (/ 18.00, 14.00 /); til(3,:) = (/ 14.00, 10.00 /);
2408 tir(1,:) = (/ 14.00, 14.00 /); tir(2,:) = (/ 14.00, 14.00 /); tir(3,:) = (/ 12.00, 8.00 /);
2412 pres_r, hr, tir, sir, ppoly_t_r, ppoly_s_r, stable_r, pol, por, kol, kor, heff)
2414 (/ 1, 1, 1, 1, 1, 1, 2, 2, 3, 3, 3, 3 /), &
2415 (/ 1, 1, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3 /), &
2416 (/ 0.00, 0.00, 0.00, 0.00, 0.00, 1.00, 0.00, 1.00, 0.00, 0.50, 1.00, 1.00 /), &
2417 (/ 0.00, 1.00, 0.00, 1.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.50, 1.00 /), &
2418 (/ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 5.00, 0.00 /), &
2419 'Deep Mixed layer on the right')
2421 til(1,:) = (/ 14.00, 14.00 /); til(2,:) = (/ 14.00, 12.00 /); til(3,:) = (/ 10.00, 8.00 /);
2422 tir(1,:) = (/ 14.00, 14.00 /); tir(2,:) = (/ 14.00, 14.00 /); tir(3,:) = (/ 14.00, 14.00 /);
2426 pres_r, hr, tir, sir, ppoly_t_r, ppoly_s_r, stable_r, pol, por, kol, kor, heff)
2428 (/ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3 /), &
2429 (/ 1, 1, 1, 1, 2, 2, 3, 3, 3, 3, 3, 3 /), &
2430 (/ 0.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00 /), &
2431 (/ 0.00, 0.00, 0.00, 1.00, 0.00, 1.00, 0.00, 1.00, 1.00, 1.00, 1.00, 1.00 /), &
2432 (/ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 /), &
2433 'Right unstratified column')
2435 til(1,:) = (/ 14.00, 14.00 /); til(2,:) = (/ 14.00, 12.00 /); til(3,:) = (/ 10.00, 8.00 /);
2436 tir(1,:) = (/ 14.00, 14.00 /); tir(2,:) = (/ 14.00, 14.00 /); tir(3,:) = (/ 12.00, 4.00 /);
2440 pres_r, hr, tir, sir, ppoly_t_r, ppoly_s_r, stable_r, pol, por, kol, kor, heff)
2442 (/ 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3 /), &
2443 (/ 1, 1, 1, 1, 2, 2, 3, 3, 3, 3, 3, 3 /), &
2444 (/ 0.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.00, 1.00, 1.00, 0.00, 1.00, 1.00 /), &
2445 (/ 0.00, 0.00, 0.00, 1.00, 0.00, 1.00, 0.00, 0.00, 0.00, 0.25, 0.50, 1.00 /), &
2446 (/ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 4.00, 0.00 /), &
2447 'Right unstratified column')
2449 til(1,:) = (/ 14.00, 14.00 /); til(2,:) = (/ 14.00, 10.00 /); til(3,:) = (/ 10.00, 2.00 /);
2450 tir(1,:) = (/ 14.00, 14.00 /); tir(2,:) = (/ 14.00, 10.00 /); tir(3,:) = (/ 10.00, 2.00 /);
2454 pres_r, hr, tir, sir, ppoly_t_r, ppoly_s_r, stable_r, pol, por, kol, kor, heff)
2456 (/ 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3 /), &
2457 (/ 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3 /), &
2458 (/ 0.00, 1.00, 1.00, 1.00, 0.00, 0.00, 1.00, 1.00, 0.00, 0.00, 1.00, 1.00 /), &
2459 (/ 0.00, 0.00, 0.00, 1.00, 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 1.00, 1.00 /), &
2460 (/ 0.00, 0.00, 0.00, 0.00, 0.00, 10.00, 0.00, 0.00, 0.00, 10.00, 0.00 /), &
2461 'Identical columns with mixed layer')
2463 til(1,:) = (/ 14.00, 12.00 /); til(2,:) = (/ 10.00, 10.00 /); til(3,:) = (/ 8.00, 2.00 /);
2464 tir(1,:) = (/ 14.00, 12.00 /); tir(2,:) = (/ 12.00, 8.00 /); tir(3,:) = (/ 8.00, 2.00 /);
2468 pres_r, hr, tir, sir, ppoly_t_r, ppoly_s_r, stable_r, pol, por, kol, kor, heff)
2470 (/ 1, 1, 1, 1, 2, 2, 3, 3, 3, 3, 3, 3 /), &
2471 (/ 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3 /), &
2472 (/ 0.00, 0.00, 1.00, 1.00, 0.00, 1.00, 0.00, 0.00, 0.00, 0.00, 1.00, 1.00 /), &
2473 (/ 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 0.00, 1.00, 1.00, 0.00, 1.00, 1.00 /), &
2474 (/ 0.00, 10.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 10.00, 0.00 /), &
2475 'Left interior unstratified')
2477 til(1,:) = (/ 12.00, 12.00 /); til(2,:) = (/ 12.00, 10.00 /); til(3,:) = (/ 10.00, 6.00 /);
2478 tir(1,:) = (/ 12.00, 10.00 /); tir(2,:) = (/ 10.00, 12.00 /); tir(3,:) = (/ 8.00, 4.00 /);
2482 pres_r, hr, tir, sir, ppoly_t_r, ppoly_s_r, stable_r, pol, por, kol, kor, heff)
2484 (/ 1, 1, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3 /), &
2485 (/ 1, 1, 1, 1, 1, 1, 2, 2, 3, 3, 3, 3 /), &
2486 (/ 0.00, 1.00, 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 0.00, 0.50, 1.00, 1.00 /), &
2487 (/ 0.00, 0.00, 0.00, 0.00, 1.00, 1.00, 0.00, 1.00, 0.00, 0.00, 0.50, 1.00 /), &
2488 (/ 0.00, 0.00, 0.00, 10.00, 0.00, 0.00, 0.00, 0.00, 0.00, 5.00, 0.00 /), &
2489 'Left mixed layer, Right unstable interior')
2491 til(1,:) = (/ 14.00, 14.00 /); til(2,:) = (/ 10.00, 10.00 /); til(3,:) = (/ 8.00, 6.00 /);
2492 tir(1,:) = (/ 10.00, 14.00 /); tir(2,:) = (/ 16.00, 16.00 /); tir(3,:) = (/ 12.00, 4.00 /);
2496 pres_r, hr, tir, sir, ppoly_t_r, ppoly_s_r, stable_r, pol, por, kol, kor, heff)
2498 (/ 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3 /), &
2499 (/ 1, 1, 1, 1, 1, 1, 2, 2, 3, 3, 3, 3 /), &
2500 (/ 0.00, 1.00, 0.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.00, 0.00, 1.00, 1.00 /), &
2501 (/ 0.00, 0.00, 0.00, 0.00, 0.00, 1.00, 0.00, 1.00, 0.00, 0.50, 0.75, 1.00 /), &
2502 (/ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 4.00, 0.00 /), &
2503 'Left thick mixed layer, Right unstable mixed')
2505 til(1,:) = (/ 8.00, 12.00 /); til(2,:) = (/ 12.00, 10.00 /); til(3,:) = (/ 8.00, 4.00 /);
2506 tir(1,:) = (/ 10.00, 14.00 /); tir(2,:) = (/ 14.00, 12.00 /); tir(3,:) = (/ 10.00, 6.00 /);
2510 pres_r, hr, tir, sir, ppoly_t_r, ppoly_s_r, stable_r, pol, por, kol, kor, heff)
2512 (/ 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3 /), &
2513 (/ 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3 /), &
2514 (/ 0.00, 1.00, 1.00, 1.00, 0.00, 0.00, 0.00, 1.00, 0.00, 0.00, 0.50, 1.00 /), &
2515 (/ 0.00, 0.00, 0.00, 1.00, 0.00, 1.00, 1.00, 0.00, 0.00, 0.50, 1.00, 1.00 /), &
2516 (/ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 5.00, 0.00 /), &
2517 'Unstable mixed layers, left cooler')
2519 call eos_manual_init(cs%EOS, form_of_eos =
eos_linear, drho_dt = -1., drho_ds = 2.)
2528 0., -0.2, 0., 10., -0.2, 0., &
2529 (/12.,-4./), (/34.,0./)),
"Temp Uniform Linearized Alpha/Beta"))
2533 0., 0., 0.8, 10., 0., 0.8, &
2534 (/12.,0./), (/34.,2./)),
"Salt Uniform Linearized Alpha/Beta"))
2538 0., -0.5, 0.5, 10., -0.5, 0.5, &
2539 (/12.,-4./), (/34.,2./)),
"Temp/salt Uniform Linearized Alpha/Beta"))
2543 0., -0.4, 0., 10., -0.6, 0., &
2544 (/12.,-4./), (/34.,0./)),
"Temp stratified Linearized Alpha/Beta"))
2548 0., 0., 1.0, 10., 0., 0.5, &
2549 (/12.,0./), (/34.,2./)),
"Salt stratified Linearized Alpha/Beta"))
2555 logical function test_fv_diff(verbose, hkm1, hk, hkp1, Skm1, Sk, Skp1, Ptrue, title)
2556 logical,
intent(in) :: verbose
2557 real,
intent(in) :: hkm1
2558 real,
intent(in) :: hk
2559 real,
intent(in) :: hkp1
2560 real,
intent(in) :: skm1
2561 real,
intent(in) :: sk
2562 real,
intent(in) :: skp1
2563 real,
intent(in) :: ptrue
2564 character(len=*),
intent(in) :: title
2570 pret =
fv_diff(hkm1, hk, hkp1, skm1, sk, skp1)
2576 write(stdunit,
'(a)') title
2578 write(stdunit,
'(2(x,a,f20.16),x,a)')
'pRet=',pret,
'pTrue=',ptrue,
'WRONG!'
2580 write(stdunit,
'(2(x,a,f20.16))')
'pRet=',pret,
'pTrue=',ptrue
2587 logical function test_fvlsq_slope(verbose, hkm1, hk, hkp1, Skm1, Sk, Skp1, Ptrue, title)
2588 logical,
intent(in) :: verbose
2589 real,
intent(in) :: hkm1
2590 real,
intent(in) :: hk
2591 real,
intent(in) :: hkp1
2592 real,
intent(in) :: skm1
2593 real,
intent(in) :: sk
2594 real,
intent(in) :: skp1
2595 real,
intent(in) :: ptrue
2596 character(len=*),
intent(in) :: title
2602 pret =
fvlsq_slope(hkm1, hk, hkp1, skm1, sk, skp1)
2608 write(stdunit,
'(a)') title
2610 write(stdunit,
'(2(x,a,f20.16),x,a)')
'pRet=',pret,
'pTrue=',ptrue,
'WRONG!'
2612 write(stdunit,
'(2(x,a,f20.16))')
'pRet=',pret,
'pTrue=',ptrue
2619 logical function test_ifndp(verbose, rhoNeg, Pneg, rhoPos, Ppos, Ptrue, title)
2620 logical,
intent(in) :: verbose
2621 real,
intent(in) :: rhoneg
2622 real,
intent(in) :: pneg
2623 real,
intent(in) :: rhopos
2624 real,
intent(in) :: ppos
2625 real,
intent(in) :: ptrue
2626 character(len=*),
intent(in) :: title
2638 write(stdunit,
'(a)') title
2640 write(stdunit,
'(4(x,a,f20.16),2(x,a,1pe22.15),x,a)') &
2641 'r1=',rhoneg,
'p1=',pneg,
'r2=',rhopos,
'p2=',ppos,
'pRet=',pret,
'pTrue=',ptrue,
'WRONG!'
2643 write(stdunit,
'(4(x,a,f20.16),2(x,a,1pe22.15))') &
2644 'r1=',rhoneg,
'p1=',pneg,
'r2=',rhopos,
'p2=',ppos,
'pRet=',pret,
'pTrue=',ptrue
2651 logical function test_data1d(verbose, nk, Po, Ptrue, title)
2652 logical,
intent(in) :: verbose
2653 integer,
intent(in) :: nk
2654 real,
dimension(nk),
intent(in) :: po
2655 real,
dimension(nk),
intent(in) :: ptrue
2656 character(len=*),
intent(in) :: title
2659 integer :: k, stdunit
2669 write(stdunit,
'(a)') title
2671 if (po(k) /= ptrue(k))
then
2673 write(stdunit,
'(a,i2,2(x,a,f20.16),x,a,1pe22.15,x,a)') &
2674 'k=',k,
'Po=',po(k),
'Ptrue=',ptrue(k),
'err=',po(k)-ptrue(k),
'WRONG!'
2677 write(stdunit,
'(a,i2,2(x,a,f20.16),x,a,1pe22.15)') &
2678 'k=',k,
'Po=',po(k),
'Ptrue=',ptrue(k),
'err=',po(k)-ptrue(k)
2686 logical function test_data1di(verbose, nk, Po, Ptrue, title)
2687 logical,
intent(in) :: verbose
2688 integer,
intent(in) :: nk
2689 integer,
dimension(nk),
intent(in) :: po
2690 integer,
dimension(nk),
intent(in) :: ptrue
2691 character(len=*),
intent(in) :: title
2694 integer :: k, stdunit
2704 write(stdunit,
'(a)') title
2706 if (po(k) /= ptrue(k))
then
2708 write(stdunit,
'(a,i2,2(x,a,i5),x,a)')
'k=',k,
'Io=',po(k),
'Itrue=',ptrue(k),
'WRONG!'
2711 write(stdunit,
'(a,i2,2(x,a,i5))')
'k=',k,
'Io=',po(k),
'Itrue=',ptrue(k)
2720 logical function test_nsp(verbose, ns, KoL, KoR, pL, pR, hEff, KoL0, KoR0, pL0, pR0, hEff0, title)
2721 logical,
intent(in) :: verbose
2722 integer,
intent(in) :: ns
2723 integer,
dimension(ns),
intent(in) :: kol
2724 integer,
dimension(ns),
intent(in) :: kor
2725 real,
dimension(ns),
intent(in) :: pl
2726 real,
dimension(ns),
intent(in) :: pr
2727 real,
dimension(ns-1),
intent(in) :: heff
2728 integer,
dimension(ns),
intent(in) :: kol0
2729 integer,
dimension(ns),
intent(in) :: kor0
2730 real,
dimension(ns),
intent(in) :: pl0
2731 real,
dimension(ns),
intent(in) :: pr0
2732 real,
dimension(ns-1),
intent(in) :: heff0
2733 character(len=*),
intent(in) :: title
2736 integer :: k, stdunit
2737 logical :: this_row_failed
2743 if (heff(k) /= heff0(k))
test_nsp = .true.
2750 write(stdunit,
'(a)') title
2752 this_row_failed =
compare_nsp_row(kol(k), kor(k), pl(k), pr(k), kol0(k), kor0(k), pl0(k), pr0(k))
2753 if (this_row_failed)
then
2754 write(stdunit,10) k,kol(k),pl(k),kor(k),pr(k),
' <-- WRONG!'
2755 write(stdunit,10) k,kol0(k),pl0(k),kor0(k),pr0(k),
' <-- should be this'
2757 write(stdunit,10) k,kol(k),pl(k),kor(k),pr(k)
2760 if (heff(k) /= heff0(k))
then
2761 write(stdunit,
'(i3,8x,"layer hEff =",2(f20.16,a))') k,heff(k),
" .neq. ",heff0(k),
' <-- WRONG!'
2763 write(stdunit,
'(i3,8x,"layer hEff =",f20.16)') k,heff(k)
2770 10
format(
"ks=",i3,
" kL=",i3,
" pL=",f20.16,
" kR=",i3,
" pR=",f20.16,a)
2774 logical function compare_nsp_row(KoL, KoR, pL, pR, KoL0, KoR0, pL0, pR0)
2775 integer,
intent(in) :: kol
2776 integer,
intent(in) :: kor
2777 real,
intent(in) :: pl
2778 real,
intent(in) :: pr
2779 integer,
intent(in) :: kol0
2780 integer,
intent(in) :: kor0
2781 real,
intent(in) :: pl0
2782 real,
intent(in) :: pr0
2792 logical function test_rnp(expected_pos, test_pos, title)
2793 real,
intent(in) :: expected_pos
2794 real,
intent(in) :: test_pos
2795 character(len=*),
intent(in) :: title
2797 integer :: stdunit = 6
2798 test_rnp = abs(expected_pos - test_pos) > 2*epsilon(test_pos)
2800 write(stdunit,
'(A, f20.16, " .neq. ", f20.16, " <-- WRONG")') title, expected_pos, test_pos
2802 write(stdunit,
'(A, f20.16, " == ", f20.16)') title, expected_pos, test_pos
2809 if (
associated(cs))
deallocate(cs)