6 use mom_cpu_clock,
only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, clock_routine
17 implicit none ;
private
19 #include <MOM_memory.h>
47 real :: cfl_limit_adjust
48 logical :: aggress_adjust
54 logical :: better_iter
57 logical :: use_visc_rem_max
59 logical :: marginal_faces
68 integer :: ish, ieh, jsh, jeh
76 subroutine continuity_ppm(u, v, hin, h, uh, vh, dt, G, GV, US, CS, uhbt, vhbt, OBC, &
77 visc_rem_u, visc_rem_v, u_cor, v_cor, BT_cont)
80 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
82 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
84 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
86 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
88 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
90 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
92 real,
intent(in) :: dt
95 real,
dimension(SZIB_(G),SZJ_(G)), &
96 optional,
intent(in) :: uhbt
98 real,
dimension(SZI_(G),SZJB_(G)), &
99 optional,
intent(in) :: vhbt
102 optional,
pointer :: obc
103 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
104 optional,
intent(in) :: visc_rem_u
110 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
111 optional,
intent(in) :: visc_rem_v
117 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
118 optional,
intent(out) :: u_cor
120 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
121 optional,
intent(out) :: v_cor
130 integer :: is, ie, js, je, nz, stencil
134 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
136 h_min = gv%Angstrom_H
138 if (.not.
associated(cs))
call mom_error(fatal, &
139 "MOM_continuity_PPM: Module must be initialized before it is used.")
140 x_first = (mod(g%first_direction,2) == 0)
142 if (
present(visc_rem_u) .neqv.
present(visc_rem_v))
call mom_error(fatal, &
143 "MOM_continuity_PPM: Either both visc_rem_u and visc_rem_v or neither"// &
144 " one must be present in call to continuity_PPM.")
146 stencil = 3 ;
if (cs%simple_2nd) stencil = 2 ;
if (cs%upwind_1st) stencil = 1
150 lb%ish = g%isc ; lb%ieh = g%iec
151 lb%jsh = g%jsc-stencil ; lb%jeh = g%jec+stencil
152 call zonal_mass_flux(u, hin, uh, dt, g, gv, us, cs, lb, uhbt, obc, visc_rem_u, u_cor, bt_cont)
156 do k=1,nz ;
do j=lb%jsh,lb%jeh ;
do i=lb%ish,lb%ieh
157 h(i,j,k) = hin(i,j,k) - dt * g%IareaT(i,j) * (uh(i,j,k) - uh(i-1,j,k))
160 enddo ;
enddo ;
enddo
163 lb%ish = g%isc ; lb%ieh = g%iec ; lb%jsh = g%jsc ; lb%jeh = g%jec
167 call meridional_mass_flux(v, h, vh, dt, g, gv, us, cs, lb, vhbt, obc, visc_rem_v, v_cor, bt_cont)
171 do k=1,nz ;
do j=lb%jsh,lb%jeh ;
do i=lb%ish,lb%ieh
172 h(i,j,k) = h(i,j,k) - dt * g%IareaT(i,j) * (vh(i,j,k) - vh(i,j-1,k))
174 if (h(i,j,k) < h_min) h(i,j,k) = h_min
175 enddo ;
enddo ;
enddo
180 lb%ish = g%isc-stencil ; lb%ieh = g%iec+stencil
181 lb%jsh = g%jsc ; lb%jeh = g%jec
183 call meridional_mass_flux(v, hin, vh, dt, g, gv, us, cs, lb, vhbt, obc, visc_rem_v, v_cor, bt_cont)
187 do k=1,nz ;
do j=lb%jsh,lb%jeh ;
do i=lb%ish,lb%ieh
188 h(i,j,k) = hin(i,j,k) - dt * g%IareaT(i,j) * (vh(i,j,k) - vh(i,j-1,k))
189 enddo ;
enddo ;
enddo
194 lb%ish = g%isc ; lb%ieh = g%iec ; lb%jsh = g%jsc ; lb%jeh = g%jec
195 call zonal_mass_flux(u, h, uh, dt, g, gv, us, cs, lb, uhbt, obc, visc_rem_u, u_cor, bt_cont)
199 do k=1,nz ;
do j=lb%jsh,lb%jeh ;
do i=lb%ish,lb%ieh
200 h(i,j,k) = h(i,j,k) - dt * g%IareaT(i,j) * (uh(i,j,k) - uh(i-1,j,k))
202 if (h(i,j,k) < h_min) h(i,j,k) = h_min
203 enddo ;
enddo ;
enddo
211 subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, LB, uhbt, OBC, &
212 visc_rem_u, u_cor, BT_cont)
215 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
217 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
219 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
222 real,
intent(in) :: dt
227 optional,
pointer :: OBC
228 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
229 optional,
intent(in) :: visc_rem_u
234 real,
dimension(SZIB_(G),SZJ_(G)), &
235 optional,
intent(in) :: uhbt
237 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
238 optional,
intent(out) :: u_cor
245 real,
dimension(SZIB_(G),SZK_(G)) :: duhdu
246 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_L, h_R
247 real,
dimension(SZIB_(G)) :: &
254 logical,
dimension(SZIB_(G)) :: do_I
255 real,
dimension(SZIB_(G),SZK_(G)) :: &
257 real,
dimension(SZIB_(G)) :: FAuI
265 integer :: i, j, k, ish, ieh, jsh, jeh, n, nz
266 logical :: local_specified_BC, use_visc_rem, set_BT_cont, any_simple_OBC
267 logical :: local_Flather_OBC, local_open_BC, is_simple
270 use_visc_rem =
present(visc_rem_u)
271 local_specified_bc = .false. ; set_bt_cont = .false. ; local_flather_obc = .false.
272 local_open_bc = .false.
273 if (
present(bt_cont)) set_bt_cont = (
associated(bt_cont))
274 if (
present(obc))
then ;
if (
associated(obc))
then
275 local_specified_bc = obc%specified_u_BCs_exist_globally
276 local_flather_obc = obc%Flather_u_BCs_exist_globally
277 local_open_bc = obc%open_u_BCs_exist_globally
279 ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh ; nz = g%ke
281 cfl_dt = cs%CFL_limit_adjust / dt
283 if (cs%aggress_adjust) cfl_dt = i_dt
289 if (cs%upwind_1st)
then
290 do j=jsh,jeh ;
do i=ish-1,ieh+1
291 h_l(i,j,k) = h_in(i,j,k) ; h_r(i,j,k) = h_in(i,j,k)
294 call ppm_reconstruction_x(h_in(:,:,k), h_l(:,:,k), h_r(:,:,k), g, lb, &
295 2.0*gv%Angstrom_H, cs%monotonic, simple_2nd=cs%simple_2nd, obc=obc)
297 do i=ish-1,ieh ; visc_rem(i,k) = 1.0 ;
enddo
309 do i=ish-1,ieh ; do_i(i) = .true. ; visc_rem_max(i) = 0.0 ;
enddo
312 if (use_visc_rem)
then ;
do i=ish-1,ieh
313 visc_rem(i,k) = visc_rem_u(i,j,k)
314 visc_rem_max(i) = max(visc_rem_max(i), visc_rem(i,k))
316 call zonal_flux_layer(u(:,j,k), h_in(:,j,k), h_l(:,j,k), h_r(:,j,k), &
317 uh(:,j,k), duhdu(:,k), visc_rem(:,k), &
318 dt, g, us, j, ish, ieh, do_i, cs%vol_CFL, obc)
319 if (local_specified_bc)
then
321 if (obc%segment(obc%segnum_u(i,j))%specified) &
322 uh(i,j,k) = obc%segment(obc%segnum_u(i,j))%normal_trans(i,j,k)
327 if ((.not.use_visc_rem).or.(.not.cs%use_visc_rem_max))
then ;
do i=ish-1,ieh
328 visc_rem_max(i) = 1.0
331 if (
present(uhbt) .or. set_bt_cont)
then
336 if (visc_rem_max(i) > 0.0) i_vrm = 1.0 / visc_rem_max(i)
338 dx_w =
ratio_max(g%areaT(i,j), g%dy_Cu(i,j), 1000.0*g%dxT(i,j))
339 dx_e =
ratio_max(g%areaT(i+1,j), g%dy_Cu(i,j), 1000.0*g%dxT(i+1,j))
340 else ; dx_w = g%dxT(i,j) ; dx_e = g%dxT(i+1,j) ;
endif
341 du_max_cfl(i) = 2.0* (cfl_dt * dx_w) * i_vrm
342 du_min_cfl(i) = -2.0 * (cfl_dt * dx_e) * i_vrm
343 uh_tot_0(i) = 0.0 ; duhdu_tot_0(i) = 0.0
345 do k=1,nz ;
do i=ish-1,ieh
346 duhdu_tot_0(i) = duhdu_tot_0(i) + duhdu(i,k)
347 uh_tot_0(i) = uh_tot_0(i) + uh(i,j,k)
349 if (use_visc_rem)
then
350 if (cs%aggress_adjust)
then
351 do k=1,nz ;
do i=ish-1,ieh
353 dx_w =
ratio_max(g%areaT(i,j), g%dy_Cu(i,j), 1000.0*g%dxT(i,j))
354 dx_e =
ratio_max(g%areaT(i+1,j), g%dy_Cu(i,j), 1000.0*g%dxT(i+1,j))
355 else ; dx_w = g%dxT(i,j) ; dx_e = g%dxT(i+1,j) ;
endif
357 du_lim = 0.499*((dx_w*i_dt - u(i,j,k)) + min(0.0,u(i-1,j,k)))
358 if (du_max_cfl(i) * visc_rem(i,k) > du_lim) &
359 du_max_cfl(i) = du_lim / visc_rem(i,k)
361 du_lim = 0.499*((-dx_e*i_dt - u(i,j,k)) + max(0.0,u(i+1,j,k)))
362 if (du_min_cfl(i) * visc_rem(i,k) < du_lim) &
363 du_min_cfl(i) = du_lim / visc_rem(i,k)
366 do k=1,nz ;
do i=ish-1,ieh
368 dx_w =
ratio_max(g%areaT(i,j), g%dy_Cu(i,j), 1000.0*g%dxT(i,j))
369 dx_e =
ratio_max(g%areaT(i+1,j), g%dy_Cu(i,j), 1000.0*g%dxT(i+1,j))
370 else ; dx_w = g%dxT(i,j) ; dx_e = g%dxT(i+1,j) ;
endif
372 if (du_max_cfl(i) * visc_rem(i,k) > dx_w*cfl_dt - u(i,j,k)) &
373 du_max_cfl(i) = (dx_w*cfl_dt - u(i,j,k)) / visc_rem(i,k)
374 if (du_min_cfl(i) * visc_rem(i,k) < -dx_e*cfl_dt - u(i,j,k)) &
375 du_min_cfl(i) = -(dx_e*cfl_dt + u(i,j,k)) / visc_rem(i,k)
379 if (cs%aggress_adjust)
then
380 do k=1,nz ;
do i=ish-1,ieh
382 dx_w =
ratio_max(g%areaT(i,j), g%dy_Cu(i,j), 1000.0*g%dxT(i,j))
383 dx_e =
ratio_max(g%areaT(i+1,j), g%dy_Cu(i,j), 1000.0*g%dxT(i+1,j))
384 else ; dx_w = g%dxT(i,j) ; dx_e = g%dxT(i+1,j) ;
endif
386 du_max_cfl(i) = min(du_max_cfl(i), 0.499 * &
387 ((dx_w*i_dt - u(i,j,k)) + min(0.0,u(i-1,j,k))) )
388 du_min_cfl(i) = max(du_min_cfl(i), 0.499 * &
389 ((-dx_e*i_dt - u(i,j,k)) + max(0.0,u(i+1,j,k))) )
392 do k=1,nz ;
do i=ish-1,ieh
394 dx_w =
ratio_max(g%areaT(i,j), g%dy_Cu(i,j), 1000.0*g%dxT(i,j))
395 dx_e =
ratio_max(g%areaT(i+1,j), g%dy_Cu(i,j), 1000.0*g%dxT(i+1,j))
396 else ; dx_w = g%dxT(i,j) ; dx_e = g%dxT(i+1,j) ;
endif
398 du_max_cfl(i) = min(du_max_cfl(i), dx_w*cfl_dt - u(i,j,k))
399 du_min_cfl(i) = max(du_min_cfl(i), -(dx_e*cfl_dt + u(i,j,k)))
404 du_max_cfl(i) = max(du_max_cfl(i),0.0)
405 du_min_cfl(i) = min(du_min_cfl(i),0.0)
408 any_simple_obc = .false.
409 if (
present(uhbt) .or. set_bt_cont)
then
410 if (local_specified_bc .or. local_flather_obc)
then ;
do i=ish-1,ieh
412 is_simple = obc%segment(obc%segnum_u(i,j))%specified
413 do_i(i) = .not.(obc%segnum_u(i,j) /= obc_none .and. is_simple)
414 any_simple_obc = any_simple_obc .or. is_simple
415 enddo ;
else ;
do i=ish-1,ieh
420 if (
present(uhbt))
then
421 call zonal_flux_adjust(u, h_in, h_l, h_r, uhbt(:,j), uh_tot_0, duhdu_tot_0, du, &
422 du_max_cfl, du_min_cfl, dt, g, us, cs, visc_rem, &
423 j, ish, ieh, do_i, .true., uh, obc=obc)
425 if (
present(u_cor))
then ;
do k=1,nz
426 do i=ish-1,ieh ; u_cor(i,j,k) = u(i,j,k) + du(i) * visc_rem(i,k) ;
enddo
427 if (local_specified_bc)
then ;
do i=ish-1,ieh
428 if (obc%segment(obc%segnum_u(i,j))%specified) &
429 u_cor(i,j,k) = obc%segment(obc%segnum_u(i,j))%normal_vel(i,j,k)
435 if (set_bt_cont)
then
437 du_max_cfl, du_min_cfl, dt, g, us, cs, visc_rem, &
438 visc_rem_max, j, ish, ieh, do_i)
439 if (any_simple_obc)
then
441 do_i(i) = obc%segment(obc%segnum_u(i,j))%specified
442 if (do_i(i)) faui(i) = gv%H_subroundoff*g%dy_Cu(i,j)
444 do k=1,nz ;
do i=ish-1,ieh ;
if (do_i(i))
then
445 if ((abs(obc%segment(obc%segnum_u(i,j))%normal_vel(i,j,k)) > 0.0) .and. &
446 (obc%segment(obc%segnum_u(i,j))%specified)) &
447 faui(i) = faui(i) + obc%segment(obc%segnum_u(i,j))%normal_trans(i,j,k) / &
448 obc%segment(obc%segnum_u(i,j))%normal_vel(i,j,k)
449 endif ;
enddo ;
enddo
450 do i=ish-1,ieh ;
if (do_i(i))
then
451 bt_cont%FA_u_W0(i,j) = faui(i) ; bt_cont%FA_u_E0(i,j) = faui(i)
452 bt_cont%FA_u_WW(i,j) = faui(i) ; bt_cont%FA_u_EE(i,j) = faui(i)
453 bt_cont%uBT_WW(i,j) = 0.0 ; bt_cont%uBT_EE(i,j) = 0.0
462 if (local_open_bc .and. set_bt_cont)
then
463 do n = 1, obc%number_of_segments
464 if (obc%segment(n)%open .and. obc%segment(n)%is_E_or_W)
then
465 i = obc%segment(n)%HI%IsdB
467 do j = obc%segment(n)%HI%Jsd, obc%segment(n)%HI%Jed
469 do k=1,nz ; fa_u = fa_u + h_in(i,j,k)*g%dy_Cu(i,j) ;
enddo
470 bt_cont%FA_u_W0(i,j) = fa_u ; bt_cont%FA_u_E0(i,j) = fa_u
471 bt_cont%FA_u_WW(i,j) = fa_u ; bt_cont%FA_u_EE(i,j) = fa_u
472 bt_cont%uBT_WW(i,j) = 0.0 ; bt_cont%uBT_EE(i,j) = 0.0
475 do j = obc%segment(n)%HI%Jsd, obc%segment(n)%HI%Jed
477 do k=1,nz ; fa_u = fa_u + h_in(i+1,j,k)*g%dy_Cu(i,j) ;
enddo
478 bt_cont%FA_u_W0(i,j) = fa_u ; bt_cont%FA_u_E0(i,j) = fa_u
479 bt_cont%FA_u_WW(i,j) = fa_u ; bt_cont%FA_u_EE(i,j) = fa_u
480 bt_cont%uBT_WW(i,j) = 0.0 ; bt_cont%uBT_EE(i,j) = 0.0
488 if (set_bt_cont)
then ;
if (
allocated(bt_cont%h_u))
then
489 if (
present(u_cor))
then
491 cs%vol_CFL, cs%marginal_faces, visc_rem_u, obc)
494 cs%vol_CFL, cs%marginal_faces, visc_rem_u, obc)
501 subroutine zonal_flux_layer(u, h, h_L, h_R, uh, duhdu, visc_rem, dt, G, US, j, &
502 ish, ieh, do_I, vol_CFL, OBC)
504 real,
dimension(SZIB_(G)),
intent(in) :: u
505 real,
dimension(SZIB_(G)),
intent(in) :: visc_rem
510 real,
dimension(SZI_(G)),
intent(in) :: h
511 real,
dimension(SZI_(G)),
intent(in) :: h_L
512 real,
dimension(SZI_(G)),
intent(in) :: h_R
513 real,
dimension(SZIB_(G)),
intent(inout) :: uh
515 real,
dimension(SZIB_(G)),
intent(inout) :: duhdu
517 real,
intent(in) :: dt
519 integer,
intent(in) :: j
520 integer,
intent(in) :: ish
521 integer,
intent(in) :: ieh
522 logical,
dimension(SZIB_(G)),
intent(in) :: do_I
523 logical,
intent(in) :: vol_CFL
532 logical :: local_open_BC
534 local_open_bc = .false.
535 if (
present(obc))
then ;
if (
associated(obc))
then
536 local_open_bc = obc%open_u_BCs_exist_globally
539 do i=ish-1,ieh ;
if (do_i(i))
then
542 if (vol_cfl)
then ; cfl = (u(i) * dt) * (g%dy_Cu(i,j) * g%IareaT(i,j))
543 else ; cfl = u(i) * dt * g%IdxT(i,j) ;
endif
544 curv_3 = h_l(i) + h_r(i) - 2.0*h(i)
545 uh(i) = g%dy_Cu(i,j) * u(i) * &
546 (h_r(i) + cfl * (0.5*(h_l(i) - h_r(i)) + curv_3*(cfl - 1.5)))
547 h_marg = h_r(i) + cfl * ((h_l(i) - h_r(i)) + 3.0*curv_3*(cfl - 1.0))
548 elseif (u(i) < 0.0)
then
549 if (vol_cfl)
then ; cfl = (-u(i) * dt) * (g%dy_Cu(i,j) * g%IareaT(i+1,j))
550 else ; cfl = -u(i) * dt * g%IdxT(i+1,j) ;
endif
551 curv_3 = h_l(i+1) + h_r(i+1) - 2.0*h(i+1)
552 uh(i) = g%dy_Cu(i,j) * u(i) * &
553 (h_l(i+1) + cfl * (0.5*(h_r(i+1)-h_l(i+1)) + curv_3*(cfl - 1.5)))
554 h_marg = h_l(i+1) + cfl * ((h_r(i+1)-h_l(i+1)) + 3.0*curv_3*(cfl - 1.0))
557 h_marg = 0.5 * (h_l(i+1) + h_r(i))
559 duhdu(i) = g%dy_Cu(i,j) * h_marg * visc_rem(i)
562 if (local_open_bc)
then
563 do i=ish-1,ieh ;
if (do_i(i))
then
564 if (obc%segment(obc%segnum_u(i,j))%open)
then
566 uh(i) = g%dy_Cu(i,j) * u(i) * h(i)
567 duhdu(i) = g%dy_Cu(i,j) * h(i) * visc_rem(i)
569 uh(i) = g%dy_Cu(i,j) * u(i) * h(i+1)
570 duhdu(i) = g%dy_Cu(i,j) * h(i+1) * visc_rem(i)
578 subroutine zonal_face_thickness(u, h, h_L, h_R, h_u, dt, G, US, LB, vol_CFL, &
579 marginal, visc_rem_u, OBC)
581 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(in) :: u
582 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
584 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h_L
586 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h_R
588 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(inout) :: h_u
589 real,
intent(in) :: dt
592 logical,
intent(in) :: vol_CFL
594 logical,
intent(in) :: marginal
596 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
597 optional,
intent(in) :: visc_rem_u
610 logical :: local_open_BC
611 integer :: i, j, k, ish, ieh, jsh, jeh, nz, n
612 ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh ; nz = g%ke
615 do k=1,nz ;
do j=jsh,jeh ;
do i=ish-1,ieh
616 if (u(i,j,k) > 0.0)
then
617 if (vol_cfl)
then ; cfl = (u(i,j,k) * dt) * (g%dy_Cu(i,j) * g%IareaT(i,j))
618 else ; cfl = u(i,j,k) * dt * g%IdxT(i,j) ;
endif
619 curv_3 = h_l(i,j,k) + h_r(i,j,k) - 2.0*h(i,j,k)
620 h_avg = h_r(i,j,k) + cfl * (0.5*(h_l(i,j,k) - h_r(i,j,k)) + curv_3*(cfl - 1.5))
621 h_marg = h_r(i,j,k) + cfl * ((h_l(i,j,k) - h_r(i,j,k)) + 3.0*curv_3*(cfl - 1.0))
622 elseif (u(i,j,k) < 0.0)
then
623 if (vol_cfl)
then ; cfl = (-u(i,j,k)*dt) * (g%dy_Cu(i,j) * g%IareaT(i+1,j))
624 else ; cfl = -u(i,j,k) * dt * g%IdxT(i+1,j) ;
endif
625 curv_3 = h_l(i+1,j,k) + h_r(i+1,j,k) - 2.0*h(i+1,j,k)
626 h_avg = h_l(i+1,j,k) + cfl * (0.5*(h_r(i+1,j,k)-h_l(i+1,j,k)) + curv_3*(cfl - 1.5))
627 h_marg = h_l(i+1,j,k) + cfl * ((h_r(i+1,j,k)-h_l(i+1,j,k)) + &
628 3.0*curv_3*(cfl - 1.0))
630 h_avg = 0.5 * (h_l(i+1,j,k) + h_r(i,j,k))
633 h_marg = 0.5 * (h_l(i+1,j,k) + h_r(i,j,k))
638 if (marginal)
then ; h_u(i,j,k) = h_marg
639 else ; h_u(i,j,k) = h_avg ;
endif
640 enddo ;
enddo ;
enddo
641 if (
present(visc_rem_u))
then
643 do k=1,nz ;
do j=jsh,jeh ;
do i=ish-1,ieh
644 h_u(i,j,k) = h_u(i,j,k) * visc_rem_u(i,j,k)
645 enddo ;
enddo ;
enddo
648 local_open_bc = .false.
649 if (
present(obc))
then ;
if (
associated(obc))
then
650 local_open_bc = obc%open_u_BCs_exist_globally
652 if (local_open_bc)
then
653 do n = 1, obc%number_of_segments
654 if (obc%segment(n)%open .and. obc%segment(n)%is_E_or_W)
then
655 i = obc%segment(n)%HI%IsdB
657 if (
present(visc_rem_u))
then ;
do k=1,nz
658 do j = obc%segment(n)%HI%jsd, obc%segment(n)%HI%jed
659 h_u(i,j,k) = h(i,j,k) * visc_rem_u(i,j,k)
661 enddo ;
else ;
do k=1,nz
662 do j = obc%segment(n)%HI%jsd, obc%segment(n)%HI%jed
663 h_u(i,j,k) = h(i,j,k)
667 if (
present(visc_rem_u))
then ;
do k=1,nz
668 do j = obc%segment(n)%HI%jsd, obc%segment(n)%HI%jed
669 h_u(i,j,k) = h(i+1,j,k) * visc_rem_u(i,j,k)
671 enddo ;
else ;
do k=1,nz
672 do j = obc%segment(n)%HI%jsd, obc%segment(n)%HI%jed
673 h_u(i,j,k) = h(i+1,j,k)
686 du, du_max_CFL, du_min_CFL, dt, G, US, CS, visc_rem, &
687 j, ish, ieh, do_I_in, full_precision, uh_3d, OBC)
689 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(in) :: u
690 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h_in
692 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h_L
694 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h_R
696 real,
dimension(SZIB_(G),SZK_(G)),
intent(in) :: visc_rem
701 real,
dimension(SZIB_(G)),
optional,
intent(in) :: uhbt
704 real,
dimension(SZIB_(G)),
intent(in) :: du_max_CFL
706 real,
dimension(SZIB_(G)),
intent(in) :: du_min_CFL
708 real,
dimension(SZIB_(G)),
intent(in) :: uh_tot_0
710 real,
dimension(SZIB_(G)),
intent(in) :: duhdu_tot_0
712 real,
dimension(SZIB_(G)),
intent(out) :: du
714 real,
intent(in) :: dt
717 integer,
intent(in) :: j
718 integer,
intent(in) :: ish
719 integer,
intent(in) :: ieh
720 logical,
dimension(SZIB_(G)),
intent(in) :: do_I_in
722 logical,
optional,
intent(in) :: full_precision
725 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
optional,
intent(inout) :: uh_3d
729 real,
dimension(SZIB_(G),SZK_(G)) :: &
732 real,
dimension(SZIB_(G)) :: &
743 integer :: i, k, nz, itt, max_itts = 20
744 logical :: full_prec, domore, do_I(SZIB_(G))
747 full_prec = .true. ;
if (
present(full_precision)) full_prec = full_precision
749 uh_aux(:,:) = 0.0 ; duhdu(:,:) = 0.0
751 if (
present(uh_3d))
then ;
do k=1,nz ;
do i=ish-1,ieh
752 uh_aux(i,k) = uh_3d(i,j,k)
753 enddo ;
enddo ;
endif
756 du(i) = 0.0 ; do_i(i) = do_i_in(i)
757 du_max(i) = du_max_cfl(i) ; du_min(i) = du_min_cfl(i)
758 uh_err(i) = uh_tot_0(i) - uhbt(i) ; duhdu_tot(i) = duhdu_tot_0(i)
759 uh_err_best(i) = abs(uh_err(i))
765 case (:1) ; tol_eta = 1e-6 * cs%tol_eta
766 case (2) ; tol_eta = 1e-4 * cs%tol_eta
767 case (3) ; tol_eta = 1e-2 * cs%tol_eta
768 case default ; tol_eta = cs%tol_eta
771 tol_eta = cs%tol_eta_aux ;
if (itt<=1) tol_eta = 1e-6 * cs%tol_eta_aux
776 if (uh_err(i) > 0.0)
then ; du_max(i) = du(i)
777 elseif (uh_err(i) < 0.0)
then ; du_min(i) = du(i)
778 else ; do_i(i) = .false. ;
endif
781 do i=ish-1,ieh ;
if (do_i(i))
then
782 if ((dt * min(g%IareaT(i,j),g%IareaT(i+1,j))*abs(uh_err(i)) > tol_eta) .or. &
783 (cs%better_iter .and. ((abs(uh_err(i)) > tol_vel * duhdu_tot(i)) .or. &
784 (abs(uh_err(i)) > uh_err_best(i))) ))
then
788 ddu = -uh_err(i) / duhdu_tot(i)
791 if (abs(ddu) < 1.0e-15*abs(du(i)))
then
793 elseif (ddu > 0.0)
then
794 if (du(i) >= du_max(i))
then
795 du(i) = 0.5*(du_prev + du_max(i))
796 if (du_max(i) - du_prev < 1.0e-15*abs(du(i))) do_i(i) = .false.
799 if (du(i) <= du_min(i))
then
800 du(i) = 0.5*(du_prev + du_min(i))
801 if (du_prev - du_min(i) < 1.0e-15*abs(du(i))) do_i(i) = .false.
806 du(i) = du(i) - uh_err(i) / duhdu_tot(i)
807 if ((du(i) >= du_max(i)) .or. (du(i) <= du_min(i))) &
808 du(i) = 0.5*(du_max(i) + du_min(i))
810 if (do_i(i)) domore = .true.
815 if (.not.domore)
exit
817 if ((itt < max_itts) .or.
present(uh_3d))
then ;
do k=1,nz
818 do i=ish-1,ieh ; u_new(i) = u(i,j,k) + du(i) * visc_rem(i,k) ;
enddo
819 call zonal_flux_layer(u_new, h_in(:,j,k), h_l(:,j,k), h_r(:,j,k), &
820 uh_aux(:,k), duhdu(:,k), visc_rem(:,k), &
821 dt, g, us, j, ish, ieh, do_i, cs%vol_CFL, obc)
824 if (itt < max_itts)
then
826 uh_err(i) = -uhbt(i) ; duhdu_tot(i) = 0.0
828 do k=1,nz ;
do i=ish-1,ieh
829 uh_err(i) = uh_err(i) + uh_aux(i,k)
830 duhdu_tot(i) = duhdu_tot(i) + duhdu(i,k)
833 uh_err_best(i) = min(uh_err_best(i), abs(uh_err(i)))
841 if (
present(uh_3d))
then ;
do k=1,nz ;
do i=ish-1,ieh
842 uh_3d(i,j,k) = uh_aux(i,k)
843 enddo ;
enddo ;
endif
849 subroutine set_zonal_bt_cont(u, h_in, h_L, h_R, BT_cont, uh_tot_0, duhdu_tot_0, &
850 du_max_CFL, du_min_CFL, dt, G, US, CS, visc_rem, &
851 visc_rem_max, j, ish, ieh, do_I)
853 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(in) :: u
854 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h_in
856 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h_L
858 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h_R
862 real,
dimension(SZIB_(G)),
intent(in) :: uh_tot_0
864 real,
dimension(SZIB_(G)),
intent(in) :: duhdu_tot_0
866 real,
dimension(SZIB_(G)),
intent(in) :: du_max_CFL
868 real,
dimension(SZIB_(G)),
intent(in) :: du_min_CFL
870 real,
intent(in) :: dt
873 real,
dimension(SZIB_(G),SZK_(G)),
intent(in) :: visc_rem
878 real,
dimension(SZIB_(G)),
intent(in) :: visc_rem_max
879 integer,
intent(in) :: j
880 integer,
intent(in) :: ish
881 integer,
intent(in) :: ieh
882 logical,
dimension(SZIB_(G)),
intent(in) :: do_I
885 real,
dimension(SZIB_(G)) :: &
918 nz = g%ke ; idt = 1.0 / dt
919 min_visc_rem = 0.1 ; cfl_min = 1e-6
922 do i=ish-1,ieh ; zeros(i) = 0.0 ;
enddo
924 du_max_cfl, du_min_cfl, dt, g, us, cs, visc_rem, &
925 j, ish, ieh, do_i, .true.)
932 if (do_i(i)) domore = .true.
933 du_cfl(i) = (cfl_min * idt) * g%dxCu(i,j)
934 dur(i) = min(0.0,du0(i) - du_cfl(i))
935 dul(i) = max(0.0,du0(i) + du_cfl(i))
936 famt_l(i) = 0.0 ; famt_r(i) = 0.0 ; famt_0(i) = 0.0
937 uhtot_l(i) = 0.0 ; uhtot_r(i) = 0.0
940 if (.not.domore)
then
941 do k=1,nz ;
do i=ish-1,ieh
942 bt_cont%FA_u_W0(i,j) = 0.0 ; bt_cont%FA_u_WW(i,j) = 0.0
943 bt_cont%FA_u_E0(i,j) = 0.0 ; bt_cont%FA_u_EE(i,j) = 0.0
944 bt_cont%uBT_WW(i,j) = 0.0 ; bt_cont%uBT_EE(i,j) = 0.0
949 do k=1,nz ;
do i=ish-1,ieh ;
if (do_i(i))
then
950 visc_rem_lim = max(visc_rem(i,k), min_visc_rem*visc_rem_max(i))
951 if (visc_rem_lim > 0.0)
then
952 if (u(i,j,k) + dur(i)*visc_rem_lim > -du_cfl(i)*visc_rem(i,k)) &
953 dur(i) = -(u(i,j,k) + du_cfl(i)*visc_rem(i,k)) / visc_rem_lim
954 if (u(i,j,k) + dul(i)*visc_rem_lim < du_cfl(i)*visc_rem(i,k)) &
955 dul(i) = -(u(i,j,k) - du_cfl(i)*visc_rem(i,k)) / visc_rem_lim
957 endif ;
enddo ;
enddo
960 do i=ish-1,ieh ;
if (do_i(i))
then
961 u_l(i) = u(i,j,k) + dul(i) * visc_rem(i,k)
962 u_r(i) = u(i,j,k) + dur(i) * visc_rem(i,k)
963 u_0(i) = u(i,j,k) + du0(i) * visc_rem(i,k)
965 call zonal_flux_layer(u_0, h_in(:,j,k), h_l(:,j,k), h_r(:,j,k), uh_0, duhdu_0, &
966 visc_rem(:,k), dt, g, us, j, ish, ieh, do_i, cs%vol_CFL)
967 call zonal_flux_layer(u_l, h_in(:,j,k), h_l(:,j,k), h_r(:,j,k), uh_l, duhdu_l, &
968 visc_rem(:,k), dt, g, us, j, ish, ieh, do_i, cs%vol_CFL)
969 call zonal_flux_layer(u_r, h_in(:,j,k), h_l(:,j,k), h_r(:,j,k), uh_r, duhdu_r, &
970 visc_rem(:,k), dt, g, us, j, ish, ieh, do_i, cs%vol_CFL)
971 do i=ish-1,ieh ;
if (do_i(i))
then
972 famt_0(i) = famt_0(i) + duhdu_0(i)
973 famt_l(i) = famt_l(i) + duhdu_l(i)
974 famt_r(i) = famt_r(i) + duhdu_r(i)
975 uhtot_l(i) = uhtot_l(i) + uh_l(i)
976 uhtot_r(i) = uhtot_r(i) + uh_r(i)
979 do i=ish-1,ieh ;
if (do_i(i))
then
980 fa_0 = famt_0(i) ; fa_avg = famt_0(i)
981 if ((dul(i) - du0(i)) /= 0.0) &
982 fa_avg = uhtot_l(i) / (dul(i) - du0(i))
983 if (fa_avg > max(fa_0, famt_l(i)))
then ; fa_avg = max(fa_0, famt_l(i))
984 elseif (fa_avg < min(fa_0, famt_l(i)))
then ; fa_0 = fa_avg ;
endif
986 bt_cont%FA_u_W0(i,j) = fa_0 ; bt_cont%FA_u_WW(i,j) = famt_l(i)
987 if (abs(fa_0-famt_l(i)) <= 1e-12*fa_0)
then ; bt_cont%uBT_WW(i,j) = 0.0 ;
else
988 bt_cont%uBT_WW(i,j) = (1.5 * (dul(i) - du0(i))) * &
989 ((famt_l(i) - fa_avg) / (famt_l(i) - fa_0))
992 fa_0 = famt_0(i) ; fa_avg = famt_0(i)
993 if ((dur(i) - du0(i)) /= 0.0) &
994 fa_avg = uhtot_r(i) / (dur(i) - du0(i))
995 if (fa_avg > max(fa_0, famt_r(i)))
then ; fa_avg = max(fa_0, famt_r(i))
996 elseif (fa_avg < min(fa_0, famt_r(i)))
then ; fa_0 = fa_avg ;
endif
998 bt_cont%FA_u_E0(i,j) = fa_0 ; bt_cont%FA_u_EE(i,j) = famt_r(i)
999 if (abs(famt_r(i) - fa_0) <= 1e-12*fa_0)
then ; bt_cont%uBT_EE(i,j) = 0.0 ;
else
1000 bt_cont%uBT_EE(i,j) = (1.5 * (dur(i) - du0(i))) * &
1001 ((famt_r(i) - fa_avg) / (famt_r(i) - fa_0))
1004 bt_cont%FA_u_W0(i,j) = 0.0 ; bt_cont%FA_u_WW(i,j) = 0.0
1005 bt_cont%FA_u_E0(i,j) = 0.0 ; bt_cont%FA_u_EE(i,j) = 0.0
1006 bt_cont%uBT_WW(i,j) = 0.0 ; bt_cont%uBT_EE(i,j) = 0.0
1012 subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, US, CS, LB, vhbt, OBC, &
1013 visc_rem_v, v_cor, BT_cont)
1016 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(in) :: v
1017 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h_in
1019 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(out) :: vh
1021 real,
intent(in) :: dt
1027 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1028 optional,
intent(in) :: visc_rem_v
1033 real,
dimension(SZI_(G),SZJB_(G)),
optional,
intent(in) :: vhbt
1035 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1036 optional,
intent(out) :: v_cor
1042 real,
dimension(SZI_(G),SZK_(G)) :: &
1044 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
1046 real,
dimension(SZI_(G)) :: &
1053 logical,
dimension(SZI_(G)) :: do_I
1054 real,
dimension(SZI_(G)) :: FAvi
1056 real,
dimension(SZI_(G),SZK_(G)) :: &
1064 integer :: i, j, k, ish, ieh, jsh, jeh, n, nz
1065 logical :: local_specified_BC, use_visc_rem, set_BT_cont, any_simple_OBC
1066 logical :: local_Flather_OBC, is_simple, local_open_BC
1069 use_visc_rem =
present(visc_rem_v)
1070 local_specified_bc = .false. ; set_bt_cont = .false. ; local_flather_obc = .false.
1071 local_open_bc = .false.
1072 if (
present(bt_cont)) set_bt_cont = (
associated(bt_cont))
1073 if (
present(obc))
then ;
if (
associated(obc))
then ;
if (obc%OBC_pe)
then
1074 local_specified_bc = obc%specified_v_BCs_exist_globally
1075 local_flather_obc = obc%Flather_v_BCs_exist_globally
1076 local_open_bc = obc%open_v_BCs_exist_globally
1077 endif ;
endif ;
endif
1078 ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh ; nz = g%ke
1080 cfl_dt = cs%CFL_limit_adjust / dt
1082 if (cs%aggress_adjust) cfl_dt = i_dt
1088 if (cs%upwind_1st)
then
1089 do j=jsh-1,jeh+1 ;
do i=ish,ieh
1090 h_l(i,j,k) = h_in(i,j,k) ; h_r(i,j,k) = h_in(i,j,k)
1093 call ppm_reconstruction_y(h_in(:,:,k), h_l(:,:,k), h_r(:,:,k), g, lb, &
1094 2.0*gv%Angstrom_H, cs%monotonic, simple_2nd=cs%simple_2nd, obc=obc)
1096 do i=ish,ieh ; visc_rem(i,k) = 1.0 ;
enddo
1109 do i=ish,ieh ; do_i(i) = .true. ; visc_rem_max(i) = 0.0 ;
enddo
1112 if (use_visc_rem)
then ;
do i=ish,ieh
1113 visc_rem(i,k) = visc_rem_v(i,j,k)
1114 visc_rem_max(i) = max(visc_rem_max(i), visc_rem(i,k))
1116 call merid_flux_layer(v(:,j,k), h_in(:,:,k), h_l(:,:,k), h_r(:,:,k), &
1117 vh(:,j,k), dvhdv(:,k), visc_rem(:,k), &
1118 dt, g, us, j, ish, ieh, do_i, cs%vol_CFL, obc)
1119 if (local_specified_bc)
then
1121 if (obc%segment(obc%segnum_v(i,j))%specified) &
1122 vh(i,j,k) = obc%segment(obc%segnum_v(i,j))%normal_trans(i,j,k)
1126 if ((.not.use_visc_rem) .or. (.not.cs%use_visc_rem_max))
then ;
do i=ish,ieh
1127 visc_rem_max(i) = 1.0
1130 if (
present(vhbt) .or. set_bt_cont)
then
1135 if (visc_rem_max(i) > 0.0) i_vrm = 1.0 / visc_rem_max(i)
1136 if (cs%vol_CFL)
then
1137 dy_s =
ratio_max(g%areaT(i,j), g%dx_Cv(i,j), 1000.0*g%dyT(i,j))
1138 dy_n =
ratio_max(g%areaT(i,j+1), g%dx_Cv(i,j), 1000.0*g%dyT(i,j+1))
1139 else ; dy_s = g%dyT(i,j) ; dy_n = g%dyT(i,j+1) ;
endif
1140 dv_max_cfl(i) = 2.0 * (cfl_dt * dy_s) * i_vrm
1141 dv_min_cfl(i) = -2.0 * (cfl_dt * dy_n) * i_vrm
1142 vh_tot_0(i) = 0.0 ; dvhdv_tot_0(i) = 0.0
1144 do k=1,nz ;
do i=ish,ieh
1145 dvhdv_tot_0(i) = dvhdv_tot_0(i) + dvhdv(i,k)
1146 vh_tot_0(i) = vh_tot_0(i) + vh(i,j,k)
1149 if (use_visc_rem)
then
1150 if (cs%aggress_adjust)
then
1151 do k=1,nz ;
do i=ish,ieh
1152 if (cs%vol_CFL)
then
1153 dy_s =
ratio_max(g%areaT(i,j), g%dx_Cv(i,j), 1000.0*g%dyT(i,j))
1154 dy_n =
ratio_max(g%areaT(i,j+1), g%dx_Cv(i,j), 1000.0*g%dyT(i,j+1))
1155 else ; dy_s = g%dyT(i,j) ; dy_n = g%dyT(i,j+1) ;
endif
1156 dv_lim = 0.499*((dy_s*i_dt - v(i,j,k)) + min(0.0,v(i,j-1,k)))
1157 if (dv_max_cfl(i) * visc_rem(i,k) > dv_lim) &
1158 dv_max_cfl(i) = dv_lim / visc_rem(i,k)
1160 dv_lim = 0.499*((-dy_n*cfl_dt - v(i,j,k)) + max(0.0,v(i,j+1,k)))
1161 if (dv_min_cfl(i) * visc_rem(i,k) < dv_lim) &
1162 dv_min_cfl(i) = dv_lim / visc_rem(i,k)
1165 do k=1,nz ;
do i=ish,ieh
1166 if (cs%vol_CFL)
then
1167 dy_s =
ratio_max(g%areaT(i,j), g%dx_Cv(i,j), 1000.0*g%dyT(i,j))
1168 dy_n =
ratio_max(g%areaT(i,j+1), g%dx_Cv(i,j), 1000.0*g%dyT(i,j+1))
1169 else ; dy_s = g%dyT(i,j) ; dy_n = g%dyT(i,j+1) ;
endif
1170 if (dv_max_cfl(i) * visc_rem(i,k) > dy_s*cfl_dt - v(i,j,k)) &
1171 dv_max_cfl(i) = (dy_s*cfl_dt - v(i,j,k)) / visc_rem(i,k)
1172 if (dv_min_cfl(i) * visc_rem(i,k) < -dy_n*cfl_dt - v(i,j,k)) &
1173 dv_min_cfl(i) = -(dy_n*cfl_dt + v(i,j,k)) / visc_rem(i,k)
1177 if (cs%aggress_adjust)
then
1178 do k=1,nz ;
do i=ish,ieh
1179 if (cs%vol_CFL)
then
1180 dy_s =
ratio_max(g%areaT(i,j), g%dx_Cv(i,j), 1000.0*g%dyT(i,j))
1181 dy_n =
ratio_max(g%areaT(i,j+1), g%dx_Cv(i,j), 1000.0*g%dyT(i,j+1))
1182 else ; dy_s = g%dyT(i,j) ; dy_n = g%dyT(i,j+1) ;
endif
1183 dv_max_cfl(i) = min(dv_max_cfl(i), 0.499 * &
1184 ((dy_s*i_dt - v(i,j,k)) + min(0.0,v(i,j-1,k))) )
1185 dv_min_cfl(i) = max(dv_min_cfl(i), 0.499 * &
1186 ((-dy_n*i_dt - v(i,j,k)) + max(0.0,v(i,j+1,k))) )
1189 do k=1,nz ;
do i=ish,ieh
1190 if (cs%vol_CFL)
then
1191 dy_s =
ratio_max(g%areaT(i,j), g%dx_Cv(i,j), 1000.0*g%dyT(i,j))
1192 dy_n =
ratio_max(g%areaT(i,j+1), g%dx_Cv(i,j), 1000.0*g%dyT(i,j+1))
1193 else ; dy_s = g%dyT(i,j) ; dy_n = g%dyT(i,j+1) ;
endif
1194 dv_max_cfl(i) = min(dv_max_cfl(i), dy_s*cfl_dt - v(i,j,k))
1195 dv_min_cfl(i) = max(dv_min_cfl(i), -(dy_n*cfl_dt + v(i,j,k)))
1200 dv_max_cfl(i) = max(dv_max_cfl(i),0.0)
1201 dv_min_cfl(i) = min(dv_min_cfl(i),0.0)
1204 any_simple_obc = .false.
1205 if (
present(vhbt) .or. set_bt_cont)
then
1206 if (local_specified_bc .or. local_flather_obc)
then ;
do i=ish,ieh
1208 is_simple = obc%segment(obc%segnum_v(i,j))%specified
1209 do_i(i) = .not.(obc%segnum_v(i,j) /= obc_none .and. is_simple)
1210 any_simple_obc = any_simple_obc .or. is_simple
1211 enddo ;
else ;
do i=ish,ieh
1216 if (
present(vhbt))
then
1218 dv_max_cfl, dv_min_cfl, dt, g, us, cs, visc_rem, &
1219 j, ish, ieh, do_i, .true., vh, obc=obc)
1221 if (
present(v_cor))
then ;
do k=1,nz
1222 do i=ish,ieh ; v_cor(i,j,k) = v(i,j,k) + dv(i) * visc_rem(i,k) ;
enddo
1223 if (local_specified_bc)
then ;
do i=ish,ieh
1224 if (obc%segment(obc%segnum_v(i,j))%specified) &
1225 v_cor(i,j,k) = obc%segment(obc%segnum_v(i,j))%normal_vel(i,j,k)
1230 if (set_bt_cont)
then
1232 dv_max_cfl, dv_min_cfl, dt, g, us, cs, visc_rem, &
1233 visc_rem_max, j, ish, ieh, do_i)
1234 if (any_simple_obc)
then
1236 do_i(i) = (obc%segment(obc%segnum_v(i,j))%specified)
1237 if (do_i(i)) favi(i) = gv%H_subroundoff*g%dx_Cv(i,j)
1239 do k=1,nz ;
do i=ish,ieh ;
if (do_i(i))
then
1240 if ((abs(obc%segment(obc%segnum_v(i,j))%normal_vel(i,j,k)) > 0.0) .and. &
1241 (obc%segment(obc%segnum_v(i,j))%specified)) &
1242 favi(i) = favi(i) + obc%segment(obc%segnum_v(i,j))%normal_trans(i,j,k) / &
1243 obc%segment(obc%segnum_v(i,j))%normal_vel(i,j,k)
1244 endif ;
enddo ;
enddo
1245 do i=ish,ieh ;
if (do_i(i))
then
1246 bt_cont%FA_v_S0(i,j) = favi(i) ; bt_cont%FA_v_N0(i,j) = favi(i)
1247 bt_cont%FA_v_SS(i,j) = favi(i) ; bt_cont%FA_v_NN(i,j) = favi(i)
1248 bt_cont%vBT_SS(i,j) = 0.0 ; bt_cont%vBT_NN(i,j) = 0.0
1257 if (local_open_bc .and. set_bt_cont)
then
1258 do n = 1, obc%number_of_segments
1259 if (obc%segment(n)%open .and. obc%segment(n)%is_N_or_S)
then
1260 j = obc%segment(n)%HI%JsdB
1262 do i = obc%segment(n)%HI%Isd, obc%segment(n)%HI%Ied
1264 do k=1,nz ; fa_v = fa_v + h_in(i,j,k)*g%dx_Cv(i,j) ;
enddo
1265 bt_cont%FA_v_S0(i,j) = fa_v ; bt_cont%FA_v_N0(i,j) = fa_v
1266 bt_cont%FA_v_SS(i,j) = fa_v ; bt_cont%FA_v_NN(i,j) = fa_v
1267 bt_cont%vBT_SS(i,j) = 0.0 ; bt_cont%vBT_NN(i,j) = 0.0
1270 do i = obc%segment(n)%HI%Isd, obc%segment(n)%HI%Ied
1272 do k=1,nz ; fa_v = fa_v + h_in(i,j+1,k)*g%dx_Cv(i,j) ;
enddo
1273 bt_cont%FA_v_S0(i,j) = fa_v ; bt_cont%FA_v_N0(i,j) = fa_v
1274 bt_cont%FA_v_SS(i,j) = fa_v ; bt_cont%FA_v_NN(i,j) = fa_v
1275 bt_cont%vBT_SS(i,j) = 0.0 ; bt_cont%vBT_NN(i,j) = 0.0
1283 if (set_bt_cont)
then ;
if (
allocated(bt_cont%h_v))
then
1284 if (
present(v_cor))
then
1286 cs%vol_CFL, cs%marginal_faces, visc_rem_v, obc)
1289 cs%vol_CFL, cs%marginal_faces, visc_rem_v, obc)
1296 subroutine merid_flux_layer(v, h, h_L, h_R, vh, dvhdv, visc_rem, dt, G, US, J, &
1297 ish, ieh, do_I, vol_CFL, OBC)
1299 real,
dimension(SZI_(G)),
intent(in) :: v
1300 real,
dimension(SZI_(G)),
intent(in) :: visc_rem
1305 real,
dimension(SZI_(G),SZJ_(G)),
intent(in) :: h
1307 real,
dimension(SZI_(G),SZJ_(G)),
intent(in) :: h_L
1309 real,
dimension(SZI_(G),SZJ_(G)),
intent(in) :: h_R
1311 real,
dimension(SZI_(G)),
intent(inout) :: vh
1313 real,
dimension(SZI_(G)),
intent(inout) :: dvhdv
1315 real,
intent(in) :: dt
1317 integer,
intent(in) :: j
1318 integer,
intent(in) :: ish
1319 integer,
intent(in) :: ieh
1320 logical,
dimension(SZI_(G)),
intent(in) :: do_I
1321 logical,
intent(in) :: vol_CFL
1330 logical :: local_open_BC
1332 local_open_bc = .false.
1333 if (
present(obc))
then ;
if (
associated(obc))
then
1334 local_open_bc = obc%open_u_BCs_exist_globally
1337 do i=ish,ieh ;
if (do_i(i))
then
1338 if (v(i) > 0.0)
then
1339 if (vol_cfl)
then ; cfl = (v(i) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j))
1340 else ; cfl = v(i) * dt * g%IdyT(i,j) ;
endif
1341 curv_3 = h_l(i,j) + h_r(i,j) - 2.0*h(i,j)
1342 vh(i) = g%dx_Cv(i,j) * v(i) * ( h_r(i,j) + cfl * &
1343 (0.5*(h_l(i,j) - h_r(i,j)) + curv_3*(cfl - 1.5)) )
1344 h_marg = h_r(i,j) + cfl * ((h_l(i,j) - h_r(i,j)) + &
1345 3.0*curv_3*(cfl - 1.0))
1346 elseif (v(i) < 0.0)
then
1347 if (vol_cfl)
then ; cfl = (-v(i) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j+1))
1348 else ; cfl = -v(i) * dt * g%IdyT(i,j+1) ;
endif
1349 curv_3 = h_l(i,j+1) + h_r(i,j+1) - 2.0*h(i,j+1)
1350 vh(i) = g%dx_Cv(i,j) * v(i) * ( h_l(i,j+1) + cfl * &
1351 (0.5*(h_r(i,j+1)-h_l(i,j+1)) + curv_3*(cfl - 1.5)) )
1352 h_marg = h_l(i,j+1) + cfl * ((h_r(i,j+1)-h_l(i,j+1)) + &
1353 3.0*curv_3*(cfl - 1.0))
1356 h_marg = 0.5 * (h_l(i,j+1) + h_r(i,j))
1358 dvhdv(i) = g%dx_Cv(i,j) * h_marg * visc_rem(i)
1361 if (local_open_bc)
then
1362 do i=ish,ieh ;
if (do_i(i))
then
1363 if (obc%segment(obc%segnum_v(i,j))%open)
then
1365 vh(i) = g%dx_Cv(i,j) * v(i) * h(i,j)
1366 dvhdv(i) = g%dx_Cv(i,j) * h(i,j) * visc_rem(i)
1368 vh(i) = g%dx_Cv(i,j) * v(i) * h(i,j+1)
1369 dvhdv(i) = g%dx_Cv(i,j) * h(i,j+1) * visc_rem(i)
1377 subroutine merid_face_thickness(v, h, h_L, h_R, h_v, dt, G, US, LB, vol_CFL, &
1378 marginal, visc_rem_v, OBC)
1380 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(in) :: v
1381 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
1383 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h_L
1385 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h_R
1387 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(inout) :: h_v
1389 real,
intent(in) :: dt
1392 logical,
intent(in) :: vol_CFL
1394 logical,
intent(in) :: marginal
1396 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
optional,
intent(in) :: visc_rem_v
1409 logical :: local_open_BC
1410 integer :: i, j, k, ish, ieh, jsh, jeh, n, nz
1411 ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh ; nz = g%ke
1414 do k=1,nz ;
do j=jsh-1,jeh ;
do i=ish,ieh
1415 if (v(i,j,k) > 0.0)
then
1416 if (vol_cfl)
then ; cfl = (v(i,j,k) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j))
1417 else ; cfl = v(i,j,k) * dt * g%IdyT(i,j) ;
endif
1418 curv_3 = h_l(i,j,k) + h_r(i,j,k) - 2.0*h(i,j,k)
1419 h_avg = h_r(i,j,k) + cfl * (0.5*(h_l(i,j,k) - h_r(i,j,k)) + curv_3*(cfl - 1.5))
1420 h_marg = h_r(i,j,k) + cfl * ((h_l(i,j,k) - h_r(i,j,k)) + &
1421 3.0*curv_3*(cfl - 1.0))
1422 elseif (v(i,j,k) < 0.0)
then
1423 if (vol_cfl)
then ; cfl = (-v(i,j,k)*dt) * (g%dx_Cv(i,j) * g%IareaT(i,j+1))
1424 else ; cfl = -v(i,j,k) * dt * g%IdyT(i,j+1) ;
endif
1425 curv_3 = h_l(i,j+1,k) + h_r(i,j+1,k) - 2.0*h(i,j+1,k)
1426 h_avg = h_l(i,j+1,k) + cfl * (0.5*(h_r(i,j+1,k)-h_l(i,j+1,k)) + curv_3*(cfl - 1.5))
1427 h_marg = h_l(i,j+1,k) + cfl * ((h_r(i,j+1,k)-h_l(i,j+1,k)) + &
1428 3.0*curv_3*(cfl - 1.0))
1430 h_avg = 0.5 * (h_l(i,j+1,k) + h_r(i,j,k))
1433 h_marg = 0.5 * (h_l(i,j+1,k) + h_r(i,j,k))
1438 if (marginal)
then ; h_v(i,j,k) = h_marg
1439 else ; h_v(i,j,k) = h_avg ;
endif
1440 enddo ;
enddo ;
enddo
1442 if (
present(visc_rem_v))
then
1444 do k=1,nz ;
do j=jsh-1,jeh ;
do i=ish,ieh
1445 h_v(i,j,k) = h_v(i,j,k) * visc_rem_v(i,j,k)
1446 enddo ;
enddo ;
enddo
1449 local_open_bc = .false.
1450 if (
present(obc))
then ;
if (
associated(obc))
then
1451 local_open_bc = obc%open_u_BCs_exist_globally
1453 if (local_open_bc)
then
1454 do n = 1, obc%number_of_segments
1455 if (obc%segment(n)%open .and. obc%segment(n)%is_N_or_S)
then
1456 j = obc%segment(n)%HI%JsdB
1458 if (
present(visc_rem_v))
then ;
do k=1,nz
1459 do i = obc%segment(n)%HI%isd, obc%segment(n)%HI%ied
1460 h_v(i,j,k) = h(i,j,k) * visc_rem_v(i,j,k)
1462 enddo ;
else ;
do k=1,nz
1463 do i = obc%segment(n)%HI%isd, obc%segment(n)%HI%ied
1464 h_v(i,j,k) = h(i,j,k)
1468 if (
present(visc_rem_v))
then ;
do k=1,nz
1469 do i = obc%segment(n)%HI%isd, obc%segment(n)%HI%ied
1470 h_v(i,j,k) = h(i,j+1,k) * visc_rem_v(i,j,k)
1472 enddo ;
else ;
do k=1,nz
1473 do i = obc%segment(n)%HI%isd, obc%segment(n)%HI%ied
1474 h_v(i,j,k) = h(i,j+1,k)
1486 dv, dv_max_CFL, dv_min_CFL, dt, G, US, CS, visc_rem, &
1487 j, ish, ieh, do_I_in, full_precision, vh_3d, OBC)
1489 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1491 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1493 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),&
1495 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1497 real,
dimension(SZI_(G),SZK_(G)),
intent(in) :: visc_rem
1503 real,
dimension(SZI_(G)), &
1504 optional,
intent(in) :: vhbt
1506 real,
dimension(SZI_(G)),
intent(in) :: dv_max_CFL
1507 real,
dimension(SZI_(G)),
intent(in) :: dv_min_CFL
1508 real,
dimension(SZI_(G)),
intent(in) :: vh_tot_0
1510 real,
dimension(SZI_(G)),
intent(in) :: dvhdv_tot_0
1512 real,
dimension(SZI_(G)),
intent(out) :: dv
1513 real,
intent(in) :: dt
1516 integer,
intent(in) :: j
1517 integer,
intent(in) :: ish
1518 integer,
intent(in) :: ieh
1519 logical,
dimension(SZI_(G)), &
1520 intent(in) :: do_I_in
1521 logical,
optional,
intent(in) :: full_precision
1523 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1524 optional,
intent(inout) :: vh_3d
1528 real,
dimension(SZI_(G),SZK_(G)) :: &
1531 real,
dimension(SZI_(G)) :: &
1542 integer :: i, k, nz, itt, max_itts = 20
1543 logical :: full_prec, domore, do_I(SZI_(G))
1546 full_prec = .true. ;
if (
present(full_precision)) full_prec = full_precision
1548 vh_aux(:,:) = 0.0 ; dvhdv(:,:) = 0.0
1550 if (
present(vh_3d))
then ;
do k=1,nz ;
do i=ish,ieh
1551 vh_aux(i,k) = vh_3d(i,j,k)
1552 enddo ;
enddo ;
endif
1555 dv(i) = 0.0 ; do_i(i) = do_i_in(i)
1556 dv_max(i) = dv_max_cfl(i) ; dv_min(i) = dv_min_cfl(i)
1557 vh_err(i) = vh_tot_0(i) - vhbt(i) ; dvhdv_tot(i) = dvhdv_tot_0(i)
1558 vh_err_best(i) = abs(vh_err(i))
1564 case (:1) ; tol_eta = 1e-6 * cs%tol_eta
1565 case (2) ; tol_eta = 1e-4 * cs%tol_eta
1566 case (3) ; tol_eta = 1e-2 * cs%tol_eta
1567 case default ; tol_eta = cs%tol_eta
1570 tol_eta = cs%tol_eta_aux ;
if (itt<=1) tol_eta = 1e-6 * cs%tol_eta_aux
1572 tol_vel = cs%tol_vel
1575 if (vh_err(i) > 0.0)
then ; dv_max(i) = dv(i)
1576 elseif (vh_err(i) < 0.0)
then ; dv_min(i) = dv(i)
1577 else ; do_i(i) = .false. ;
endif
1580 do i=ish,ieh ;
if (do_i(i))
then
1581 if ((dt * min(g%IareaT(i,j),g%IareaT(i,j+1))*abs(vh_err(i)) > tol_eta) .or. &
1582 (cs%better_iter .and. ((abs(vh_err(i)) > tol_vel * dvhdv_tot(i)) .or. &
1583 (abs(vh_err(i)) > vh_err_best(i))) ))
then
1587 ddv = -vh_err(i) / dvhdv_tot(i)
1590 if (abs(ddv) < 1.0e-15*abs(dv(i)))
then
1592 elseif (ddv > 0.0)
then
1593 if (dv(i) >= dv_max(i))
then
1594 dv(i) = 0.5*(dv_prev + dv_max(i))
1595 if (dv_max(i) - dv_prev < 1.0e-15*abs(dv(i))) do_i(i) = .false.
1598 if (dv(i) <= dv_min(i))
then
1599 dv(i) = 0.5*(dv_prev + dv_min(i))
1600 if (dv_prev - dv_min(i) < 1.0e-15*abs(dv(i))) do_i(i) = .false.
1605 dv(i) = dv(i) - vh_err(i) / dvhdv_tot(i)
1606 if ((dv(i) >= dv_max(i)) .or. (dv(i) <= dv_min(i))) &
1607 dv(i) = 0.5*(dv_max(i) + dv_min(i))
1609 if (do_i(i)) domore = .true.
1614 if (.not.domore)
exit
1616 if ((itt < max_itts) .or.
present(vh_3d))
then ;
do k=1,nz
1617 do i=ish,ieh ; v_new(i) = v(i,j,k) + dv(i) * visc_rem(i,k) ;
enddo
1618 call merid_flux_layer(v_new, h_in(:,:,k), h_l(:,:,k), h_r(:,:,k), &
1619 vh_aux(:,k), dvhdv(:,k), visc_rem(:,k), &
1620 dt, g, us, j, ish, ieh, do_i, cs%vol_CFL, obc)
1623 if (itt < max_itts)
then
1625 vh_err(i) = -vhbt(i) ; dvhdv_tot(i) = 0.0
1627 do k=1,nz ;
do i=ish,ieh
1628 vh_err(i) = vh_err(i) + vh_aux(i,k)
1629 dvhdv_tot(i) = dvhdv_tot(i) + dvhdv(i,k)
1632 vh_err_best(i) = min(vh_err_best(i), abs(vh_err(i)))
1640 if (
present(vh_3d))
then ;
do k=1,nz ;
do i=ish,ieh
1641 vh_3d(i,j,k) = vh_aux(i,k)
1642 enddo ;
enddo ;
endif
1648 subroutine set_merid_bt_cont(v, h_in, h_L, h_R, BT_cont, vh_tot_0, dvhdv_tot_0, &
1649 dv_max_CFL, dv_min_CFL, dt, G, US, CS, visc_rem, &
1650 visc_rem_max, j, ish, ieh, do_I)
1652 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(in) :: v
1653 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h_in
1655 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h_L
1657 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h_R
1661 real,
dimension(SZI_(G)),
intent(in) :: vh_tot_0
1663 real,
dimension(SZI_(G)),
intent(in) :: dvhdv_tot_0
1665 real,
dimension(SZI_(G)),
intent(in) :: dv_max_CFL
1667 real,
dimension(SZI_(G)),
intent(in) :: dv_min_CFL
1669 real,
intent(in) :: dt
1672 real,
dimension(SZI_(G),SZK_(G)),
intent(in) :: visc_rem
1677 real,
dimension(SZI_(G)),
intent(in) :: visc_rem_max
1678 integer,
intent(in) :: j
1679 integer,
intent(in) :: ish
1680 integer,
intent(in) :: ieh
1681 logical,
dimension(SZI_(G)),
intent(in) :: do_I
1684 real,
dimension(SZI_(G)) :: &
1704 real :: visc_rem_lim
1707 real :: min_visc_rem
1717 nz = g%ke ; idt = 1.0 / dt
1718 min_visc_rem = 0.1 ; cfl_min = 1e-6
1721 do i=ish,ieh ; zeros(i) = 0.0 ;
enddo
1723 dv_max_cfl, dv_min_cfl, dt, g, us, cs, visc_rem, &
1724 j, ish, ieh, do_i, .true.)
1730 do i=ish,ieh ;
if (do_i(i))
then
1732 dv_cfl(i) = (cfl_min * idt) * g%dyCv(i,j)
1733 dvr(i) = min(0.0,dv0(i) - dv_cfl(i))
1734 dvl(i) = max(0.0,dv0(i) + dv_cfl(i))
1735 famt_l(i) = 0.0 ; famt_r(i) = 0.0 ; famt_0(i) = 0.0
1736 vhtot_l(i) = 0.0 ; vhtot_r(i) = 0.0
1739 if (.not.domore)
then
1740 do k=1,nz ;
do i=ish,ieh
1741 bt_cont%FA_v_S0(i,j) = 0.0 ; bt_cont%FA_v_SS(i,j) = 0.0
1742 bt_cont%vBT_SS(i,j) = 0.0
1743 bt_cont%FA_v_N0(i,j) = 0.0 ; bt_cont%FA_v_NN(i,j) = 0.0
1744 bt_cont%vBT_NN(i,j) = 0.0
1749 do k=1,nz ;
do i=ish,ieh ;
if (do_i(i))
then
1750 visc_rem_lim = max(visc_rem(i,k), min_visc_rem*visc_rem_max(i))
1751 if (visc_rem_lim > 0.0)
then
1752 if (v(i,j,k) + dvr(i)*visc_rem_lim > -dv_cfl(i)*visc_rem(i,k)) &
1753 dvr(i) = -(v(i,j,k) + dv_cfl(i)*visc_rem(i,k)) / visc_rem_lim
1754 if (v(i,j,k) + dvl(i)*visc_rem_lim < dv_cfl(i)*visc_rem(i,k)) &
1755 dvl(i) = -(v(i,j,k) - dv_cfl(i)*visc_rem(i,k)) / visc_rem_lim
1757 endif ;
enddo ;
enddo
1759 do i=ish,ieh ;
if (do_i(i))
then
1760 v_l(i) = v(i,j,k) + dvl(i) * visc_rem(i,k)
1761 v_r(i) = v(i,j,k) + dvr(i) * visc_rem(i,k)
1762 v_0(i) = v(i,j,k) + dv0(i) * visc_rem(i,k)
1764 call merid_flux_layer(v_0, h_in(:,:,k), h_l(:,:,k), h_r(:,:,k), vh_0, dvhdv_0, &
1765 visc_rem(:,k), dt, g, us, j, ish, ieh, do_i, cs%vol_CFL)
1766 call merid_flux_layer(v_l, h_in(:,:,k), h_l(:,:,k), h_r(:,:,k), vh_l, dvhdv_l, &
1767 visc_rem(:,k), dt, g, us, j, ish, ieh, do_i, cs%vol_CFL)
1768 call merid_flux_layer(v_r, h_in(:,:,k), h_l(:,:,k), h_r(:,:,k), vh_r, dvhdv_r, &
1769 visc_rem(:,k), dt, g, us, j, ish, ieh, do_i, cs%vol_CFL)
1770 do i=ish,ieh ;
if (do_i(i))
then
1771 famt_0(i) = famt_0(i) + dvhdv_0(i)
1772 famt_l(i) = famt_l(i) + dvhdv_l(i)
1773 famt_r(i) = famt_r(i) + dvhdv_r(i)
1774 vhtot_l(i) = vhtot_l(i) + vh_l(i)
1775 vhtot_r(i) = vhtot_r(i) + vh_r(i)
1778 do i=ish,ieh ;
if (do_i(i))
then
1779 fa_0 = famt_0(i) ; fa_avg = famt_0(i)
1780 if ((dvl(i) - dv0(i)) /= 0.0) &
1781 fa_avg = vhtot_l(i) / (dvl(i) - dv0(i))
1782 if (fa_avg > max(fa_0, famt_l(i)))
then ; fa_avg = max(fa_0, famt_l(i))
1783 elseif (fa_avg < min(fa_0, famt_l(i)))
then ; fa_0 = fa_avg ;
endif
1784 bt_cont%FA_v_S0(i,j) = fa_0 ; bt_cont%FA_v_SS(i,j) = famt_l(i)
1785 if (abs(fa_0-famt_l(i)) <= 1e-12*fa_0)
then ; bt_cont%vBT_SS(i,j) = 0.0 ;
else
1786 bt_cont%vBT_SS(i,j) = (1.5 * (dvl(i) - dv0(i))) * &
1787 ((famt_l(i) - fa_avg) / (famt_l(i) - fa_0))
1790 fa_0 = famt_0(i) ; fa_avg = famt_0(i)
1791 if ((dvr(i) - dv0(i)) /= 0.0) &
1792 fa_avg = vhtot_r(i) / (dvr(i) - dv0(i))
1793 if (fa_avg > max(fa_0, famt_r(i)))
then ; fa_avg = max(fa_0, famt_r(i))
1794 elseif (fa_avg < min(fa_0, famt_r(i)))
then ; fa_0 = fa_avg ;
endif
1795 bt_cont%FA_v_N0(i,j) = fa_0 ; bt_cont%FA_v_NN(i,j) = famt_r(i)
1796 if (abs(famt_r(i) - fa_0) <= 1e-12*fa_0)
then ; bt_cont%vBT_NN(i,j) = 0.0 ;
else
1797 bt_cont%vBT_NN(i,j) = (1.5 * (dvr(i) - dv0(i))) * &
1798 ((famt_r(i) - fa_avg) / (famt_r(i) - fa_0))
1801 bt_cont%FA_v_S0(i,j) = 0.0 ; bt_cont%FA_v_SS(i,j) = 0.0
1802 bt_cont%FA_v_N0(i,j) = 0.0 ; bt_cont%FA_v_NN(i,j) = 0.0
1803 bt_cont%vBT_SS(i,j) = 0.0 ; bt_cont%vBT_NN(i,j) = 0.0
1811 real,
dimension(SZI_(G),SZJ_(G)),
intent(in) :: h_in
1812 real,
dimension(SZI_(G),SZJ_(G)),
intent(out) :: h_L
1814 real,
dimension(SZI_(G),SZJ_(G)),
intent(out) :: h_R
1817 real,
intent(in) :: h_min
1819 logical,
optional,
intent(in) :: monotonic
1822 logical,
optional,
intent(in) :: simple_2nd
1828 real,
dimension(SZI_(G),SZJ_(G)) :: slp
1829 real,
parameter :: oneSixth = 1./6.
1830 real :: h_ip1, h_im1
1832 logical :: use_CW84, use_2nd
1833 character(len=256) :: mesg
1834 integer :: i, j, isl, iel, jsl, jel, n, stencil
1835 logical :: local_open_BC
1838 use_cw84 = .false. ;
if (
present(monotonic)) use_cw84 = monotonic
1839 use_2nd = .false. ;
if (
present(simple_2nd)) use_2nd = simple_2nd
1841 local_open_bc = .false.
1842 if (
present(obc))
then ;
if (
associated(obc))
then
1843 local_open_bc = obc%open_u_BCs_exist_globally
1846 isl = lb%ish-1 ; iel = lb%ieh+1 ; jsl = lb%jsh ; jel = lb%jeh
1849 stencil = 2 ;
if (use_2nd) stencil = 1
1851 if ((isl-stencil < g%isd) .or. (iel+stencil > g%ied))
then
1852 write(mesg,
'("In MOM_continuity_PPM, PPM_reconstruction_x called with a ", &
1853 & "x-halo that needs to be increased by ",i2,".")') &
1854 stencil + max(g%isd-isl,iel-g%ied)
1857 if ((jsl < g%jsd) .or. (jel > g%jed))
then
1858 write(mesg,
'("In MOM_continuity_PPM, PPM_reconstruction_x called with a ", &
1859 & "y-halo that needs to be increased by ",i2,".")') &
1860 max(g%jsd-jsl,jel-g%jed)
1865 do j=jsl,jel ;
do i=isl,iel
1866 h_im1 = g%mask2dT(i-1,j) * h_in(i-1,j) + (1.0-g%mask2dT(i-1,j)) * h_in(i,j)
1867 h_ip1 = g%mask2dT(i+1,j) * h_in(i+1,j) + (1.0-g%mask2dT(i+1,j)) * h_in(i,j)
1868 h_l(i,j) = 0.5*( h_im1 + h_in(i,j) )
1869 h_r(i,j) = 0.5*( h_ip1 + h_in(i,j) )
1872 do j=jsl,jel ;
do i=isl-1,iel+1
1873 if ((g%mask2dT(i-1,j) * g%mask2dT(i,j) * g%mask2dT(i+1,j)) == 0.0)
then
1877 slp(i,j) = 0.5 * (h_in(i+1,j) - h_in(i-1,j))
1879 dmx = max(h_in(i+1,j), h_in(i-1,j), h_in(i,j)) - h_in(i,j)
1880 dmn = h_in(i,j) - min(h_in(i+1,j), h_in(i-1,j), h_in(i,j))
1881 slp(i,j) = sign(1.,slp(i,j)) * min(abs(slp(i,j)), 2. * min(dmx, dmn))
1886 if (local_open_bc)
then
1887 do n=1, obc%number_of_segments
1888 segment => obc%segment(n)
1889 if (.not. segment%on_pe) cycle
1893 do j=segment%HI%jsd,segment%HI%jed
1901 do j=jsl,jel ;
do i=isl,iel
1906 h_im1 = g%mask2dT(i-1,j) * h_in(i-1,j) + (1.0-g%mask2dT(i-1,j)) * h_in(i,j)
1907 h_ip1 = g%mask2dT(i+1,j) * h_in(i+1,j) + (1.0-g%mask2dT(i+1,j)) * h_in(i,j)
1909 h_l(i,j) = 0.5*( h_im1 + h_in(i,j) ) + onesixth*( slp(i-1,j) - slp(i,j) )
1910 h_r(i,j) = 0.5*( h_ip1 + h_in(i,j) ) + onesixth*( slp(i,j) - slp(i+1,j) )
1914 if (local_open_bc)
then
1915 do n=1, obc%number_of_segments
1916 segment => obc%segment(n)
1917 if (.not. segment%on_pe) cycle
1920 do j=segment%HI%jsd,segment%HI%jed
1921 h_l(i+1,j) = h_in(i,j)
1922 h_r(i+1,j) = h_in(i,j)
1923 h_l(i,j) = h_in(i,j)
1924 h_r(i,j) = h_in(i,j)
1928 do j=segment%HI%jsd,segment%HI%jed
1929 h_l(i,j) = h_in(i+1,j)
1930 h_r(i,j) = h_in(i+1,j)
1931 h_l(i+1,j) = h_in(i+1,j)
1932 h_r(i+1,j) = h_in(i+1,j)
1941 call ppm_limit_pos(h_in, h_l, h_r, h_min, g, isl, iel, jsl, jel)
1950 real,
dimension(SZI_(G),SZJ_(G)),
intent(in) :: h_in
1951 real,
dimension(SZI_(G),SZJ_(G)),
intent(out) :: h_L
1953 real,
dimension(SZI_(G),SZJ_(G)),
intent(out) :: h_R
1956 real,
intent(in) :: h_min
1958 logical,
optional,
intent(in) :: monotonic
1961 logical,
optional,
intent(in) :: simple_2nd
1967 real,
dimension(SZI_(G),SZJ_(G)) :: slp
1968 real,
parameter :: oneSixth = 1./6.
1969 real :: h_jp1, h_jm1
1971 logical :: use_CW84, use_2nd
1972 character(len=256) :: mesg
1973 integer :: i, j, isl, iel, jsl, jel, n, stencil
1974 logical :: local_open_BC
1977 use_cw84 = .false. ;
if (
present(monotonic)) use_cw84 = monotonic
1978 use_2nd = .false. ;
if (
present(simple_2nd)) use_2nd = simple_2nd
1980 local_open_bc = .false.
1981 if (
present(obc))
then ;
if (
associated(obc))
then
1982 local_open_bc = obc%open_u_BCs_exist_globally
1985 isl = lb%ish ; iel = lb%ieh ; jsl = lb%jsh-1 ; jel = lb%jeh+1
1988 stencil = 2 ;
if (use_2nd) stencil = 1
1990 if ((isl < g%isd) .or. (iel > g%ied))
then
1991 write(mesg,
'("In MOM_continuity_PPM, PPM_reconstruction_y called with a ", &
1992 & "x-halo that needs to be increased by ",i2,".")') &
1993 max(g%isd-isl,iel-g%ied)
1996 if ((jsl-stencil < g%jsd) .or. (jel+stencil > g%jed))
then
1997 write(mesg,
'("In MOM_continuity_PPM, PPM_reconstruction_y called with a ", &
1998 & "y-halo that needs to be increased by ",i2,".")') &
1999 stencil + max(g%jsd-jsl,jel-g%jed)
2004 do j=jsl,jel ;
do i=isl,iel
2005 h_jm1 = g%mask2dT(i,j-1) * h_in(i,j-1) + (1.0-g%mask2dT(i,j-1)) * h_in(i,j)
2006 h_jp1 = g%mask2dT(i,j+1) * h_in(i,j+1) + (1.0-g%mask2dT(i,j+1)) * h_in(i,j)
2007 h_l(i,j) = 0.5*( h_jm1 + h_in(i,j) )
2008 h_r(i,j) = 0.5*( h_jp1 + h_in(i,j) )
2011 do j=jsl-1,jel+1 ;
do i=isl,iel
2012 if ((g%mask2dT(i,j-1) * g%mask2dT(i,j) * g%mask2dT(i,j+1)) == 0.0)
then
2016 slp(i,j) = 0.5 * (h_in(i,j+1) - h_in(i,j-1))
2018 dmx = max(h_in(i,j+1), h_in(i,j-1), h_in(i,j)) - h_in(i,j)
2019 dmn = h_in(i,j) - min(h_in(i,j+1), h_in(i,j-1), h_in(i,j))
2020 slp(i,j) = sign(1.,slp(i,j)) * min(abs(slp(i,j)), 2. * min(dmx, dmn))
2025 if (local_open_bc)
then
2026 do n=1, obc%number_of_segments
2027 segment => obc%segment(n)
2028 if (.not. segment%on_pe) cycle
2032 do i=segment%HI%isd,segment%HI%ied
2040 do j=jsl,jel ;
do i=isl,iel
2043 h_jm1 = g%mask2dT(i,j-1) * h_in(i,j-1) + (1.0-g%mask2dT(i,j-1)) * h_in(i,j)
2044 h_jp1 = g%mask2dT(i,j+1) * h_in(i,j+1) + (1.0-g%mask2dT(i,j+1)) * h_in(i,j)
2046 h_l(i,j) = 0.5*( h_jm1 + h_in(i,j) ) + onesixth*( slp(i,j-1) - slp(i,j) )
2047 h_r(i,j) = 0.5*( h_jp1 + h_in(i,j) ) + onesixth*( slp(i,j) - slp(i,j+1) )
2051 if (local_open_bc)
then
2052 do n=1, obc%number_of_segments
2053 segment => obc%segment(n)
2054 if (.not. segment%on_pe) cycle
2057 do i=segment%HI%isd,segment%HI%ied
2058 h_l(i,j+1) = h_in(i,j)
2059 h_r(i,j+1) = h_in(i,j)
2060 h_l(i,j) = h_in(i,j)
2061 h_r(i,j) = h_in(i,j)
2065 do i=segment%HI%isd,segment%HI%ied
2066 h_l(i,j) = h_in(i,j+1)
2067 h_r(i,j) = h_in(i,j+1)
2068 h_l(i,j+1) = h_in(i,j+1)
2069 h_r(i,j+1) = h_in(i,j+1)
2078 call ppm_limit_pos(h_in, h_l, h_r, h_min, g, isl, iel, jsl, jel)
2088 subroutine ppm_limit_pos(h_in, h_L, h_R, h_min, G, iis, iie, jis, jie)
2090 real,
dimension(SZI_(G),SZJ_(G)),
intent(in) :: h_in
2091 real,
dimension(SZI_(G),SZJ_(G)),
intent(inout) :: h_L
2092 real,
dimension(SZI_(G),SZJ_(G)),
intent(inout) :: h_R
2093 real,
intent(in) :: h_min
2095 integer,
intent(in) :: iis
2096 integer,
intent(in) :: iie
2097 integer,
intent(in) :: jis
2098 integer,
intent(in) :: jie
2101 real :: curv, dh, scale
2102 character(len=256) :: mesg
2105 do j=jis,jie ;
do i=iis,iie
2108 curv = 3.0*(h_l(i,j) + h_r(i,j) - 2.0*h_in(i,j))
2109 if (curv > 0.0)
then
2110 dh = h_r(i,j) - h_l(i,j)
2111 if (abs(dh) < curv)
then
2112 if (h_in(i,j) <= h_min)
then
2113 h_l(i,j) = h_in(i,j) ; h_r(i,j) = h_in(i,j)
2114 elseif (12.0*curv*(h_in(i,j) - h_min) < (curv**2 + 3.0*dh**2))
then
2117 scale = 12.0*curv*(h_in(i,j) - h_min) / (curv**2 + 3.0*dh**2)
2118 h_l(i,j) = h_in(i,j) + scale*(h_l(i,j) - h_in(i,j))
2119 h_r(i,j) = h_in(i,j) + scale*(h_r(i,j) - h_in(i,j))
2131 real,
dimension(SZI_(G),SZJ_(G)),
intent(in) :: h_in
2132 real,
dimension(SZI_(G),SZJ_(G)),
intent(inout) :: h_L
2134 real,
dimension(SZI_(G),SZJ_(G)),
intent(inout) :: h_R
2136 integer,
intent(in) :: iis
2137 integer,
intent(in) :: iie
2138 integer,
intent(in) :: jis
2139 integer,
intent(in) :: jie
2142 real :: h_i, RLdiff, RLdiff2, RLmean, FunFac
2143 character(len=256) :: mesg
2146 do j=jis,jie ;
do i=iis,iie
2150 if ( ( h_r(i,j) - h_i ) * ( h_i - h_l(i,j) ) <= 0. )
then
2151 h_l(i,j) = h_i ; h_r(i,j) = h_i
2153 rldiff = h_r(i,j) - h_l(i,j)
2154 rlmean = 0.5 * ( h_r(i,j) + h_l(i,j) )
2155 funfac = 6. * rldiff * ( h_i - rlmean )
2156 rldiff2 = rldiff * rldiff
2157 if ( funfac > rldiff2 ) h_l(i,j) = 3. * h_i - 2. * h_r(i,j)
2158 if ( funfac < -rldiff2 ) h_r(i,j) = 3. * h_i - 2. * h_l(i,j)
2166 function ratio_max(a, b, maxrat)
result(ratio)
2167 real,
intent(in) :: a
2168 real,
intent(in) :: b
2169 real,
intent(in) :: maxrat
2172 if (abs(a) > abs(maxrat*b))
then
2181 type(time_type),
target,
intent(in) :: time
2187 type(
diag_ctrl),
target,
intent(inout) :: diag
2191 #include "version_variable.h"
2193 character(len=40) :: mdl =
"MOM_continuity_PPM"
2195 if (
associated(cs))
then
2196 call mom_error(warning,
"continuity_PPM_init called with associated control structure.")
2203 call get_param(param_file, mdl,
"MONOTONIC_CONTINUITY", cs%monotonic, &
2204 "If true, CONTINUITY_PPM uses the Colella and Woodward "//&
2205 "monotonic limiter. The default (false) is to use a "//&
2206 "simple positive definite limiter.", default=.false.)
2207 call get_param(param_file, mdl,
"SIMPLE_2ND_PPM_CONTINUITY", cs%simple_2nd, &
2208 "If true, CONTINUITY_PPM uses a simple 2nd order "//&
2209 "(arithmetic mean) interpolation of the edge values. "//&
2210 "This may give better PV conservation properties. While "//&
2211 "it formally reduces the accuracy of the continuity "//&
2212 "solver itself in the strongly advective limit, it does "//&
2213 "not reduce the overall order of accuracy of the dynamic "//&
2214 "core.", default=.false.)
2215 call get_param(param_file, mdl,
"UPWIND_1ST_CONTINUITY", cs%upwind_1st, &
2216 "If true, CONTINUITY_PPM becomes a 1st-order upwind "//&
2217 "continuity solver. This scheme is highly diffusive "//&
2218 "but may be useful for debugging or in single-column "//&
2219 "mode where its minimal stencil is useful.", default=.false.)
2220 call get_param(param_file, mdl,
"ETA_TOLERANCE", cs%tol_eta, &
2221 "The tolerance for the differences between the "//&
2222 "barotropic and baroclinic estimates of the sea surface "//&
2223 "height due to the fluxes through each face. The total "//&
2224 "tolerance for SSH is 4 times this value. The default "//&
2225 "is 0.5*NK*ANGSTROM, and this should not be set less "//&
2226 "than about 10^-15*MAXIMUM_DEPTH.", units=
"m", scale=gv%m_to_H, &
2227 default=0.5*g%ke*gv%Angstrom_m, unscaled=tol_eta_m)
2230 call get_param(param_file, mdl,
"ETA_TOLERANCE_AUX", cs%tol_eta_aux, &
2231 "The tolerance for free-surface height discrepancies "//&
2232 "between the barotropic solution and the sum of the "//&
2233 "layer thicknesses when calculating the auxiliary "//&
2234 "corrected velocities. By default, this is the same as "//&
2235 "ETA_TOLERANCE, but can be made larger for efficiency.", &
2236 units=
"m", default=tol_eta_m, scale=gv%m_to_H)
2237 call get_param(param_file, mdl,
"VELOCITY_TOLERANCE", cs%tol_vel, &
2238 "The tolerance for barotropic velocity discrepancies "//&
2239 "between the barotropic solution and the sum of the "//&
2240 "layer thicknesses.", units=
"m s-1", default=3.0e8, scale=us%m_s_to_L_T)
2243 call get_param(param_file, mdl,
"CONT_PPM_AGGRESS_ADJUST", cs%aggress_adjust,&
2244 "If true, allow the adjusted velocities to have a "//&
2245 "relative CFL change up to 0.5.", default=.false.)
2246 cs%vol_CFL = cs%aggress_adjust
2247 call get_param(param_file, mdl,
"CONT_PPM_VOLUME_BASED_CFL", cs%vol_CFL, &
2248 "If true, use the ratio of the open face lengths to the "//&
2249 "tracer cell areas when estimating CFL numbers. The "//&
2250 "default is set by CONT_PPM_AGGRESS_ADJUST.", &
2251 default=cs%aggress_adjust, do_not_read=cs%aggress_adjust)
2252 call get_param(param_file, mdl,
"CONTINUITY_CFL_LIMIT", cs%CFL_limit_adjust, &
2253 "The maximum CFL of the adjusted velocities.", units=
"nondim", &
2255 call get_param(param_file, mdl,
"CONT_PPM_BETTER_ITER", cs%better_iter, &
2256 "If true, stop corrective iterations using a velocity "//&
2257 "based criterion and only stop if the iteration is "//&
2258 "better than all predecessors.", default=.true.)
2259 call get_param(param_file, mdl,
"CONT_PPM_USE_VISC_REM_MAX", &
2260 cs%use_visc_rem_max, &
2261 "If true, use more appropriate limiting bounds for "//&
2262 "corrections in strongly viscous columns.", default=.true.)
2263 call get_param(param_file, mdl,
"CONT_PPM_MARGINAL_FACE_AREAS", cs%marginal_faces, &
2264 "If true, use the marginal face areas from the continuity "//&
2265 "solver for use as the weights in the barotropic solver. "//&
2266 "Otherwise use the transport averaged areas.", default=.true.)
2270 id_clock_update = cpu_clock_id(
'(Ocean continuity update)', grain=clock_routine)
2271 id_clock_correct = cpu_clock_id(
'(Ocean continuity correction)', grain=clock_routine)
2280 stencil = 3 ;
if (cs%simple_2nd) stencil = 2 ;
if (cs%upwind_1st) stencil = 1