6 use mom_cpu_clock,
only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
10 use mom_domains,
only : sum_across_pes, max_across_pes
21 implicit none ;
private
23 #include <MOM_memory.h>
37 type(group_pass_type) :: pass_uhr_vhr_t_hprev
50 subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, &
51 h_prev_opt, max_iter_in, x_first_in, uhr_out, vhr_out, h_out)
54 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
56 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
58 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
61 real,
intent(in) :: dt
65 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
66 optional,
intent(in) :: h_prev_opt
67 integer,
optional,
intent(in) :: max_iter_in
68 logical,
optional,
intent(in) :: x_first_in
70 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
71 optional,
intent(out) :: uhr_out
73 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
74 optional,
intent(out) :: vhr_out
76 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
77 optional,
intent(out) :: h_out
80 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
82 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: &
84 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)) :: &
86 real :: uh_neglect(szib_(g),szj_(g))
87 real :: vh_neglect(szi_(g),szjb_(g))
92 logical :: domore_u(szj_(g),szk_(g))
93 logical :: domore_v(szjb_(g),szk_(g))
97 integer :: domore_k(szk_(g))
100 integer :: i, j, k, m, is, ie, js, je, isd, ied, jsd, jed, nz, itt, ntr, do_any
101 integer :: isv, iev, jsv, jev
102 integer :: isdb, iedb, jsdb, jedb
104 domore_u(:,:) = .false.
105 domore_v(:,:) = .false.
106 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
107 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
108 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
109 landvolfill = 1.0e-20
112 if (.not.
associated(cs))
call mom_error(fatal,
"MOM_tracer_advect: "// &
113 "tracer_advect_init must be called before advect_tracer.")
114 if (.not.
associated(reg))
call mom_error(fatal,
"MOM_tracer_advect: "// &
115 "register_tracer must be called before advect_tracer.")
116 if (reg%ntr==0)
return
118 x_first = (mod(g%first_direction,2) == 0)
121 if (cs%usePPM .and. .not. cs%useHuynh) stencil = 3
124 do m=1,ntr ; tr(m) = reg%Tr(m) ;
enddo
127 max_iter = 2*int(ceiling(dt/cs%dt)) + 1
129 if (
present(max_iter_in)) max_iter = max_iter_in
130 if (
present(x_first_in)) x_first = x_first_in
148 do j=jsd,jed ;
do i=isdb,iedb ; uhr(i,j,k) = 0.0 ;
enddo ;
enddo
149 do j=jsdb,jedb ;
do i=isd,ied ; vhr(i,j,k) = 0.0 ;
enddo ;
enddo
150 do j=jsd,jed ;
do i=isd,ied ; hprev(i,j,k) = 0.0 ;
enddo ;
enddo
153 do j=js,je ;
do i=is-1,ie ; uhr(i,j,k) = uhtr(i,j,k) ;
enddo ;
enddo
154 do j=js-1,je ;
do i=is,ie ; vhr(i,j,k) = vhtr(i,j,k) ;
enddo ;
enddo
155 if (.not.
present(h_prev_opt))
then
159 do i=is,ie ;
do j=js,je
160 hprev(i,j,k) = max(0.0, g%areaT(i,j)*h_end(i,j,k) + &
161 ((uhr(i,j,k) - uhr(i-1,j,k)) + (vhr(i,j,k) - vhr(i,j-1,k))))
165 hprev(i,j,k) = hprev(i,j,k) + &
166 max(0.0, 1.0e-13*hprev(i,j,k) - g%areaT(i,j)*h_end(i,j,k))
169 do i=is,ie ;
do j=js,je
170 hprev(i,j,k) = h_prev_opt(i,j,k)
177 do j=jsd,jed ;
do i=isd,ied-1
178 uh_neglect(i,j) = gv%H_subroundoff*min(g%areaT(i,j),g%areaT(i+1,j))
181 do j=jsd,jed-1 ;
do i=isd,ied
182 vh_neglect(i,j) = gv%H_subroundoff*min(g%areaT(i,j),g%areaT(i,j+1))
188 if (
associated(tr(m)%ad_x))
then
189 do k=1,nz ;
do j=jsd,jed ;
do i=isd,ied
190 tr(m)%ad_x(i,j,k) = 0.0
191 enddo ;
enddo ;
enddo
193 if (
associated(tr(m)%ad_y))
then
194 do k=1,nz ;
do j=jsd,jed ;
do i=isd,ied
195 tr(m)%ad_y(i,j,k) = 0.0
196 enddo ;
enddo ;
enddo
198 if (
associated(tr(m)%advection_xy))
then
199 do k=1,nz ;
do j=jsd,jed ;
do i=isd,ied
200 tr(m)%advection_xy(i,j,k) = 0.0
201 enddo ;
enddo ;
enddo
203 if (
associated(tr(m)%ad2d_x))
then
204 do j=jsd,jed ;
do i=isd,ied ; tr(m)%ad2d_x(i,j) = 0.0 ;
enddo ;
enddo
206 if (
associated(tr(m)%ad2d_y))
then
207 do j=jsd,jed ;
do i=isd,ied ; tr(m)%ad2d_y(i,j) = 0.0 ;
enddo ;
enddo
212 isv = is ; iev = ie ; jsv = js ; jev = je
216 if (isv > is-stencil)
then
219 nsten_halo = min(is-isd,ied-ie,js-jsd,jed-je)/stencil
220 isv = is-nsten_halo*stencil ; jsv = js-nsten_halo*stencil
221 iev = ie+nsten_halo*stencil ; jev = je+nsten_halo*stencil
224 if ((nsten_halo > 1) .or. (itt==1))
then
227 do k=1,nz ;
if (domore_k(k) > 0)
then
228 do j=jsv,jev ;
if (.not.domore_u(j,k))
then
229 do i=isv+stencil-1,iev-stencil;
if (uhr(i,j,k) /= 0.0)
then
230 domore_u(j,k) = .true. ;
exit
233 do j=jsv+stencil-1,jev-stencil ;
if (.not.domore_v(j,k))
then
234 do i=isv+stencil,iev-stencil;
if (vhr(i,j,k) /= 0.0)
then
235 domore_v(j,k) = .true. ;
exit
242 do j=jsv,jev ;
if (domore_u(j,k)) domore_k(k) = 1 ;
enddo
243 do j=jsv+stencil-1,jev-stencil ;
if (domore_v(j,k)) domore_k(k) = 1 ;
enddo
250 isv = isv + stencil ; iev = iev - stencil
251 jsv = jsv + stencil ; jev = jev - stencil
262 do k=1,nz ;
if (domore_k(k) > 0)
then
267 call advect_x(tr, hprev, uhr, uh_neglect, obc, domore_u, ntr, idt, &
268 isv, iev, jsv-stencil, jev+stencil, k, g, gv, us, cs%usePPM, cs%useHuynh)
271 call advect_y(tr, hprev, vhr, vh_neglect, obc, domore_v, ntr, idt, &
272 isv, iev, jsv, jev, k, g, gv, us, cs%usePPM, cs%useHuynh)
275 do j=jsv-stencil,jev+stencil ;
if (domore_u(j,k)) domore_k(k) = 1 ;
enddo
276 do j=jsv-1,jev ;
if (domore_v(j,k)) domore_k(k) = 1 ;
enddo
281 call advect_y(tr, hprev, vhr, vh_neglect, obc, domore_v, ntr, idt, &
282 isv-stencil, iev+stencil, jsv, jev, k, g, gv, us, cs%usePPM, cs%useHuynh)
285 call advect_x(tr, hprev, uhr, uh_neglect, obc, domore_u, ntr, idt, &
286 isv, iev, jsv, jev, k, g, gv, us, cs%usePPM, cs%useHuynh)
289 do j=jsv,jev ;
if (domore_u(j,k)) domore_k(k) = 1 ;
enddo
290 do j=jsv-1,jev ;
if (domore_v(j,k)) domore_k(k) = 1 ;
enddo
298 if (itt >= max_iter)
then
303 if (isv > is-stencil)
then
306 call sum_across_pes(domore_k(:), nz)
308 do k=1,nz ; do_any = do_any + domore_k(k) ;
enddo
309 if (do_any == 0)
then
317 if (
present(uhr_out)) uhr_out(:,:,:) = uhr(:,:,:)
318 if (
present(vhr_out)) vhr_out(:,:,:) = vhr(:,:,:)
319 if (
present(h_out)) h_out(:,:,:) = hprev(:,:,:)
328 subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &
329 is, ie, js, je, k, G, GV, US, usePPM, useHuynh)
332 type(
tracer_type),
dimension(ntr),
intent(inout) :: Tr
333 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(inout) :: hprev
335 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(inout) :: uhr
337 real,
dimension(SZIB_(G),SZJ_(G)),
intent(inout) :: uh_neglect
340 logical,
dimension(SZJ_(G),SZK_(G)),
intent(inout) :: domore_u
342 real,
intent(in) :: Idt
343 integer,
intent(in) :: ntr
344 integer,
intent(in) :: is
345 integer,
intent(in) :: ie
346 integer,
intent(in) :: js
347 integer,
intent(in) :: je
348 integer,
intent(in) :: k
350 logical,
intent(in) :: usePPM
351 logical,
intent(in) :: useHuynh
354 real,
dimension(SZI_(G),ntr) :: &
356 real,
dimension(SZIB_(G),ntr) :: &
358 real,
dimension(SZI_(G),ntr) :: &
367 real :: uhh(SZIB_(G))
369 real,
dimension(SZIB_(G)) :: &
377 logical :: do_i(SZIB_(G))
379 integer :: i, j, m, n, i_up, stencil
380 real :: aR, aL, dMx, dMn, Tp, Tc, Tm, dA, mA, a6
381 real :: fac1,u_L_in,u_L_out
383 logical :: usePLMslope
385 useplmslope = .not. (useppm .and. usehuynh)
388 if (useppm .and. .not. usehuynh) stencil = 2
390 min_h = 0.1*gv%Angstrom_H
391 h_neglect = gv%H_subroundoff
394 do i=is-1,ie ; cfl(i) = 0.0 ;
enddo
396 do j=js,je ;
if (domore_u(j,k))
then
397 domore_u(j,k) = .false.
401 if (useplmslope)
then
402 do m=1,ntr ;
do i=is-stencil,ie+stencil
417 tp = tr(m)%t(i+1,j,k) ; tc = tr(m)%t(i,j,k) ; tm = tr(m)%t(i-1,j,k)
418 dmx = max( tp, tc, tm ) - tc
419 dmn= tc - min( tp, tc, tm )
420 slope_x(i,m) = g%mask2dCu(i,j)*g%mask2dCu(i-1,j) * &
421 sign( min(0.5*abs(tp-tm), 2.0*dmx, 2.0*dmn), tp-tm )
428 t_tmp(i,m) = tr(m)%t(i,j,k)
432 if (
associated(obc))
then ;
if (obc%OBC_pe)
then
433 do n=1,obc%number_of_segments
434 segment=>obc%segment(n)
435 if (.not.
associated(segment%tr_Reg)) cycle
436 if (segment%is_E_or_W)
then
437 if (j>=segment%HI%jsd .and. j<=segment%HI%jed)
then
440 if (
associated(segment%tr_Reg%Tr(m)%tres))
then
441 if (segment%direction == obc_direction_w)
then
442 t_tmp(i,m) = segment%tr_Reg%Tr(m)%tres(i,j,k)
444 t_tmp(i+1,m) = segment%tr_Reg%Tr(m)%tres(i,j,k)
447 if (segment%direction == obc_direction_w)
then
448 t_tmp(i,m) = segment%tr_Reg%Tr(m)%OBC_inflow_conc
450 t_tmp(i+1,m) = segment%tr_Reg%Tr(m)%OBC_inflow_conc
455 do i=segment%HI%IsdB-1,segment%HI%IsdB+1
456 tp = t_tmp(i+1,m) ; tc = t_tmp(i,m) ; tm = t_tmp(i-1,m)
457 dmx = max( tp, tc, tm ) - tc
458 dmn= tc - min( tp, tc, tm )
459 slope_x(i,m) = g%mask2dCu(i,j)*g%mask2dCu(i-1,j) * &
460 sign( min(0.5*abs(tp-tm), 2.0*dmx, 2.0*dmn), tp-tm )
475 if (uhr(i,j,k) == 0.0)
then
478 elseif (uhr(i,j,k) < 0.0)
then
479 hup = hprev(i+1,j,k) - g%areaT(i+1,j)*min_h
480 hlos = max(0.0,uhr(i+1,j,k))
481 if ((((hup - hlos) + uhr(i,j,k)) < 0.0) .and. &
482 ((0.5*hup + uhr(i,j,k)) < 0.0))
then
483 uhh(i) = min(-0.5*hup,-hup+hlos,0.0)
484 domore_u(j,k) = .true.
489 cfl(i) = - uhh(i) / (hprev(i+1,j,k) + h_neglect*g%areaT(i+1,j))
491 hup = hprev(i,j,k) - g%areaT(i,j)*min_h
492 hlos = max(0.0,-uhr(i-1,j,k))
493 if ((((hup - hlos) - uhr(i,j,k)) < 0.0) .and. &
494 ((0.5*hup - uhr(i,j,k)) < 0.0))
then
495 uhh(i) = max(0.5*hup,hup-hlos,0.0)
496 domore_u(j,k) = .true.
501 cfl(i) = uhh(i) / (hprev(i,j,k) + h_neglect*g%areaT(i,j))
507 do m=1,ntr ;
do i=is-1,ie
509 if (uhh(i) >= 0.0)
then
516 tp = t_tmp(i_up+1,m) ; tc = t_tmp(i_up,m) ; tm = t_tmp(i_up-1,m)
519 al = ( 5.*tc + ( 2.*tm - tp ) )/6.
520 al = max( min(tc,tm), al) ; al = min( max(tc,tm), al)
521 ar = ( 5.*tc + ( 2.*tp - tm ) )/6.
522 ar = max( min(tc,tp), ar) ; ar = min( max(tc,tp), ar)
524 al = 0.5 * ((tm + tc) + (slope_x(i_up-1,m) - slope_x(i_up,m)) / 3.)
525 ar = 0.5 * ((tc + tp) + (slope_x(i_up,m) - slope_x(i_up+1,m)) / 3.)
528 da = ar - al ; ma = 0.5*( ar + al )
529 if (g%mask2dCu(i_up,j)*g%mask2dCu(i_up-1,j)*(tp-tc)*(tc-tm) <= 0.)
then
531 elseif ( da*(tc-ma) > (da*da)/6. )
then
533 elseif ( da*(tc-ma) < - (da*da)/6. )
then
537 a6 = 6.*tc - 3. * (ar + al)
539 if (uhh(i) >= 0.0)
then
540 flux_x(i,m) = uhh(i)*( ar - 0.5 * cfl(i) * ( &
541 ( ar - al ) - a6 * ( 1. - 2./3. * cfl(i) ) ) )
543 flux_x(i,m) = uhh(i)*( al + 0.5 * cfl(i) * ( &
544 ( ar - al ) + a6 * ( 1. - 2./3. * cfl(i) ) ) )
548 do m=1,ntr ;
do i=is-1,ie
549 if (uhh(i) >= 0.0)
then
559 flux_x(i,m) = uhh(i)*( tc + 0.5 * slope_x(i,m) * ( 1. - cfl(i) ) )
572 flux_x(i,m) = uhh(i)*( tc - 0.5 * slope_x(i+1,m) * ( 1. - cfl(i) ) )
580 if (
associated(obc))
then ;
if (obc%OBC_pe)
then
581 if (obc%specified_u_BCs_exist_globally .or. obc%open_u_BCs_exist_globally)
then
582 do n=1,obc%number_of_segments
583 segment=>obc%segment(n)
584 if (.not.
associated(segment%tr_Reg)) cycle
585 if (segment%is_E_or_W)
then
586 if (j>=segment%HI%jsd .and. j<=segment%HI%jed)
then
590 if ((uhr(i,j,k) > 0.0) .and. (segment%direction == obc_direction_w) .or. &
591 (uhr(i,j,k) < 0.0) .and. (segment%direction == obc_direction_e))
then
595 if (
associated(segment%tr_Reg%Tr(m)%tres))
then
596 flux_x(i,m) = uhh(i)*segment%tr_Reg%Tr(m)%tres(i,j,k)
597 else ; flux_x(i,m) = uhh(i)*segment%tr_Reg%Tr(m)%OBC_inflow_conc ;
endif
605 if (obc%open_u_BCs_exist_globally)
then
606 do n=1,obc%number_of_segments
607 segment=>obc%segment(n)
609 if (segment%is_E_or_W .and. (j >= segment%HI%jsd .and. j<= segment%HI%jed))
then
610 if (segment%specified) cycle
611 if (.not.
associated(segment%tr_Reg)) cycle
614 if ((uhr(i,j,k) > 0.0) .and. (g%mask2dT(i,j) < 0.5) .or. &
615 (uhr(i,j,k) < 0.0) .and. (g%mask2dT(i+1,j) < 0.5))
then
618 if (
associated(segment%tr_Reg%Tr(m)%tres))
then
619 flux_x(i,m) = uhh(i)*segment%tr_Reg%Tr(m)%tres(i,j,k)
620 else; flux_x(i,m) = uhh(i)*segment%tr_Reg%Tr(m)%OBC_inflow_conc;
endif
631 uhr(i,j,k) = uhr(i,j,k) - uhh(i)
632 if (abs(uhr(i,j,k)) < uh_neglect(i,j)) uhr(i,j,k) = 0.0
635 if ((uhh(i) /= 0.0) .or. (uhh(i-1) /= 0.0))
then
637 hlst(i) = hprev(i,j,k)
638 hprev(i,j,k) = hprev(i,j,k) - (uhh(i) - uhh(i-1))
639 if (hprev(i,j,k) <= 0.0)
then ; do_i(i) = .false.
640 elseif (hprev(i,j,k) < h_neglect*g%areaT(i,j))
then
641 hlst(i) = hlst(i) + (h_neglect*g%areaT(i,j) - hprev(i,j,k))
642 ihnew(i) = 1.0 / (h_neglect*g%areaT(i,j))
643 else ; ihnew(i) = 1.0 / hprev(i,j,k) ;
endif
655 if (ihnew(i) > 0.0)
then
656 tr(m)%t(i,j,k) = (tr(m)%t(i,j,k) * hlst(i) - &
657 (flux_x(i,m) - flux_x(i-1,m))) * ihnew(i)
663 if (
associated(tr(m)%ad_x))
then ;
do i=is,ie ;
if (do_i(i))
then
664 tr(m)%ad_x(i,j,k) = tr(m)%ad_x(i,j,k) + flux_x(i,m)*idt
665 endif ;
enddo ;
endif
666 if (
associated(tr(m)%ad2d_x))
then ;
do i=is,ie ;
if (do_i(i))
then
667 tr(m)%ad2d_x(i,j) = tr(m)%ad2d_x(i,j) + flux_x(i,m)*idt
668 endif ;
enddo ;
endif
672 if (
associated(tr(m)%advection_xy))
then
673 do i=is,ie ;
if (do_i(i))
then
674 tr(m)%advection_xy(i,j,k) = tr(m)%advection_xy(i,j,k) - (flux_x(i,m) - flux_x(i-1,m)) * &
687 subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, &
688 is, ie, js, je, k, G, GV, US, usePPM, useHuynh)
691 type(
tracer_type),
dimension(ntr),
intent(inout) :: Tr
692 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(inout) :: hprev
694 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(inout) :: vhr
696 real,
dimension(SZI_(G),SZJB_(G)),
intent(inout) :: vh_neglect
699 logical,
dimension(SZJB_(G),SZK_(G)),
intent(inout) :: domore_v
701 real,
intent(in) :: Idt
702 integer,
intent(in) :: ntr
703 integer,
intent(in) :: is
704 integer,
intent(in) :: ie
705 integer,
intent(in) :: js
706 integer,
intent(in) :: je
707 integer,
intent(in) :: k
709 logical,
intent(in) :: usePPM
710 logical,
intent(in) :: useHuynh
713 real,
dimension(SZI_(G),ntr,SZJ_(G)) :: &
715 real,
dimension(SZI_(G),ntr,SZJB_(G)) :: &
717 real,
dimension(SZI_(G),ntr,SZJB_(G)) :: &
721 real :: vhh(SZI_(G),SZJB_(G))
727 real,
dimension(SZIB_(G)) :: &
735 logical :: do_j_tr(SZJ_(G))
736 logical :: do_i(SZIB_(G))
738 integer :: i, j, j2, m, n, j_up, stencil
739 real :: aR, aL, dMx, dMn, Tp, Tc, Tm, dA, mA, a6
740 real :: fac1,v_L_in,v_L_out
742 logical :: usePLMslope
744 useplmslope = .not. (useppm .and. usehuynh)
747 if (useppm .and. .not. usehuynh) stencil = 2
749 min_h = 0.1*gv%Angstrom_H
750 h_neglect = gv%H_subroundoff
764 do j=js-1,je ;
if (domore_v(j,k))
then ;
do j2=1-stencil,stencil ; do_j_tr(j+j2) = .true. ;
enddo ;
endif ;
enddo
768 if (useplmslope)
then
769 do j=js-stencil,je+stencil ;
if (do_j_tr(j))
then ;
do m=1,ntr ;
do i=is,ie
784 tp = tr(m)%t(i,j+1,k) ; tc = tr(m)%t(i,j,k) ; tm = tr(m)%t(i,j-1,k)
785 dmx = max( tp, tc, tm ) - tc
786 dmn= tc - min( tp, tc, tm )
787 slope_y(i,m,j) = g%mask2dCv(i,j)*g%mask2dCv(i,j-1) * &
788 sign( min(0.5*abs(tp-tm), 2.0*dmx, 2.0*dmn), tp-tm )
789 enddo ;
enddo ;
endif ;
enddo
795 do j=g%jsd,g%jed;
do m=1,ntr;
do i=g%isd,g%ied
796 t_tmp(i,m,j) = tr(m)%t(i,j,k)
797 enddo ;
enddo ;
enddo
800 if (
associated(obc))
then ;
if (obc%OBC_pe)
then
801 do n=1,obc%number_of_segments
802 segment=>obc%segment(n)
803 if (.not.
associated(segment%tr_Reg)) cycle
805 if (segment%is_N_or_S)
then
806 if (i>=segment%HI%isd .and. i<=segment%HI%ied)
then
809 if (
associated(segment%tr_Reg%Tr(m)%tres))
then
810 if (segment%direction == obc_direction_s)
then
811 t_tmp(i,m,j) = segment%tr_Reg%Tr(m)%tres(i,j,k)
813 t_tmp(i,m,j+1) = segment%tr_Reg%Tr(m)%tres(i,j,k)
816 if (segment%direction == obc_direction_s)
then
817 t_tmp(i,m,j) = segment%tr_Reg%Tr(m)%OBC_inflow_conc
819 t_tmp(i,m,j+1) = segment%tr_Reg%Tr(m)%OBC_inflow_conc
824 do j=segment%HI%JsdB-1,segment%HI%JsdB+1
825 tp = t_tmp(i,m,j+1) ; tc = t_tmp(i,m,j) ; tm = t_tmp(i,m,j-1)
826 dmx = max( tp, tc, tm ) - tc
827 dmn= tc - min( tp, tc, tm )
828 slope_y(i,m,j) = g%mask2dCv(i,j)*g%mask2dCv(i,j-1) * &
829 sign( min(0.5*abs(tp-tm), 2.0*dmx, 2.0*dmn), tp-tm )
842 do j=js-1,je ;
if (domore_v(j,k))
then
843 domore_v(j,k) = .false.
846 if (vhr(i,j,k) == 0.0)
then
849 elseif (vhr(i,j,k) < 0.0)
then
850 hup = hprev(i,j+1,k) - g%areaT(i,j+1)*min_h
851 hlos = max(0.0,vhr(i,j+1,k))
852 if ((((hup - hlos) + vhr(i,j,k)) < 0.0) .and. &
853 ((0.5*hup + vhr(i,j,k)) < 0.0))
then
854 vhh(i,j) = min(-0.5*hup,-hup+hlos,0.0)
855 domore_v(j,k) = .true.
857 vhh(i,j) = vhr(i,j,k)
860 cfl(i) = - vhh(i,j) / (hprev(i,j+1,k) + h_neglect*g%areaT(i,j+1))
862 hup = hprev(i,j,k) - g%areaT(i,j)*min_h
863 hlos = max(0.0,-vhr(i,j-1,k))
864 if ((((hup - hlos) - vhr(i,j,k)) < 0.0) .and. &
865 ((0.5*hup - vhr(i,j,k)) < 0.0))
then
866 vhh(i,j) = max(0.5*hup,hup-hlos,0.0)
867 domore_v(j,k) = .true.
869 vhh(i,j) = vhr(i,j,k)
872 cfl(i) = vhh(i,j) / (hprev(i,j,k) + h_neglect*g%areaT(i,j))
877 do m=1,ntr ;
do i=is,ie
879 if (vhh(i,j) >= 0.0)
then
886 tp = tr(m)%t(i,j_up+1,k) ; tc = tr(m)%t(i,j_up,k) ; tm = tr(m)%t(i,j_up-1,k)
889 al = ( 5.*tc + ( 2.*tm - tp ) )/6.
890 al = max( min(tc,tm), al) ; al = min( max(tc,tm), al)
891 ar = ( 5.*tc + ( 2.*tp - tm ) )/6.
892 ar = max( min(tc,tp), ar) ; ar = min( max(tc,tp), ar)
894 al = 0.5 * ((tm + tc) + (slope_y(i,m,j_up-1) - slope_y(i,m,j_up)) / 3.)
895 ar = 0.5 * ((tc + tp) + (slope_y(i,m,j_up) - slope_y(i,m,j_up+1)) / 3.)
898 da = ar - al ; ma = 0.5*( ar + al )
899 if (g%mask2dCv(i,j_up)*g%mask2dCv(i,j_up-1)*(tp-tc)*(tc-tm) <= 0.)
then
901 elseif ( da*(tc-ma) > (da*da)/6. )
then
903 elseif ( da*(tc-ma) < - (da*da)/6. )
then
907 a6 = 6.*tc - 3. * (ar + al)
909 if (vhh(i,j) >= 0.0)
then
910 flux_y(i,m,j) = vhh(i,j)*( ar - 0.5 * cfl(i) * ( &
911 ( ar - al ) - a6 * ( 1. - 2./3. * cfl(i) ) ) )
913 flux_y(i,m,j) = vhh(i,j)*( al + 0.5 * cfl(i) * ( &
914 ( ar - al ) + a6 * ( 1. - 2./3. * cfl(i) ) ) )
918 do m=1,ntr ;
do i=is,ie
919 if (vhh(i,j) >= 0.0)
then
929 flux_y(i,m,j) = vhh(i,j)*( tc + 0.5 * slope_y(i,m,j) * ( 1. - cfl(i) ) )
941 tc = tr(m)%t(i,j+1,k)
942 flux_y(i,m,j) = vhh(i,j)*( tc - 0.5 * slope_y(i,m,j+1) * ( 1. - cfl(i) ) )
949 if (
associated(obc))
then ;
if (obc%OBC_pe)
then
950 if (obc%specified_v_BCs_exist_globally .or. obc%open_v_BCs_exist_globally)
then
951 do n=1,obc%number_of_segments
952 segment=>obc%segment(n)
953 if (.not. segment%specified) cycle
954 if (.not.
associated(segment%tr_Reg)) cycle
955 if (obc%segment(n)%is_N_or_S)
then
956 if (j >= segment%HI%JsdB .and. j<= segment%HI%JedB)
then
957 do i=segment%HI%isd,segment%HI%ied
960 if ((vhr(i,j,k) > 0.0) .and. (segment%direction == obc_direction_s) .or. &
961 (vhr(i,j,k) < 0.0) .and. (segment%direction == obc_direction_n))
then
962 vhh(i,j) = vhr(i,j,k)
964 if (
associated(segment%tr_Reg%Tr(m)%t))
then
965 flux_y(i,m,j) = vhh(i,j)*obc%segment(n)%tr_Reg%Tr(m)%tres(i,j,k)
966 else ; flux_y(i,m,j) = vhh(i,j)*obc%segment(n)%tr_Reg%Tr(m)%OBC_inflow_conc ;
endif
975 if (obc%open_v_BCs_exist_globally)
then
976 do n=1,obc%number_of_segments
977 segment=>obc%segment(n)
978 if (segment%specified) cycle
979 if (.not.
associated(segment%tr_Reg)) cycle
980 if (segment%is_N_or_S .and. (j >= segment%HI%JsdB .and. j<= segment%HI%JedB))
then
981 do i=segment%HI%isd,segment%HI%ied
983 if ((vhr(i,j,k) > 0.0) .and. (g%mask2dT(i,j) < 0.5) .or. &
984 (vhr(i,j,k) < 0.0) .and. (g%mask2dT(i,j+1) < 0.5))
then
985 vhh(i,j) = vhr(i,j,k)
987 if (
associated(segment%tr_Reg%Tr(m)%t))
then
988 flux_y(i,m,j) = vhh(i,j)*segment%tr_Reg%Tr(m)%tres(i,j,k)
989 else ; flux_y(i,m,j) = vhh(i,j)*segment%tr_Reg%Tr(m)%OBC_inflow_conc ;
endif
999 do i=is,ie ; vhh(i,j) = 0.0 ;
enddo
1000 do m=1,ntr ;
do i=is,ie ; flux_y(i,m,j) = 0.0 ;
enddo ;
enddo
1003 do j=js-1,je ;
do i=is,ie
1004 vhr(i,j,k) = vhr(i,j,k) - vhh(i,j)
1005 if (abs(vhr(i,j,k)) < vh_neglect(i,j)) vhr(i,j,k) = 0.0
1010 do j=js,je ;
if (do_j_tr(j))
then
1012 if ((vhh(i,j) /= 0.0) .or. (vhh(i,j-1) /= 0.0))
then
1014 hlst(i) = hprev(i,j,k)
1015 hprev(i,j,k) = max(hprev(i,j,k) - (vhh(i,j) - vhh(i,j-1)), 0.0)
1016 if (hprev(i,j,k) <= 0.0)
then ; do_i(i) = .false.
1017 elseif (hprev(i,j,k) < h_neglect*g%areaT(i,j))
then
1018 hlst(i) = hlst(i) + (h_neglect*g%areaT(i,j) - hprev(i,j,k))
1019 ihnew(i) = 1.0 / (h_neglect*g%areaT(i,j))
1020 else ; ihnew(i) = 1.0 / hprev(i,j,k) ;
endif
1021 else ; do_i(i) = .false. ;
endif
1026 do i=is,ie ;
if (do_i(i))
then
1027 tr(m)%t(i,j,k) = (tr(m)%t(i,j,k) * hlst(i) - &
1028 (flux_y(i,m,j) - flux_y(i,m,j-1))) * ihnew(i)
1032 if (
associated(tr(m)%ad_y))
then ;
do i=is,ie ;
if (do_i(i))
then
1033 tr(m)%ad_y(i,j,k) = tr(m)%ad_y(i,j,k) + flux_y(i,m,j)*idt
1034 endif ;
enddo ;
endif
1035 if (
associated(tr(m)%ad2d_y))
then ;
do i=is,ie ;
if (do_i(i))
then
1036 tr(m)%ad2d_y(i,j) = tr(m)%ad2d_y(i,j) + flux_y(i,m,j)*idt
1037 endif ;
enddo ;
endif
1041 if (
associated(tr(m)%advection_xy))
then
1042 do i=is,ie ;
if (do_i(i))
then
1043 tr(m)%advection_xy(i,j,k) = tr(m)%advection_xy(i,j,k) - (flux_y(i,m,j) - flux_y(i,m,j-1))* idt * &
1055 type(time_type),
target,
intent(in) :: time
1059 type(
diag_ctrl),
target,
intent(inout) :: diag
1062 integer,
save :: init_calls = 0
1065 # include "version_variable.h"
1066 character(len=40) :: mdl =
"MOM_tracer_advect"
1067 character(len=256) :: mesg
1069 if (
associated(cs))
then
1070 call mom_error(warning,
"tracer_advect_init called with associated control structure.")
1079 call get_param(param_file, mdl,
"DT", cs%dt, fail_if_missing=.true., &
1080 desc=
"The (baroclinic) dynamics time step.", units=
"s", scale=us%s_to_T)
1081 call get_param(param_file, mdl,
"DEBUG", cs%debug, default=.false.)
1082 call get_param(param_file, mdl,
"TRACER_ADVECTION_SCHEME", mesg, &
1083 desc=
"The horizontal transport scheme for tracers:\n"//&
1084 " PLM - Piecewise Linear Method\n"//&
1085 " PPM:H3 - Piecewise Parabolic Method (Huyhn 3rd order)\n"// &
1086 " PPM - Piecewise Parabolic Method (Colella-Woodward)" &
1088 select case (trim(mesg))
1093 cs%useHuynh = .true.
1096 cs%useHuynh = .false.
1098 call mom_error(fatal,
"MOM_tracer_advect, tracer_advect_init: "//&
1099 "Unknown TRACER_ADVECTION_SCHEME = "//trim(mesg))
1102 id_clock_advect = cpu_clock_id(
'(Ocean advect tracer)', grain=clock_module)
1103 id_clock_pass = cpu_clock_id(
'(Ocean tracer halo updates)', grain=clock_routine)
1104 id_clock_sync = cpu_clock_id(
'(Ocean tracer global synch)', grain=clock_routine)
1112 if (
associated(cs))
deallocate(cs)