56
57
64#ifdef NESTING
66#endif
67#ifdef WEC
69#endif
71
72
73
74 integer, intent(in) :: ng, tile
75 integer, intent(in) :: LBi, UBi, LBj, UBj
76 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
77 integer, intent(in) :: krhs, kstp, kout
78
79#ifdef ASSUMED_SHAPE
80 real(r8), intent(in) :: vbar(LBi:,LBj:,:)
81 real(r8), intent(in) :: zeta(LBi:,LBj:,:)
82
83 real(r8), intent(inout) :: ubar(LBi:,LBj:,:)
84#else
85 real(r8), intent(in) :: vbar(LBi:UBi,LBj:UBj,:)
86 real(r8), intent(in) :: zeta(LBi:UBi,LBj:UBj,:)
87
88 real(r8), intent(inout) :: ubar(LBi:UBi,LBj:UBj,:)
89#endif
90
91
92
93 integer :: Imin, Imax
94 integer :: i, j, know
95#ifdef NESTING
96 integer :: Idg, Jdg, cr, dg, m, rg, tnew, told
97#endif
98
99 real(r8), parameter :: eps = 1.0e-20_r8
100
101 real(r8) :: Ce, Cx, Zx
102 real(r8) :: bry_pgr, bry_cor, bry_str, bry_val
103 real(r8) :: cff, cff1, cff2, cff3, dt2d, dUde, dUdt, dUdx
104 real(r8) :: obc_in, obc_out, phi, tau
105#if defined ATM_PRESS && defined PRESS_COMPENSATE
106 real(r8) :: OneAtm, fac
107#endif
108
109 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: grad
110
111#include "set_bounds.h"
112
113
114
115
116
117#if defined STEP2D_FB_AB3_AM4 || defined STEP2D_FB_LF_AM3
118 know=kstp
120#else
121 IF (first_2d_step) THEN
122 know=krhs
125 know=krhs
127 ELSE
128 know=kstp
130 END IF
131#endif
132#if defined ATM_PRESS && defined PRESS_COMPENSATE
133 oneatm=1013.25_r8
134 fac=100.0_r8/(
g*
rho0)
135#endif
136
137
138
139
140
141 IF (
domain(ng)%Western_Edge(tile))
THEN
142
143
144
146 DO j=jstr,jend+1
147 grad(istr ,j)=ubar(istr ,j ,know)- &
148 & ubar(istr ,j-1,know)
149 grad(istr+1,j)=ubar(istr+1,j ,know)- &
150 & ubar(istr+1,j-1,know)
151 END DO
152 DO j=jstr,jend
154 dudt=ubar(istr+1,j,know)-ubar(istr+1,j,kout)
155 dudx=ubar(istr+1,j,kout)-ubar(istr+2,j,kout)
156
159 obc_out=0.5_r8* &
160 & (
clima(ng)%M2nudgcof(istr-1,j)+ &
161 &
clima(ng)%M2nudgcof(istr ,j))
162 obc_in =
obcfac(ng)*obc_out
163 ELSE
166 END IF
167 IF ((dudt*dudx).lt.0.0_r8) THEN
168 tau=obc_in
169 ELSE
170 tau=obc_out
171 END IF
172#ifdef IMPLICIT_NUDGING
173 IF (tau.gt.0.0_r8) tau=1.0_r8/tau
174#else
175 tau=tau*dt2d
176#endif
177 END IF
178
179 IF ((dudt*dudx).lt.0.0_r8) dudt=0.0_r8
180 IF ((dudt*(grad(istr+1,j )+ &
181 & grad(istr+1,j+1))).gt.0.0_r8) THEN
182 dude=grad(istr+1,j )
183 ELSE
184 dude=grad(istr+1,j+1)
185 END IF
186 cff=max(dudx*dudx+dude*dude,eps)
187 cx=dudt*dudx
188#ifdef RADIATION_2D
189 ce=min(cff,max(dudt*dude,-cff))
190#else
191 ce=0.0_r8
192#endif
193#if defined CELERITY_WRITE && defined FORWARD_WRITE
197#endif
198 ubar(istr,j,kout)=(cff*ubar(istr ,j,know)+ &
199 & cx *ubar(istr+1,j,kout)- &
200 & max(ce,0.0_r8)*grad(istr,j )- &
201 & min(ce,0.0_r8)*grad(istr,j+1))/ &
202 & (cff+cx)
203
205#ifdef IMPLICIT_NUDGING
206 phi=
dt(ng)/(tau+
dt(ng))
207 ubar(istr,j,kout)=(1.0_r8-phi)*ubar(istr,j,kout)+ &
209#else
210 ubar(istr,j,kout)=ubar(istr,j,kout)+ &
212 & ubar(istr,j,know))
213#endif
214 END IF
215#ifdef MASKING
216 ubar(istr,j,kout)=ubar(istr,j,kout)* &
217 &
grid(ng)%umask(istr,j)
218#endif
219 END IF
220 END DO
221
222
223
225 DO j=jstr,jend
227#if defined SSH_TIDES && !defined UV_TIDES
229 bry_pgr=-
g*(zeta(istr,j,know)- &
231 & 0.5_r8*
grid(ng)%pm(istr,j)
232 ELSE
233 bry_pgr=-
g*(zeta(istr ,j,know)- &
234 & zeta(istr-1,j,know))* &
235 & 0.5_r8*(
grid(ng)%pm(istr-1,j)+ &
236 &
grid(ng)%pm(istr ,j))
237 END IF
238# ifdef UV_COR
239 bry_cor=0.125_r8*(vbar(istr-1,j ,know)+ &
240 & vbar(istr-1,j+1,know)+ &
241 & vbar(istr ,j ,know)+ &
242 & vbar(istr ,j+1,know))* &
243 & (
grid(ng)%f(istr-1,j)+ &
244 &
grid(ng)%f(istr ,j))
245# else
246 bry_cor=0.0_r8
247# endif
248 cff1=1.0_r8/(0.5_r8*(
grid(ng)%h(istr-1,j)+ &
249 & zeta(istr-1,j,know)+ &
250 &
grid(ng)%h(istr ,j)+ &
251 & zeta(istr ,j,know)))
252 bry_str=cff1*(
forces(ng)%sustr(istr,j)- &
253 &
forces(ng)%bustr(istr,j))
254 cx=1.0_r8/sqrt(
g*0.5_r8*(
grid(ng)%h(istr-1,j)+ &
255 & zeta(istr-1,j,know)+ &
256 &
grid(ng)%h(istr ,j)+ &
257 & zeta(istr ,j,know)))
258 cff2=
grid(ng)%om_u(istr,j)*cx
259
260 bry_val=ubar(istr+1,j,know)+ &
261 & cff2*(bry_pgr+ &
262 & bry_cor+ &
263 & bry_str)
264#else
266#endif
267 cff=1.0_r8/(0.5_r8*(
grid(ng)%h(istr-1,j)+ &
268 & zeta(istr-1,j,know)+ &
269 &
grid(ng)%h(istr ,j)+ &
270 & zeta(istr ,j,know)))
272#if defined ATM_PRESS && defined PRESS_COMPENSATE
273 ubar(istr,j,kout)=bry_val- &
274 & cx*(0.5_r8* &
275 & (zeta(istr-1,j,know)+ &
276 & zeta(istr ,j,know)+ &
277 & fac*(
forces(ng)%Pair(istr-1,j)+ &
278 &
forces(ng)%Pair(istr ,j)- &
279 & 2.0_r8*oneatm))- &
281#else
282 ubar(istr,j,kout)=bry_val- &
283 & cx*(0.5_r8*(zeta(istr-1,j,know)+ &
284 & zeta(istr ,j,know))- &
286#endif
287#ifdef MASKING
288 ubar(istr,j,kout)=ubar(istr,j,kout)* &
289 &
grid(ng)%umask(istr,j)
290#endif
291 END IF
292 END DO
293
294
295
297 DO j=jstr,jend
299#if defined SSH_TIDES && !defined UV_TIDES
301 bry_pgr=-
g*(zeta(istr,j,know)- &
303 & 0.5_r8*
grid(ng)%pm(istr,j)
304 ELSE
305 bry_pgr=-
g*(zeta(istr ,j,know)- &
306 & zeta(istr-1,j,know))* &
307 & 0.5_r8*(
grid(ng)%pm(istr-1,j)+ &
308 &
grid(ng)%pm(istr ,j))
309 END IF
310# ifdef UV_COR
311 bry_cor=0.125_r8*(vbar(istr-1,j ,know)+ &
312 & vbar(istr-1,j+1,know)+ &
313 & vbar(istr ,j ,know)+ &
314 & vbar(istr ,j+1,know))* &
315 & (
grid(ng)%f(istr-1,j)+ &
316 &
grid(ng)%f(istr ,j))
317# else
318 bry_cor=0.0_r8
319# endif
320 cff1=1.0_r8/(0.5_r8*(
grid(ng)%h(istr-1,j)+ &
321 & zeta(istr-1,j,know)+ &
322 &
grid(ng)%h(istr ,j)+ &
323 & zeta(istr ,j,know)))
324 bry_str=cff1*(
forces(ng)%sustr(istr,j)- &
325 &
forces(ng)%bustr(istr,j))
326 cx=1.0_r8/sqrt(
g*0.5_r8*(
grid(ng)%h(istr-1,j)+ &
327 & zeta(istr-1,j,know)+ &
328 &
grid(ng)%h(istr ,j)+ &
329 & zeta(istr ,j,know)))
330 cff2=
grid(ng)%om_u(istr,j)*cx
331
332 bry_val=ubar(istr+1,j,know)+ &
333 & cff2*(bry_pgr+ &
334 & bry_cor+ &
335 & bry_str)
336#else
338#endif
339#ifdef WET_DRY
340 cff=0.5_r8*(
grid(ng)%h(istr-1,j)+ &
341 & zeta(istr-1,j,know)+ &
342 &
grid(ng)%h(istr ,j)+ &
343 & zeta(istr ,j,know))
344#else
345 cff=0.5_r8*(
grid(ng)%h(istr-1,j)+ &
346 &
grid(ng)%h(istr ,j))
347#endif
349 cx=dt2d*cff1*cff*0.5_r8*(
grid(ng)%pm(istr-1,j)+ &
350 &
grid(ng)%pm(istr ,j))
351 zx=(0.5_r8+cx)*zeta(istr ,j,know)+ &
352 & (0.5_r8-cx)*zeta(istr-1,j,know)
354 cff2=(1.0_r8-
co/cx)**2
355 cff3=zeta(istr,j,kout)+ &
356 & cx*zeta(istr-1,j,know)- &
357 & (1.0_r8+cx)*zeta(istr,j,know)
358 zx=zx+cff2*cff3
359 END IF
360 ubar(istr,j,kout)=0.5_r8* &
361 & ((1.0_r8-cx)*ubar(istr,j,know)+ &
362 & cx*ubar(istr+1,j,know)+ &
363 & bry_val- &
364 & cff1*(zx-
boundary(ng)%zeta_west(j)))
365#ifdef MASKING
366 ubar(istr,j,kout)=ubar(istr,j,kout)* &
367 &
grid(ng)%umask(istr,j)
368#endif
369 END IF
370 END DO
371
372
373
375 DO j=jstr,jend
377 ubar(istr,j,kout)=
boundary(ng)%ubar_west(j)
378#ifdef MASKING
379 ubar(istr,j,kout)=ubar(istr,j,kout)* &
380 &
grid(ng)%umask(istr,j)
381#endif
382 END IF
383 END DO
384
385
386
388 DO j=jstr,jend
390 ubar(istr,j,kout)=ubar(istr+1,j,kout)
391#ifdef MASKING
392 ubar(istr,j,kout)=ubar(istr,j,kout)* &
393 &
grid(ng)%umask(istr,j)
394#endif
395 END IF
396 END DO
397
398
399
401 DO j=jstr,jend
404 bry_pgr=-
g*(zeta(istr,j,know)- &
406 & 0.5_r8*
grid(ng)%pm(istr,j)
407 ELSE
408 bry_pgr=-
g*(zeta(istr ,j,know)- &
409 & zeta(istr-1,j,know))* &
410 & 0.5_r8*(
grid(ng)%pm(istr-1,j)+ &
411 &
grid(ng)%pm(istr ,j))
412 END IF
413#ifdef UV_COR
414 bry_cor=0.125_r8*(vbar(istr-1,j ,know)+ &
415 & vbar(istr-1,j+1,know)+ &
416 & vbar(istr ,j ,know)+ &
417 & vbar(istr ,j+1,know))* &
418 & (
grid(ng)%f(istr-1,j)+ &
419 &
grid(ng)%f(istr ,j))
420#else
421 bry_cor=0.0_r8
422#endif
423 cff=1.0_r8/(0.5_r8*(
grid(ng)%h(istr-1,j)+ &
424 & zeta(istr-1,j,know)+ &
425 &
grid(ng)%h(istr ,j)+ &
426 & zeta(istr ,j,know)))
427 bry_str=cff*(
forces(ng)%sustr(istr,j)- &
428 &
forces(ng)%bustr(istr,j))
429 ubar(istr,j,kout)=ubar(istr,j,know)+ &
430 & dt2d*(bry_pgr+ &
431 & bry_cor+ &
432 & bry_str)
433#ifdef MASKING
434 ubar(istr,j,kout)=ubar(istr,j,kout)* &
435 &
grid(ng)%umask(istr,j)
436#endif
437 END IF
438 END DO
439
440
441
443 DO j=jstr,jend
445 ubar(istr,j,kout)=0.0_r8
446 END IF
447 END DO
448
449#ifdef NESTING
450
451
452
453
459 & (rg.eq.ng).and.(
dxmax(dg).gt.
dxmax(rg)))
THEN
462 DO j=jstr,jend
466 cff=0.5_r8*
grid(ng)%on_u(istr,j)* &
467 & (
grid(ng)%h(istr-1,j)+zeta(istr-1,j,kout)+ &
468 &
grid(ng)%h(istr ,j)+zeta(istr ,j,kout))
470 bry_val=cff1*
refined(cr)%U2d_flux(1,m,tnew)/cff
471# ifdef WEC
472 bry_val=bry_val-
ocean(ng)%ubar_stokes(istr,j)
473# endif
474# ifdef MASKING
475 bry_val=bry_val*
grid(ng)%umask(istr,j)
476# endif
477# ifdef NESTING_DEBUG
479# endif
480 ubar(istr,j,kout)=bry_val
481 END DO
482 END IF
483 END DO
484#endif
485 END IF
486 END IF
487
488
489
490
491
492 IF (
domain(ng)%Eastern_Edge(tile))
THEN
493
494
495
497 DO j=jstr,jend+1
498 grad(iend ,j)=ubar(iend ,j ,know)- &
499 & ubar(iend ,j-1,know)
500 grad(iend+1,j)=ubar(iend+1,j ,know)- &
501 & ubar(iend+1,j-1,know)
502 END DO
503 DO j=jstr,jend
505 dudt=ubar(iend,j,know)-ubar(iend ,j,kout)
506 dudx=ubar(iend,j,kout)-ubar(iend-1,j,kout)
507
510 obc_out=0.5_r8* &
511 & (
clima(ng)%M2nudgcof(iend ,j)+ &
512 &
clima(ng)%M2nudgcof(iend+1,j))
513 obc_in =
obcfac(ng)*obc_out
514 ELSE
517 END IF
518 IF ((dudt*dudx).lt.0.0_r8) THEN
519 tau=obc_in
520 ELSE
521 tau=obc_out
522 END IF
523#ifdef IMPLICIT_NUDGING
524 IF (tau.gt.0.0_r8) tau=1.0_r8/tau
525#else
526 tau=tau*dt2d
527#endif
528 END IF
529
530 IF ((dudt*dudx).lt.0.0_r8) dudt=0.0_r8
531 IF ((dudt*(grad(iend,j )+ &
532 & grad(iend,j+1))).gt.0.0_r8) THEN
533 dude=grad(iend,j)
534 ELSE
535 dude=grad(iend,j+1)
536 END IF
537 cff=max(dudx*dudx+dude*dude,eps)
538 cx=dudt*dudx
539#ifdef RADIATION_2D
540 ce=min(cff,max(dudt*dude,-cff))
541#else
542 ce=0.0_r8
543#endif
544#if defined CELERITY_WRITE && defined FORWARD_WRITE
548#endif
549 ubar(iend+1,j,kout)=(cff*ubar(iend+1,j,know)+ &
550 & cx *ubar(iend ,j,kout)- &
551 & max(ce,0.0_r8)*grad(iend+1,j )- &
552 & min(ce,0.0_r8)*grad(iend+1,j+1))/ &
553 & (cff+cx)
554
556#ifdef IMPLICIT_NUDGING
557 phi=
dt(ng)/(tau+
dt(ng))
558 ubar(iend+1,j,kout)=(1.0_r8-phi)*ubar(iend+1,j,kout)+ &
560#else
561 ubar(iend+1,j,kout)=ubar(iend+1,j,kout)+ &
563 & ubar(iend+1,j,know))
564#endif
565 END IF
566#ifdef MASKING
567 ubar(iend+1,j,kout)=ubar(iend+1,j,kout)* &
568 &
grid(ng)%umask(iend+1,j)
569#endif
570 END IF
571 END DO
572
573
574
576 DO j=jstr,jend
578#if defined SSH_TIDES && !defined UV_TIDES
581 & zeta(iend,j,know))* &
582 & 0.5_r8*
grid(ng)%pm(iend,j)
583 ELSE
584 bry_pgr=-
g*(zeta(iend+1,j,know)- &
585 & zeta(iend ,j,know))* &
586 & 0.5_r8*(
grid(ng)%pm(iend ,j)+ &
587 &
grid(ng)%pm(iend+1,j))
588 END IF
589# ifdef UV_COR
590 bry_cor=0.125_r8*(vbar(iend ,j ,know)+ &
591 & vbar(iend ,j+1,know)+ &
592 & vbar(iend+1,j ,know)+ &
593 & vbar(iend+1,j+1,know))* &
594 & (
grid(ng)%f(iend ,j)+ &
595 &
grid(ng)%f(iend+1,j))
596# else
597 bry_cor=0.0_r8
598# endif
599 cff1=1.0_r8/(0.5_r8*(
grid(ng)%h(iend ,j)+ &
600 & zeta(iend ,j,know)+ &
601 &
grid(ng)%h(iend+1,j)+ &
602 & zeta(iend+1,j,know)))
603 bry_str=cff1*(
forces(ng)%sustr(iend+1,j)- &
604 &
forces(ng)%bustr(iend+1,j))
605 cx=1.0_r8/sqrt(
g*0.5_r8*(
grid(ng)%h(iend+1,j)+ &
606 & zeta(iend+1,j,know)+ &
607 &
grid(ng)%h(iend ,j)+ &
608 & zeta(iend ,j,know)))
609 cff2=
grid(ng)%om_u(iend+1,j)*cx
610
611 bry_val=ubar(iend,j,know)+ &
612 & cff2*(bry_pgr+ &
613 & bry_cor+ &
614 & bry_str)
615#else
617#endif
618 cff=1.0_r8/(0.5_r8*(
grid(ng)%h(iend ,j)+ &
619 & zeta(iend ,j,know)+ &
620 &
grid(ng)%h(iend+1,j)+ &
621 & zeta(iend+1,j,know)))
623#if defined ATM_PRESS && defined PRESS_COMPENSATE
624 ubar(iend+1,j,kout)=bry_val+ &
625 & cx*(0.5_r8* &
626 & (zeta(iend ,j,know)+ &
627 & zeta(iend+1,j,know)+ &
628 & fac*(
forces(ng)%Pair(iend ,j)+ &
629 &
forces(ng)%Pair(iend+1,j)- &
630 & 2.0_r8*oneatm))- &
632#else
633 ubar(iend+1,j,kout)=bry_val+ &
634 & cx*(0.5_r8*(zeta(iend ,j,know)+ &
635 & zeta(iend+1,j,know))- &
637#endif
638#ifdef MASKING
639 ubar(iend+1,j,kout)=ubar(iend+1,j,kout)* &
640 &
grid(ng)%umask(iend+1,j)
641#endif
642 END IF
643 END DO
644
645
646
648 DO j=jstr,jend
650#if defined SSH_TIDES && !defined UV_TIDES
653 & zeta(iend,j,know))* &
654 & 0.5_r8*
grid(ng)%pm(iend,j)
655 ELSE
656 bry_pgr=-
g*(zeta(iend+1,j,know)- &
657 & zeta(iend ,j,know))* &
658 & 0.5_r8*(
grid(ng)%pm(iend ,j)+ &
659 &
grid(ng)%pm(iend+1,j))
660 END IF
661# ifdef UV_COR
662 bry_cor=0.125_r8*(vbar(iend ,j ,know)+ &
663 & vbar(iend ,j+1,know)+ &
664 & vbar(iend+1,j ,know)+ &
665 & vbar(iend+1,j+1,know))* &
666 & (
grid(ng)%f(iend ,j)+ &
667 &
grid(ng)%f(iend+1,j))
668# else
669 bry_cor=0.0_r8
670# endif
671 cff1=1.0_r8/(0.5_r8*(
grid(ng)%h(iend ,j)+ &
672 & zeta(iend ,j,know)+ &
673 &
grid(ng)%h(iend+1,j)+ &
674 & zeta(iend+1,j,know)))
675 bry_str=cff1*(
forces(ng)%sustr(iend+1,j)- &
676 &
forces(ng)%bustr(iend+1,j))
677 cx=1.0_r8/sqrt(
g*0.5_r8*(
grid(ng)%h(iend+1,j)+ &
678 & zeta(iend+1,j,know)+ &
679 &
grid(ng)%h(iend ,j)+ &
680 & zeta(iend ,j,know)))
681 cff2=
grid(ng)%om_u(iend+1,j)*cx
682
683 bry_val=ubar(iend,j,know)+ &
684 & cff2*(bry_pgr+ &
685 & bry_cor+ &
686 & bry_str)
687#else
689#endif
690#ifdef WET_DRY
691 cff=0.5_r8*(
grid(ng)%h(iend ,j)+ &
692 & zeta(iend ,j,know)+ &
693 &
grid(ng)%h(iend+1,j)+ &
694 & zeta(iend+1,j,know))
695#else
696 cff=0.5_r8*(
grid(ng)%h(iend ,j)+ &
697 &
grid(ng)%h(iend+1,j))
698#endif
700 cx=dt2d*cff1*cff*0.5_r8*(
grid(ng)%pm(iend ,j)+ &
701 &
grid(ng)%pm(iend+1,j))
702 zx=(0.5_r8+cx)*zeta(iend ,j,know)+ &
703 & (0.5_r8-cx)*zeta(iend+1,j,know)
705 cff2=(1.0_r8-
co/cx)**2
706 cff3=zeta(iend,j,kout)+ &
707 & cx*zeta(iend+1,j,know)- &
708 & (1.0_r8+cx)*zeta(iend,j,know)
709 zx=zx+cff2*cff3
710 END IF
711 ubar(iend+1,j,kout)=0.5_r8* &
712 & ((1.0_r8-cx)*ubar(iend+1,j,know)+ &
713 & cx*ubar(iend,j,know)+ &
714 & bry_val+ &
715 & cff1*(zx-
boundary(ng)%zeta_east(j)))
716#ifdef MASKING
717 ubar(iend+1,j,kout)=ubar(iend+1,j,kout)* &
718 &
grid(ng)%umask(iend+1,j)
719#endif
720 END IF
721 END DO
722
723
724
726 DO j=jstr,jend
728 ubar(iend+1,j,kout)=
boundary(ng)%ubar_east(j)
729#ifdef MASKING
730 ubar(iend+1,j,kout)=ubar(iend+1,j,kout)* &
731 &
grid(ng)%umask(iend+1,j)
732#endif
733 END IF
734 END DO
735
736
737
739 DO j=jstr,jend
741 ubar(iend+1,j,kout)=ubar(iend,j,kout)
742#ifdef MASKING
743 ubar(iend+1,j,kout)=ubar(iend+1,j,kout)* &
744 &
grid(ng)%umask(iend+1,j)
745#endif
746 END IF
747 END DO
748
749
750
752 DO j=jstr,jend
756 & zeta(iend,j,know))* &
757 & 0.5_r8*
grid(ng)%pm(iend,j)
758 ELSE
759 bry_pgr=-
g*(zeta(iend+1,j,know)- &
760 & zeta(iend ,j,know))* &
761 & 0.5_r8*(
grid(ng)%pm(iend ,j)+ &
762 &
grid(ng)%pm(iend+1,j))
763 END IF
764#ifdef UV_COR
765 bry_cor=0.125_r8*(vbar(iend ,j ,know)+ &
766 & vbar(iend ,j+1,know)+ &
767 & vbar(iend+1,j ,know)+ &
768 & vbar(iend+1,j+1,know))* &
769 & (
grid(ng)%f(iend ,j)+ &
770 &
grid(ng)%f(iend+1,j))
771#else
772 bry_cor=0.0_r8
773#endif
774 cff=1.0_r8/(0.5_r8*(
grid(ng)%h(iend ,j)+ &
775 & zeta(iend ,j,know)+ &
776 &
grid(ng)%h(iend+1,j)+ &
777 & zeta(iend+1,j,know)))
778 bry_str=cff*(
forces(ng)%sustr(iend+1,j)- &
779 &
forces(ng)%bustr(iend+1,j))
780 ubar(iend+1,j,kout)=ubar(iend+1,j,know)+ &
781 & dt2d*(bry_pgr+ &
782 & bry_cor+ &
783 & bry_str)
784#ifdef MASKING
785 ubar(iend+1,j,kout)=ubar(iend+1,j,kout)* &
786 &
grid(ng)%umask(iend+1,j)
787#endif
788 END IF
789 END DO
790
791
792
794 DO j=jstr,jend
796 ubar(iend+1,j,kout)=0.0_r8
797 END IF
798 END DO
799
800#ifdef NESTING
801
802
803
804
810 & (rg.eq.ng).and.(
dxmax(dg).gt.
dxmax(rg)))
THEN
813 DO j=jstr,jend
817 cff=0.5_r8*
grid(ng)%on_u(iend+1,j)* &
818 & (
grid(ng)%h(iend+1,j)+zeta(iend+1,j,kout)+ &
819 &
grid(ng)%h(iend ,j)+zeta(iend ,j,kout))
821 bry_val=cff1*
refined(cr)%U2d_flux(1,m,tnew)/cff
822# ifdef WEC
823 bry_val=bry_val-
ocean(ng)%ubar_stokes(iend+1,j)
824# endif
825# ifdef MASKING
826 bry_val=bry_val*
grid(ng)%umask(iend+1,j)
827# endif
828# ifdef NESTING_DEBUG
830# endif
831 ubar(iend+1,j,kout)=bry_val
832 END DO
833 END IF
834 END DO
835#endif
836 END IF
837 END IF
838
839
840
841
842
843 IF (
domain(ng)%Southern_Edge(tile))
THEN
844
845
846
848 DO i=istru-1,iend
849 grad(i,jstr-1)=ubar(i+1,jstr-1,know)- &
850 & ubar(i ,jstr-1,know)
851 grad(i,jstr )=ubar(i+1,jstr ,know)- &
852 & ubar(i ,jstr ,know)
853 END DO
854 DO i=istru,iend
856 dudt=ubar(i,jstr,know)-ubar(i,jstr ,kout)
857 dude=ubar(i,jstr,kout)-ubar(i,jstr+1,kout)
858
861 obc_out=0.5_r8* &
862 & (
clima(ng)%M2nudgcof(i-1,jstr-1)+ &
863 &
clima(ng)%M2nudgcof(i ,jstr-1))
864 obc_in =
obcfac(ng)*obc_out
865 ELSE
868 END IF
869 IF ((dudt*dude).lt.0.0_r8) THEN
870 tau=obc_in
871 ELSE
872 tau=obc_out
873 END IF
874#ifdef IMPLICIT_NUDGING
875 IF (tau.gt.0.0_r8) tau=1.0_r8/tau
876#else
877 tau=tau*dt2d
878#endif
879 END IF
880
881 IF ((dudt*dude).lt.0.0_r8) dudt=0.0_r8
882 IF ((dudt*(grad(i-1,jstr)+ &
883 & grad(i ,jstr))).gt.0.0_r8) THEN
884 dudx=grad(i-1,jstr)
885 ELSE
886 dudx=grad(i ,jstr)
887 END IF
888 cff=max(dudx*dudx+dude*dude,eps)
889#ifdef RADIATION_2D
890 cx=min(cff,max(dudt*dudx,-cff))
891#else
892 cx=0.0_r8
893#endif
894 ce=dudt*dude
895#if defined CELERITY_WRITE && defined FORWARD_WRITE
899#endif
900 ubar(i,jstr-1,kout)=(cff*ubar(i,jstr-1,know)+ &
901 & ce *ubar(i,jstr ,kout)- &
902 & max(cx,0.0_r8)*grad(i-1,jstr-1)- &
903 & min(cx,0.0_r8)*grad(i ,jstr-1))/ &
904 & (cff+ce)
905
907#ifdef IMPLICIT_NUDGING
908 phi=
dt(ng)/(tau+
dt(ng))
909 ubar(i,jstr-1,kout)=(1.0_r8-phi)*ubar(i,jstr-1,kout)+ &
911#else
912 ubar(i,jstr-1,kout)=ubar(i,jstr-1,kout)+ &
913 & tau*(
boundary(ng)%ubar_south(i)- &
914 & ubar(i,jstr-1,know))
915#endif
916 END IF
917#ifdef MASKING
918 ubar(i,jstr-1,kout)=ubar(i,jstr-1,kout)* &
919 &
grid(ng)%umask(i,jstr-1)
920#endif
921 END IF
922 END DO
923
924
925
929 DO i=istru,iend
931 cff=dt2d*0.5_r8*(
grid(ng)%pn(i-1,jstr)+ &
932 &
grid(ng)%pn(i ,jstr))
933 cff1=sqrt(
g*0.5_r8*(
grid(ng)%h(i-1,jstr)+ &
934 & zeta(i-1,jstr,know)+ &
935 &
grid(ng)%h(i ,jstr)+ &
936 & zeta(i ,jstr,know)))
937 ce=cff*cff1
938 cff2=1.0_r8/(1.0_r8+ce)
939 ubar(i,jstr-1,kout)=cff2*(ubar(i,jstr-1,know)+ &
940 & ce*ubar(i,jstr,kout))
941#ifdef MASKING
942 ubar(i,jstr-1,kout)=ubar(i,jstr-1,kout)* &
943 &
grid(ng)%umask(i,jstr-1)
944#endif
945 END IF
946 END DO
947
948
949
951 DO i=istru,iend
953 ubar(i,jstr-1,kout)=
boundary(ng)%ubar_south(i)
954#ifdef MASKING
955 ubar(i,jstr-1,kout)=ubar(i,jstr-1,kout)* &
956 &
grid(ng)%umask(i,jstr-1)
957#endif
958 END IF
959 END DO
960
961
962
964 DO i=istru,iend
966 ubar(i,jstr-1,kout)=ubar(i,jstr,kout)
967#ifdef MASKING
968 ubar(i,jstr-1,kout)=ubar(i,jstr-1,kout)* &
969 &
grid(ng)%umask(i,jstr-1)
970#endif
971 END IF
972 END DO
973
974
975
976
979 imin=istru
980 imax=iend
981 ELSE
982 imin=istr
983 imax=iendr
984 END IF
985 DO i=imin,imax
987 ubar(i,jstr-1,kout)=
gamma2(ng)*ubar(i,jstr,kout)
988#ifdef MASKING
989 ubar(i,jstr-1,kout)=ubar(i,jstr-1,kout)* &
990 &
grid(ng)%umask(i,jstr-1)
991#endif
992 END IF
993 END DO
994 END IF
995 END IF
996
997
998
999
1000
1001 IF (
domain(ng)%Northern_Edge(tile))
THEN
1002
1003
1004
1006 DO i=istru-1,iend
1007 grad(i,jend )=ubar(i+1,jend ,know)- &
1008 & ubar(i ,jend ,know)
1009 grad(i,jend+1)=ubar(i+1,jend+1,know)- &
1010 & ubar(i ,jend+1,know)
1011 END DO
1012 DO i=istru,iend
1014 dudt=ubar(i,jend,know)-ubar(i,jend ,kout)
1015 dude=ubar(i,jend,kout)-ubar(i,jend-1,kout)
1016
1019 obc_out=0.5_r8* &
1020 & (
clima(ng)%M2nudgcof(i-1,jend+1)+ &
1021 &
clima(ng)%M2nudgcof(i ,jend+1))
1022 obc_in =
obcfac(ng)*obc_out
1023 ELSE
1026 END IF
1027 IF ((dudt*dude).lt.0.0_r8) THEN
1028 tau=obc_in
1029 ELSE
1030 tau=obc_out
1031 END IF
1032#ifdef IMPLICIT_NUDGING
1033 IF (tau.gt.0.0_r8) tau=1.0_r8/tau
1034#else
1035 tau=tau*dt2d
1036#endif
1037 END IF
1038
1039 IF ((dudt*dude).lt.0.0_r8) dudt=0.0_r8
1040 IF ((dudt*(grad(i-1,jend)+ &
1041 & grad(i ,jend))).gt.0.0_r8) THEN
1042 dudx=grad(i-1,jend)
1043 ELSE
1044 dudx=grad(i ,jend)
1045 END IF
1046 cff=max(dudx*dudx+dude*dude,eps)
1047#ifdef RADIATION_2D
1048 cx=min(cff,max(dudt*dudx,-cff))
1049#else
1050 cx=0.0_r8
1051#endif
1052 ce=dudt*dude
1053#if defined CELERITY_WRITE && defined FORWARD_WRITE
1057#endif
1058 ubar(i,jend+1,kout)=(cff*ubar(i,jend+1,know)+ &
1059 & ce *ubar(i,jend ,kout)- &
1060 & max(cx,0.0_r8)*grad(i-1,jend+1)- &
1061 & min(cx,0.0_r8)*grad(i ,jend+1))/ &
1062 & (cff+ce)
1063
1065#ifdef IMPLICIT_NUDGING
1066 phi=
dt(ng)/(tau+
dt(ng))
1067 ubar(i,jend+1,kout)=(1.0_r8-phi)*ubar(i,jend+1,kout)+ &
1069#else
1070 ubar(i,jend+1,kout)=ubar(i,jend+1,kout)+ &
1071 & tau*(
boundary(ng)%ubar_north(i)- &
1072 & ubar(i,jend+1,know))
1073#endif
1074 END IF
1075#ifdef MASKING
1076 ubar(i,jend+1,kout)=ubar(i,jend+1,kout)* &
1077 &
grid(ng)%umask(i,jend+1)
1078#endif
1079 END IF
1080 END DO
1081
1082
1083
1087 DO i=istru,iend
1089 cff=dt2d*0.5_r8*(
grid(ng)%pn(i-1,jend)+ &
1090 &
grid(ng)%pn(i ,jend))
1091 cff1=sqrt(
g*0.5_r8*(
grid(ng)%h(i-1,jend)+ &
1092 & zeta(i-1,jend,know)+ &
1093 &
grid(ng)%h(i ,jend)+ &
1094 & zeta(i ,jend,know)))
1095 ce=cff*cff1
1096 cff2=1.0_r8/(1.0_r8+ce)
1097 ubar(i,jend+1,kout)=cff2*(ubar(i,jend+1,know)+ &
1098 & ce*ubar(i,jend,kout))
1099#ifdef MASKING
1100 ubar(i,jend+1,kout)=ubar(i,jend+1,kout)* &
1101 &
grid(ng)%umask(i,jend+1)
1102#endif
1103 END IF
1104 END DO
1105
1106
1107
1109 DO i=istru,iend
1111 ubar(i,jend+1,kout)=
boundary(ng)%ubar_north(i)
1112#ifdef MASKING
1113 ubar(i,jend+1,kout)=ubar(i,jend+1,kout)* &
1114 &
grid(ng)%umask(i,jend+1)
1115#endif
1116 END IF
1117 END DO
1118
1119
1120
1122 DO i=istru,iend
1124 ubar(i,jend+1,kout)=ubar(i,jend,kout)
1125#ifdef MASKING
1126 ubar(i,jend+1,kout)=ubar(i,jend+1,kout)* &
1127 &
grid(ng)%umask(i,jend+1)
1128#endif
1129 END IF
1130 END DO
1131
1132
1133
1134
1137 imin=istru
1138 imax=iend
1139 ELSE
1140 imin=istr
1141 imax=iendr
1142 END IF
1143 DO i=imin,imax
1145 ubar(i,jend+1,kout)=
gamma2(ng)*ubar(i,jend,kout)
1146#ifdef MASKING
1147 ubar(i,jend+1,kout)=ubar(i,jend+1,kout)* &
1148 &
grid(ng)%umask(i,jend+1)
1149#endif
1150 END IF
1151 END DO
1152 END IF
1153 END IF
1154
1155
1156
1157
1158
1160 IF (
domain(ng)%SouthWest_Corner(tile))
THEN
1163 ubar(istr,jstr-1,kout)=0.5_r8*(ubar(istr+1,jstr-1,kout)+ &
1164 & ubar(istr ,jstr ,kout))
1165 END IF
1166 END IF
1167 IF (
domain(ng)%SouthEast_Corner(tile))
THEN
1170 ubar(iend+1,jstr-1,kout)=0.5_r8*(ubar(iend ,jstr-1,kout)+ &
1171 & ubar(iend+1,jstr ,kout))
1172 END IF
1173 END IF
1174 IF (
domain(ng)%NorthWest_Corner(tile))
THEN
1177 ubar(istr,jend+1,kout)=0.5_r8*(ubar(istr ,jend ,kout)+ &
1178 & ubar(istr+1,jend+1,kout))
1179 END IF
1180 END IF
1181 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
1184 ubar(iend+1,jend+1,kout)=0.5_r8*(ubar(iend+1,jend ,kout)+ &
1185 & ubar(iend ,jend+1,kout))
1186 END IF
1187 END IF
1188 END IF
1189
1190#if defined WET_DRY
1191
1192
1193
1194
1195
1197 IF (
domain(ng)%Western_Edge(tile))
THEN
1198 DO j=jstr,jend
1201 cff1=abs(abs(
grid(ng)%umask_wet(istr,j))-1.0_r8)
1202 cff2=0.5_r8+dsign(0.5_r8,ubar(istr,j,kout))* &
1203 &
grid(ng)%umask_wet(istr,j)
1204 cff=0.5_r8*
grid(ng)%umask_wet(istr,j)*cff1+ &
1205 & cff2*(1.0_r8-cff1)
1206 ubar(istr,j,kout)=ubar(istr,j,kout)*cff
1207 END IF
1208 END DO
1209 END IF
1210 IF (
domain(ng)%Eastern_Edge(tile))
THEN
1211 DO j=jstr,jend
1214 cff1=abs(abs(
grid(ng)%umask_wet(iend+1,j))-1.0_r8)
1215 cff2=0.5_r8+dsign(0.5_r8,ubar(iend+1,j,kout))* &
1216 &
grid(ng)%umask_wet(iend+1,j)
1217 cff=0.5_r8*
grid(ng)%umask_wet(iend+1,j)*cff1+ &
1218 & cff2*(1.0_r8-cff1)
1219 ubar(iend+1,j,kout)=ubar(iend+1,j,kout)*cff
1220 END IF
1221 END DO
1222 END IF
1223 END IF
1224
1226 IF (
domain(ng)%Southern_Edge(tile))
THEN
1227 DO i=istru,iend
1230 cff1=abs(abs(
grid(ng)%umask_wet(i,jstr-1))-1.0_r8)
1231 cff2=0.5_r8+dsign(0.5_r8,ubar(i,jstr-1,kout))* &
1232 &
grid(ng)%umask_wet(i,jstr-1)
1233 cff=0.5_r8*
grid(ng)%umask_wet(i,jstr-1)*cff1+ &
1234 & cff2*(1.0_r8-cff1)
1235 ubar(i,jstr-1,kout)=ubar(i,jstr-1,kout)*cff
1236 END IF
1237 END DO
1238 END IF
1239 IF (
domain(ng)%Northern_Edge(tile))
THEN
1240 DO i=istr,iend
1243 cff1=abs(abs(
grid(ng)%umask_wet(i,jend+1))-1.0_r8)
1244 cff2=0.5_r8+dsign(0.5_r8,ubar(i,jend+1,kout))* &
1245 &
grid(ng)%umask_wet(i,jend+1)
1246 cff=0.5_r8*
grid(ng)%umask_wet(i,jend+1)*cff1+ &
1247 & cff2*(1.0_r8-cff1)
1248 ubar(i,jend+1,kout)=ubar(i,jend+1,kout)*cff
1249 END IF
1250 END DO
1251 END IF
1252 END IF
1253
1255 IF (
domain(ng)%SouthWest_Corner(tile))
THEN
1260 cff1=abs(abs(
grid(ng)%umask_wet(istr,jstr-1))-1.0_r8)
1261 cff2=0.5_r8+dsign(0.5_r8,ubar(istr,jstr-1,kout))* &
1262 &
grid(ng)%umask_wet(istr,jstr-1)
1263 cff=0.5_r8*
grid(ng)%umask_wet(istr,jstr-1)*cff1+ &
1264 & cff2*(1.0_r8-cff1)
1265 ubar(istr,jstr-1,kout)=ubar(istr,jstr-1,kout)*cff
1266 END IF
1267 END IF
1268 IF (
domain(ng)%SouthEast_Corner(tile))
THEN
1273 cff1=abs(abs(
grid(ng)%umask_wet(iend+1,jstr-1))-1.0_r8)
1274 cff2=0.5_r8+dsign(0.5_r8,ubar(iend+1,jstr-1,kout))* &
1275 &
grid(ng)%umask_wet(iend+1,jstr-1)
1276 cff=0.5_r8*
grid(ng)%umask_wet(iend+1,jstr-1)*cff1+ &
1277 & cff2*(1.0_r8-cff1)
1278 ubar(iend+1,jstr-1,kout)=ubar(iend+1,jstr-1,kout)*cff
1279 END IF
1280 END IF
1281 IF (
domain(ng)%NorthWest_Corner(tile))
THEN
1286 cff1=abs(abs(
grid(ng)%umask_wet(istr,jend+1))-1.0_r8)
1287 cff2=0.5_r8+dsign(0.5_r8,ubar(istr,jend+1,kout))* &
1288 &
grid(ng)%umask_wet(istr,jend+1)
1289 cff=0.5_r8*
grid(ng)%umask_wet(istr,jend+1)*cff1+ &
1290 & cff2*(1.0_r8-cff1)
1291 ubar(istr,jend+1,kout)=ubar(istr,jend+1,kout)*cff
1292 END IF
1293 END IF
1294 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
1299 cff1=abs(abs(
grid(ng)%umask_wet(iend+1,jend+1))-1.0_r8)
1300 cff2=0.5_r8+dsign(0.5_r8,ubar(iend+1,jend+1,kout))* &
1301 &
grid(ng)%umask_wet(iend+1,jend+1)
1302 cff=0.5_r8*
grid(ng)%umask_wet(iend+1,jend+1)*cff1+ &
1303 & cff2*(1.0_r8-cff1)
1304 ubar(iend+1,jend+1,kout)=ubar(iend+1,jend+1,kout)*cff
1305 END IF
1306 END IF
1307 END IF
1308#endif
1309
1310 RETURN
type(t_boundary), dimension(:), allocatable boundary
type(t_apply), dimension(:), allocatable lbc_apply
type(t_clima), dimension(:), allocatable clima
type(t_forces), dimension(:), allocatable forces
type(t_grid), dimension(:), allocatable grid
type(t_bcp), dimension(:,:), allocatable bry_contact
integer, dimension(:), allocatable rollingindex
type(t_refined), dimension(:), allocatable refined
type(t_ngc), dimension(:), allocatable ucontact
type(t_lbc), dimension(:,:,:), allocatable lbc
type(t_domain), dimension(:), allocatable domain
logical, dimension(:), allocatable lnudgem2clm
real(dp), dimension(:), allocatable dt
logical, dimension(:), allocatable ewperiodic
logical, dimension(:), allocatable nsperiodic
logical, dimension(:), allocatable predictor_2d_step
real(dp), dimension(:), allocatable obcfac
real(r8), dimension(:), allocatable gamma2
integer, parameter isouth
real(dp), dimension(:), allocatable dxmax
real(dp), dimension(:), allocatable dtfast
real(dp), dimension(:,:), allocatable m2obc_out
logical, dimension(:), allocatable refinedgrid
integer, parameter inorth
real(dp), dimension(:,:), allocatable m2obc_in