61 & LBi, UBi, LBj, UBj, &
62 & IminS, ImaxS, JminS, JmaxS, &
69 integer,
intent(in) :: ng, tile
70 integer,
intent(in) :: lbi, ubi, lbj, ubj
71 integer,
intent(in) :: imins, imaxs, jmins, jmaxs
72 integer,
intent(in) ::
krhs,
kstp, kout
75 real(r8),
intent(inout) :: zeta(lbi:,lbj:,:)
77 real(r8),
intent(inout) :: zeta(lbi:ubi,lbj:ubj,:)
84 real(r8),
parameter :: eps =1.0e-20_r8
87 real(r8) :: cff, cff1, cff2, dt2d, dzde, dzdt, dzdx, tau
89 real(r8),
dimension(IminS:ImaxS,JminS:JmaxS) :: grad
91#include "set_bounds.h"
97#if defined STEP2D_FB_AB3_AM4 || defined STEP2D_FB_LF_AM3
101 IF (first_2d_step)
THEN
117 IF (
domain(ng)%Western_Edge(tile))
THEN
123 grad(istr-1,j)=zeta(istr-1,j ,know)- &
124 & zeta(istr-1,j-1,know)
126 grad(istr-1,j)=grad(istr-1,j)*
grid(ng)%vmask(istr-1,j)
128 grad(istr,j)=zeta(istr,j ,know)- &
129 & zeta(istr,j-1,know)
131 grad(istr,j)=grad(istr,j)*
grid(ng)%vmask(istr,j)
136 dzdt=zeta(istr,j,know)-zeta(istr ,j,kout)
137 dzdx=zeta(istr,j,kout)-zeta(istr+1,j,kout)
140 IF ((dzdt*dzdx).lt.0.0_r8)
THEN
148 IF ((dzdt*dzdx).lt.0.0_r8) dzdt=0.0_r8
149 IF ((dzdt*(grad(istr,j)+grad(istr,j+1))).gt.0.0_r8)
THEN
154 cff=max(dzdx*dzdx+dzde*dzde,eps)
157 ce=min(cff,max(dzdt*dzde,-cff))
161#if defined CELERITY_WRITE && defined FORWARD_WRITE
166 zeta(istr-1,j,kout)=(cff*zeta(istr-1,j,know)+ &
167 & cx *zeta(istr ,j,kout)- &
168 & max(ce,0.0_r8)*grad(istr-1,j )- &
169 & min(ce,0.0_r8)*grad(istr-1,j+1))/ &
173 zeta(istr-1,j,kout)=zeta(istr-1,j,kout)+ &
175 & zeta(istr-1,j,know))
178 zeta(istr-1,j,kout)=zeta(istr-1,j,kout)* &
179 &
grid(ng)%rmask(istr-1,j)
189 cff=dt2d*
grid(ng)%pm(istr,j)
191 cff1=sqrt(
g*(max(
grid(ng)%h(istr,j)+ &
192 & zeta(istr,j,know),
dcrit(ng))))
194 cff1=sqrt(
g*(
grid(ng)%h(istr,j)+ &
195 & zeta(istr,j,know)))
198 zeta(istr-1,j,kout)=(1.0_r8-cx)*zeta(istr-1,j,know)+ &
199 & cx*zeta(istr,j,know)
201 zeta(istr-1,j,kout)=zeta(istr-1,j,kout)* &
202 &
grid(ng)%rmask(istr-1,j)
212 cff=dt2d*
grid(ng)%pm(istr,j)
214 cff1=sqrt(
g*(max(
grid(ng)%h(istr,j)+ &
215 & zeta(istr,j,know),
dcrit(ng))))
217 cff1=sqrt(
g*(
grid(ng)%h(istr,j)+ &
218 & zeta(istr,j,know)))
221 cff2=1.0_r8/(1.0_r8+cx)
222 zeta(istr-1,j,kout)=cff2*(zeta(istr-1,j,know)+ &
223 & cx*zeta(istr,j,kout))
225 zeta(istr-1,j,kout)=zeta(istr-1,j,kout)* &
226 &
grid(ng)%rmask(istr-1,j)
236 zeta(istr-1,j,kout)=
boundary(ng)%zeta_west(j)
238 zeta(istr-1,j,kout)=zeta(istr-1,j,kout)* &
239 &
grid(ng)%rmask(istr-1,j)
249 zeta(istr-1,j,kout)=zeta(istr,j,kout)
251 zeta(istr-1,j,kout)=zeta(istr-1,j,kout)* &
252 &
grid(ng)%rmask(istr-1,j)
262 zeta(istr-1,j,kout)=zeta(istr,j,kout)
264 zeta(istr-1,j,kout)=zeta(istr-1,j,kout)* &
265 &
grid(ng)%rmask(istr-1,j)
276 IF (
domain(ng)%Eastern_Edge(tile))
THEN
282 grad(iend ,j)=zeta(iend ,j ,know)- &
283 & zeta(iend ,j-1,know)
285 grad(iend ,j)=grad(iend ,j)*
grid(ng)%vmask(iend ,j)
287 grad(iend+1,j)=zeta(iend+1,j ,know)- &
288 & zeta(iend+1,j-1,know)
290 grad(iend+1,j)=grad(iend+1,j)*
grid(ng)%vmask(iend+1,j)
295 dzdt=zeta(iend,j,know)-zeta(iend ,j,kout)
296 dzdx=zeta(iend,j,kout)-zeta(iend-1,j,kout)
299 IF ((dzdt*dzdx).lt.0.0_r8)
THEN
307 IF ((dzdt*dzdx).lt.0.0_r8) dzdt=0.0_r8
308 IF ((dzdt*(grad(iend,j)+grad(iend,j+1))).gt.0.0_r8)
THEN
313 cff=max(dzdx*dzdx+dzde*dzde,eps)
316 ce=min(cff,max(dzdt*dzde,-cff))
320#if defined CELERITY_WRITE && defined FORWARD_WRITE
325 zeta(iend+1,j,kout)=(cff*zeta(iend+1,j,know)+ &
326 & cx *zeta(iend ,j,kout)- &
327 & max(ce,0.0_r8)*grad(iend+1,j )- &
328 & min(ce,0.0_r8)*grad(iend+1,j+1))/ &
332 zeta(iend+1,j,kout)=zeta(iend+1,j,kout)+ &
334 & zeta(iend+1,j,know))
337 zeta(iend+1,j,kout)=zeta(iend+1,j,kout)* &
338 &
grid(ng)%rmask(iend+1,j)
348 cff=dt2d*
grid(ng)%pm(iend,j)
350 cff1=sqrt(
g*(max(
grid(ng)%h(iend,j)+ &
351 & zeta(iend,j,know),
dcrit(ng))))
353 cff1=sqrt(
g*(
grid(ng)%h(iend,j)+ &
354 & zeta(iend,j,know)))
357 zeta(iend+1,j,kout)=(1.0_r8-cx)*zeta(iend+1,j,know)+ &
358 & cx*zeta(iend,j,know)
360 zeta(iend+1,j,kout)=zeta(iend+1,j,kout)* &
361 &
grid(ng)%rmask(iend+1,j)
371 cff=dt2d*
grid(ng)%pm(iend,j)
373 cff1=sqrt(
g*(max(
grid(ng)%h(iend,j)+ &
374 & zeta(iend,j,know),
dcrit(ng))))
376 cff1=sqrt(
g*(
grid(ng)%h(iend,j)+ &
377 & zeta(iend,j,know)))
380 cff2=1.0_r8/(1.0_r8+cx)
381 zeta(iend+1,j,kout)=cff2*(zeta(iend+1,j,know)+ &
382 & cx*zeta(iend,j,kout))
384 zeta(iend+1,j,kout)=zeta(iend+1,j,kout)* &
385 &
grid(ng)%rmask(iend+1,j)
395 zeta(iend+1,j,kout)=
boundary(ng)%zeta_east(j)
397 zeta(iend+1,j,kout)=zeta(iend+1,j,kout)* &
398 &
grid(ng)%rmask(iend+1,j)
408 zeta(iend+1,j,kout)=zeta(iend,j,kout)
410 zeta(iend+1,j,kout)=zeta(iend+1,j,kout)* &
411 &
grid(ng)%rmask(iend+1,j)
421 zeta(iend+1,j,kout)=zeta(iend,j,kout)
423 zeta(iend+1,j,kout)=zeta(iend+1,j,kout)* &
424 &
grid(ng)%rmask(iend+1,j)
435 IF (
domain(ng)%Southern_Edge(tile))
THEN
441 grad(i,jstr )=zeta(i ,jstr,know)- &
442 & zeta(i-1,jstr,know)
444 grad(i,jstr )=grad(i,jstr )*
grid(ng)%umask(i,jstr )
446 grad(i,jstr-1)=zeta(i ,jstr-1,know)- &
447 & zeta(i-1,jstr-1,know)
449 grad(i,jstr-1)=grad(i,jstr-1)*
grid(ng)%umask(i,jstr-1)
454 dzdt=zeta(i,jstr,know)-zeta(i,jstr ,kout)
455 dzde=zeta(i,jstr,kout)-zeta(i,jstr-1,kout)
458 IF ((dzdt*dzde).lt.0.0_r8)
THEN
466 IF ((dzdt*dzde).lt.0.0_r8) dzdt=0.0_r8
467 IF ((dzdt*(grad(i,jstr)+grad(i+1,jstr))).gt.0.0_r8)
THEN
472 cff=max(dzdx*dzdx+dzde*dzde,eps)
474 cx=min(cff,max(dzdt*dzdx,-cff))
479#if defined CELERITY_WRITE && defined FORWARD_WRITE
484 zeta(i,jstr-1,kout)=(cff*zeta(i,jstr-1,know)+ &
485 & ce *zeta(i,jstr ,kout)- &
486 & max(cx,0.0_r8)*grad(i ,jstr)- &
487 & min(cx,0.0_r8)*grad(i+1,jstr))/ &
491 zeta(i,jstr-1,kout)=zeta(i,jstr-1,kout)+ &
492 & tau*(
boundary(ng)%zeta_south(i)- &
493 & zeta(i,jstr-1,know))
496 zeta(i,jstr-1,kout)=zeta(i,jstr-1,kout)* &
497 &
grid(ng)%rmask(i,jstr-1)
507 cff=dt2d*
grid(ng)%pn(i,jstr)
509 cff1=sqrt(
g*(max(
grid(ng)%h(i,jstr)+ &
510 & zeta(i,jstr,know),
dcrit(ng))))
512 cff1=sqrt(
g*(
grid(ng)%h(i,jstr)+ &
513 & zeta(i,jstr,know)))
516 zeta(i,jstr-1,kout)=(1.0_r8-ce)*zeta(i,jstr-1,know)+ &
517 & ce*zeta(i,jstr,know)
519 zeta(i,jstr-1,kout)=zeta(i,jstr-1,kout)* &
520 &
grid(ng)%rmask(i,jstr-1)
530 cff=dt2d*
grid(ng)%pn(i,jstr)
532 cff1=sqrt(
g*(max(
grid(ng)%h(i,jstr)+ &
533 & zeta(i,jstr,know),
dcrit(ng))))
535 cff1=sqrt(
g*(
grid(ng)%h(i,jstr)+ &
536 & zeta(i,jstr,know)))
539 cff2=1.0_r8/(1.0_r8+ce)
540 zeta(i,jstr-1,kout)=cff2*(zeta(i,jstr-1,know)+ &
541 & ce*zeta(i,jstr,kout))
543 zeta(i,jstr-1,kout)=zeta(i,jstr-1,kout)* &
544 &
grid(ng)%rmask(i,jstr-1)
554 zeta(i,jstr-1,kout)=
boundary(ng)%zeta_south(i)
556 zeta(i,jstr-1,kout)=zeta(i,jstr-1,kout)* &
557 &
grid(ng)%rmask(i,jstr-1)
567 zeta(i,jstr-1,kout)=zeta(i,jstr,kout)
569 zeta(i,jstr-1,kout)=zeta(i,jstr-1,kout)* &
570 &
grid(ng)%rmask(i,jstr-1)
580 zeta(i,jstr-1,kout)=zeta(i,jstr,kout)
582 zeta(i,jstr-1,kout)=zeta(i,jstr-1,kout)* &
583 &
grid(ng)%rmask(i,jstr-1)
594 IF (
domain(ng)%Northern_Edge(tile))
THEN
600 grad(i,jend )=zeta(i ,jend ,know)- &
601 & zeta(i-1,jend ,know)
603 grad(i,jend )=grad(i,jend )*
grid(ng)%umask(i,jend )
605 grad(i,jend+1)=zeta(i ,jend+1,know)- &
606 & zeta(i-1,jend+1,know)
608 grad(i,jend+1)=grad(i,jend+1)*
grid(ng)%umask(i,jend+1)
613 dzdt=zeta(i,jend,know)-zeta(i,jend ,kout)
614 dzde=zeta(i,jend,kout)-zeta(i,jend-1,kout)
617 IF ((dzdt*dzde).lt.0.0_r8)
THEN
625 IF ((dzdt*dzde).lt.0.0_r8) dzdt=0.0_r8
626 IF ((dzdt*(grad(i,jend)+grad(i+1,jend))).gt.0.0_r8)
THEN
631 cff=max(dzdx*dzdx+dzde*dzde,eps)
633 cx=min(cff,max(dzdt*dzdx,-cff))
638#if defined CELERITY_WRITE && defined FORWARD_WRITE
643 zeta(i,jend+1,kout)=(cff*zeta(i,jend+1,know)+ &
644 & ce *zeta(i,jend ,kout)- &
645 & max(cx,0.0_r8)*grad(i ,jend+1)- &
646 & min(cx,0.0_r8)*grad(i+1,jend+1))/ &
650 zeta(i,jend+1,kout)=zeta(i,jend+1,kout)+ &
651 & tau*(
boundary(ng)%zeta_north(i)- &
652 & zeta(i,jend+1,know))
655 zeta(i,jend+1,kout)=zeta(i,jend+1,kout)* &
656 &
grid(ng)%rmask(i,jend+1)
666 cff=dt2d*
grid(ng)%pn(i,jend)
668 cff1=sqrt(
g*(max(
grid(ng)%h(i,jend)+ &
669 & zeta(i,jend,know),
dcrit(ng))))
671 cff1=sqrt(
g*(
grid(ng)%h(i,jend)+ &
672 & zeta(i,jend,know)))
675 zeta(i,jend+1,kout)=(1.0_r8-ce)*zeta(i,jend+1,know)+ &
676 & ce*zeta(i,jend,know)
678 zeta(i,jend+1,kout)=zeta(i,jend+1,kout)* &
679 &
grid(ng)%rmask(i,jend+1)
689 cff=dt2d*
grid(ng)%pn(i,jend)
691 cff1=sqrt(
g*(max(
grid(ng)%h(i,jend)+ &
692 & zeta(i,jend,know),
dcrit(ng))))
694 cff1=sqrt(
g*(
grid(ng)%h(i,jend)+ &
695 & zeta(i,jend,know)))
698 cff2=1.0_r8/(1.0_r8+ce)
699 zeta(i,jend+1,kout)=cff2*(zeta(i,jend+1,know)+ &
700 & ce*zeta(i,jend,kout))
702 zeta(i,jend+1,kout)=zeta(i,jend+1,kout)* &
703 &
grid(ng)%rmask(i,jend+1)
713 zeta(i,jend+1,kout)=
boundary(ng)%zeta_north(i)
715 zeta(i,jend+1,kout)=zeta(i,jend+1,kout)* &
716 &
grid(ng)%rmask(i,jend+1)
726 zeta(i,jend+1,kout)=zeta(i,jend,kout)
728 zeta(i,jend+1,kout)=zeta(i,jend+1,kout)* &
729 &
grid(ng)%rmask(i,jend+1)
739 zeta(i,jend+1,kout)=zeta(i,jend,kout)
741 zeta(i,jend+1,kout)=zeta(i,jend+1,kout)* &
742 &
grid(ng)%rmask(i,jend+1)
754 IF (
domain(ng)%SouthWest_Corner(tile))
THEN
757 zeta(istr-1,jstr-1,kout)=0.5_r8*(zeta(istr ,jstr-1,kout)+ &
758 & zeta(istr-1,jstr ,kout))
761 IF (
domain(ng)%SouthEast_Corner(tile))
THEN
764 zeta(iend+1,jstr-1,kout)=0.5_r8*(zeta(iend ,jstr-1,kout)+ &
765 & zeta(iend+1,jstr ,kout))
768 IF (
domain(ng)%NorthWest_Corner(tile))
THEN
771 zeta(istr-1,jend+1,kout)=0.5_r8*(zeta(istr-1,jend ,kout)+ &
772 & zeta(istr ,jend+1,kout))
775 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
778 zeta(iend+1,jend+1,kout)=0.5_r8*(zeta(iend+1,jend ,kout)+ &
779 & zeta(iend ,jend+1,kout))
792 IF (
domain(ng)%Western_Edge(tile))
THEN
795 IF (zeta(istr-1,j,kout).le. &
797 zeta(istr-1,j,kout)=cff-
grid(ng)%h(istr-1,j)
802 IF (
domain(ng)%Eastern_Edge(tile))
THEN
805 IF (zeta(iend+1,j,kout).le. &
807 zeta(iend+1,j,kout)=cff-
grid(ng)%h(iend+1,j)
815 IF (
domain(ng)%Southern_Edge(tile))
THEN
818 IF (zeta(i,jstr-1,kout).le. &
820 zeta(i,jstr-1,kout)=cff-
grid(ng)%h(i,jstr-1)
825 IF (
domain(ng)%Northern_Edge(tile))
THEN
828 IF (zeta(i,jend+1,kout).le. &
830 zeta(i,jend+1,kout)=cff-
grid(ng)%h(i,jend+1)
838 IF (
domain(ng)%SouthWest_Corner(tile))
THEN
841 IF (zeta(istr-1,jstr-1,kout).le. &
842 & (
dcrit(ng)-
grid(ng)%h(istr-1,jstr-1)))
THEN
843 zeta(istr-1,jstr-1,kout)=cff-
grid(ng)%h(istr-1,jstr-1)
847 IF (
domain(ng)%SouthEast_Corner(tile))
THEN
850 IF (zeta(iend+1,jstr-1,kout).le. &
851 & (
dcrit(ng)-
grid(ng)%h(iend+1,jstr-1)))
THEN
852 zeta(iend+1,jstr-1,kout)=cff-
grid(ng)%h(iend+1,jstr-1)
856 IF (
domain(ng)%NorthWest_Corner(tile))
THEN
859 IF (zeta(istr-1,jend+1,kout).le. &
860 & (
dcrit(ng)-
grid(ng)%h(istr-1,jend+1)))
THEN
861 zeta(istr-1,jend+1,kout)=cff-
grid(ng)%h(istr-1,jend+1)
865 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
868 IF (zeta(iend+1,jend+1,kout).le. &
869 & (
dcrit(ng)-
grid(ng)%h(iend+1,jend+1)))
THEN
870 zeta(iend+1,jend+1,kout)=cff-
grid(ng)%h(iend+1,jend+1)