57
58
65#ifdef NESTING
67#endif
68#ifdef WEC
70#endif
72
73
74
75 integer, intent(in) :: ng, tile
76 integer, intent(in) :: LBi, UBi, LBj, UBj
77 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
78 integer, intent(in) :: krhs, kstp, kout
79
80#ifdef ASSUMED_SHAPE
81 real(r8), intent(in) :: ubar(LBi:,LBj:,:)
82 real(r8), intent(in) :: zeta(LBi:,LBj:,:)
83
84 real(r8), intent(inout) :: vbar(LBi:,LBj:,:)
85#else
86 real(r8), intent(in) :: ubar(LBi:UBi,LBj:UBj,:)
87 real(r8), intent(in) :: zeta(LBi:UBi,LBj:UBj,:)
88
89 real(r8), intent(inout) :: vbar(LBi:UBi,LBj:UBj,:)
90#endif
91
92
93
94 integer :: Jmin, Jmax
95 integer :: i, j, know
96#ifdef NESTING
97 integer :: Idg, Jdg, cr, dg, m, rg, tnew, told
98#endif
99
100 real(r8), parameter :: eps = 1.0e-20_r8
101
102 real(r8) :: Ce, Cx, Ze
103 real(r8) :: bry_pgr, bry_cor, bry_str, bry_wec, bry_val
104 real(r8) :: cff, cff1, cff2, cff3, dt2d, dVde, dVdt, dVdx
105 real(r8) :: obc_in, obc_out, phi, tau
106#ifdef WEC_VF
107 real(r8), parameter :: Lwave_min = 1.0_r8
108 real(r8) :: sigma, osigma, waven, waveny
109#endif
110#if defined ATM_PRESS && defined PRESS_COMPENSATE
111 real(r8) :: OneAtm, fac
112#endif
113
114 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: grad
115
116#include "set_bounds.h"
117
118
119
120
121
122#if defined STEP2D_FB_AB3_AM4 || defined STEP2D_FB_LF_AM3
123 know=kstp
125#else
126 IF (first_2d_step) THEN
127 know=krhs
130 know=krhs
132 ELSE
133 know=kstp
135 END IF
136#endif
137#if defined ATM_PRESS && defined PRESS_COMPENSATE
138 oneatm=1013.25_r8
139 fac=100.0_r8/(
g*
rho0)
140#endif
141
142
143
144
145
146 IF (
domain(ng)%Southern_Edge(tile))
THEN
147
148
149
151 DO i=istr,iend+1
152 grad(i,jstr )=vbar(i ,jstr ,know)- &
153 & vbar(i-1,jstr ,know)
154 grad(i,jstr+1)=vbar(i ,jstr+1,know)- &
155 & vbar(i-1,jstr+1,know)
156 END DO
157 DO i=istr,iend
159 dvdt=vbar(i,jstr+1,know)-vbar(i,jstr+1,kout)
160 dvde=vbar(i,jstr+1,kout)-vbar(i,jstr+2,kout)
161
164 obc_out=0.5_r8* &
165 & (
clima(ng)%M2nudgcof(i,jstr-1)+ &
166 &
clima(ng)%M2nudgcof(i,jstr ))
167 obc_in =
obcfac(ng)*obc_out
168 ELSE
171 END IF
172 IF ((dvdt*dvde).lt.0.0_r8) THEN
173 tau=obc_in
174 ELSE
175 tau=obc_out
176 END IF
177#ifdef IMPLICIT_NUDGING
178 IF (tau.gt.0.0_r8) tau=1.0_r8/tau
179#else
180 tau=tau*dt2d
181#endif
182 END IF
183
184 IF ((dvdt*dvde).lt.0.0_r8) dvdt=0.0_r8
185 IF ((dvdt*(grad(i ,jstr+1)+ &
186 & grad(i+1,jstr+1))).gt.0.0_r8) THEN
187 dvdx=grad(i ,jstr+1)
188 ELSE
189 dvdx=grad(i+1,jstr+1)
190 END IF
191 cff=max(dvdx*dvdx+dvde*dvde,eps)
192#ifdef RADIATION_2D
193 cx=min(cff,max(dvdt*dvdx,-cff))
194#else
195 cx=0.0_r8
196#endif
197 ce=dvdt*dvde
198#if defined CELERITY_WRITE && defined FORWARD_WRITE
202#endif
203 vbar(i,jstr,kout)=(cff*vbar(i,jstr ,know)+ &
204 & ce *vbar(i,jstr+1,kout)- &
205 & max(cx,0.0_r8)*grad(i ,jstr)- &
206 & min(cx,0.0_r8)*grad(i+1,jstr))/ &
207 & (cff+ce)
208
210#ifdef IMPLICIT_NUDGING
211 phi=
dt(ng)/(tau+
dt(ng))
212 vbar(i,jstr,kout)=(1.0_r8-phi)*vbar(i,jstr,kout)+ &
214
215#else
216 vbar(i,jstr,kout)=vbar(i,jstr,kout)+ &
217 & tau*(
boundary(ng)%vbar_south(i)- &
218 & vbar(i,jstr,know))
219#endif
220 END IF
221#ifdef MASKING
222 vbar(i,jstr,kout)=vbar(i,jstr,kout)* &
223 &
grid(ng)%vmask(i,jstr)
224#endif
225 END IF
226 END DO
227
228
229
231 DO i=istr,iend
233#if defined SSH_TIDES && !defined UV_TIDES
235 bry_pgr=-
g*(zeta(i,jstr,know)- &
237 & 0.5_r8*
grid(ng)%pn(i,jstr)
238 ELSE
239 bry_pgr=-
g*(zeta(i,jstr ,know)- &
240 & zeta(i,jstr-1,know))* &
241 & 0.5_r8*(
grid(ng)%pn(i,jstr-1)+ &
242 &
grid(ng)%pn(i,jstr ))
243 END IF
244# ifdef UV_COR
245 bry_cor=-0.125_r8*(ubar(i ,jstr-1,know)+ &
246 & ubar(i+1,jstr-1,know)+ &
247 & ubar(i ,jstr ,know)+ &
248 & ubar(i+1,jstr ,know))* &
249 & (
grid(ng)%f(i,jstr-1)+ &
250 &
grid(ng)%f(i,jstr ))
251# else
252 bry_cor=0.0_r8
253# endif
254 cff1=1.0_r8/(0.5_r8*(
grid(ng)%h(i,jstr-1)+ &
255 & zeta(i,jstr-1,know)+ &
256 &
grid(ng)%h(i,jstr )+ &
257 & zeta(i,jstr ,know)))
258 bry_str=cff1*(
forces(ng)%svstr(i,jstr)- &
259 &
forces(ng)%bvstr(i,jstr))
260 ce=1.0_r8/sqrt(
g*0.5_r8*(
grid(ng)%h(i,jstr-1)+ &
261 & zeta(i,jstr-1,know)+ &
262 &
grid(ng)%h(i,jstr )+ &
263 & zeta(i,jstr ,know)))
264# ifdef WEC_VF
265
266 bry_pgr=min(0.5_r8*(
grid(ng)%h(i,jstr)+ &
267 & zeta(i,jstr,know)), 1.0_r8)*bry_pgr
268 cff=1.0_r8/(0.5_r8*(
grid(ng)%h(i,jstr-1)+ &
269 & zeta(i,jstr-1,know)+ &
270 &
grid(ng)%h(i,jstr )+ &
271 & zeta(i,jstr ,know)))
272 cff1=1.5_r8*
pi-
forces(ng)%Dwave(i,jstr)- &
273 &
grid(ng)%angler(i,jstr)
274 waven=2.0_r8*
pi/max(
forces(ng)%Lwave(i,jstr),lwave_min)
275 waveny=waven*sin(cff1)
276 sigma=sqrt(max(
g*waven*tanh(waven/cff),eps))
277 osigma=1.0_r8/sigma
279 &
forces(ng)%Dissip_break(i,jstr)*waveny
280 bry_wec=cff1
281# ifdef WEC_ROLLER
282 bry_wec=bry_wec+osigma*
forces(ng)%Dissip_roller(i,jstr)* &
283 & waveny
284# endif
285
286# else
287 bry_wec=0.0_r8
288# endif
289 cff2=
grid(ng)%on_v(i,jstr)*ce
290
291 bry_val=vbar(i,jstr+1,know)+ &
292 & cff2*(bry_pgr+ &
293 & bry_cor+ &
294 & bry_wec+ &
295 & bry_str)
296#else
298#endif
299 cff=1.0_r8/(0.5_r8*(
grid(ng)%h(i,jstr-1)+ &
300 & zeta(i,jstr-1,know)+ &
301 &
grid(ng)%h(i,jstr )+ &
302 & zeta(i,jstr ,know)))
304#if defined ATM_PRESS && defined PRESS_COMPENSATE
305 vbar(i,jstr,kout)=bry_val- &
306 & ce*(0.5_r8* &
307 & (zeta(i,jstr-1,know)+ &
308 & zeta(i,jstr ,know)+ &
309 & fac*(
forces(ng)%Pair(i,jstr-1)+ &
310 &
forces(ng)%Pair(i,jstr )- &
311 & 2.0_r8*oneatm))- &
313#else
314 vbar(i,jstr,kout)=bry_val- &
315 & ce*(0.5_r8*(zeta(i,jstr-1,know)+ &
316 & zeta(i,jstr ,know))- &
318#endif
319#ifdef MASKING
320 vbar(i,jstr,kout)=vbar(i,jstr,kout)* &
321 &
grid(ng)%vmask(i,jstr)
322#endif
323 END IF
324 END DO
325
326
327
329 DO i=istr,iend
331#if defined SSH_TIDES && !defined UV_TIDES
333 bry_pgr=-
g*(zeta(i,jstr,know)- &
335 & 0.5_r8*
grid(ng)%pn(i,jstr)
336 ELSE
337 bry_pgr=-
g*(zeta(i,jstr ,know)- &
338 & zeta(i,jstr-1,know))* &
339 & 0.5_r8*(
grid(ng)%pn(i,jstr-1)+ &
340 &
grid(ng)%pn(i,jstr ))
341 END IF
342# ifdef UV_COR
343 bry_cor=-0.125_r8*(ubar(i ,jstr-1,know)+ &
344 & ubar(i+1,jstr-1,know)+ &
345 & ubar(i ,jstr ,know)+ &
346 & ubar(i+1,jstr ,know))* &
347 & (
grid(ng)%f(i,jstr-1)+ &
348 &
grid(ng)%f(i,jstr ))
349# else
350 bry_cor=0.0_r8
351# endif
352 cff1=1.0_r8/(0.5_r8*(
grid(ng)%h(i,jstr-1)+ &
353 & zeta(i,jstr-1,know)+ &
354 &
grid(ng)%h(i,jstr )+ &
355 & zeta(i,jstr ,know)))
356 bry_str=cff1*(
forces(ng)%svstr(i,jstr)- &
357 &
forces(ng)%bvstr(i,jstr))
358 ce=1.0_r8/sqrt(
g*0.5_r8*(
grid(ng)%h(i,jstr-1)+ &
359 & zeta(i,jstr-1,know)+ &
360 &
grid(ng)%h(i,jstr )+ &
361 & zeta(i,jstr ,know)))
362 cff2=
grid(ng)%on_v(i,jstr)*ce
363
364 bry_val=vbar(i,jstr+1,know)+ &
365 & cff2*(bry_pgr+ &
366 & bry_cor+ &
367 & bry_str)
368#else
370#endif
371#ifdef WET_DRY
372 cff=0.5_r8*(
grid(ng)%h(i,jstr-1)+ &
373 & zeta(i,jstr-1,know)+ &
374 &
grid(ng)%h(i,jstr )+ &
375 & zeta(i,jstr ,know))
376#else
377 cff=0.5_r8*(
grid(ng)%h(i,jstr-1)+ &
378 &
grid(ng)%h(i,jstr ))
379#endif
381 ce=dt2d*cff1*cff*0.5_r8*(
grid(ng)%pn(i,jstr-1)+ &
382 &
grid(ng)%pn(i,jstr ))
383 ze=(0.5_r8+ce)*zeta(i,jstr ,know)+ &
384 & (0.5_r8-ce)*zeta(i,jstr-1,know)
386 cff2=(1.0_r8-
co/ce)**2
387 cff3=zeta(i,jstr,kout)+ &
388 & ce*zeta(i,jstr-1,know)- &
389 & (1.0_r8+ce)*zeta(i,jstr,know)
390 ze=ze+cff2*cff3
391 END IF
392 vbar(i,jstr,kout)=0.5_r8* &
393 & ((1.0_r8-ce)*vbar(i,jstr,know)+ &
394 & ce*vbar(i,jstr+1,know)+ &
395 & bry_val- &
396 & cff1*(ze-
boundary(ng)%zeta_south(i)))
397#ifdef MASKING
398 vbar(i,jstr,kout)=vbar(i,jstr,kout)* &
399 &
grid(ng)%vmask(i,jstr)
400#endif
401 END IF
402 END DO
403
404
405
407 DO i=istr,iend
409 vbar(i,jstr,kout)=
boundary(ng)%vbar_south(i)
410#ifdef MASKING
411 vbar(i,jstr,kout)=vbar(i,jstr,kout)* &
412 &
grid(ng)%vmask(i,jstr)
413#endif
414 END IF
415 END DO
416
417
418
420 DO i=istr,iend
422 vbar(i,jstr,kout)=vbar(i,jstr+1,kout)
423#ifdef MASKING
424 vbar(i,jstr,kout)=vbar(i,jstr,kout)* &
425 &
grid(ng)%vmask(i,jstr)
426#endif
427 END IF
428 END DO
429
430
431
433 DO i=istr,iend
436 bry_pgr=-
g*(zeta(i,jstr,know)- &
438 & 0.5_r8*
grid(ng)%pn(i,jstr)
439 ELSE
440 bry_pgr=-
g*(zeta(i,jstr ,know)- &
441 & zeta(i,jstr-1,know))* &
442 & 0.5_r8*(
grid(ng)%pn(i,jstr-1)+ &
443 &
grid(ng)%pn(i,jstr ))
444 END IF
445#ifdef UV_COR
446 bry_cor=-0.125_r8*(ubar(i ,jstr-1,know)+ &
447 & ubar(i+1,jstr-1,know)+ &
448 & ubar(i ,jstr ,know)+ &
449 & ubar(i+1,jstr ,know))* &
450 & (
grid(ng)%f(i,jstr-1)+ &
451 &
grid(ng)%f(i,jstr ))
452#else
453 bry_cor=0.0_r8
454#endif
455 cff=1.0_r8/(0.5_r8*(
grid(ng)%h(i,jstr-1)+ &
456 & zeta(i,jstr-1,know)+ &
457 &
grid(ng)%h(i,jstr )+ &
458 & zeta(i,jstr ,know)))
459 bry_str=cff*(
forces(ng)%svstr(i,jstr)- &
460 &
forces(ng)%bvstr(i,jstr))
461 vbar(i,jstr,kout)=vbar(i,jstr,know)+ &
462 & dt2d*(bry_pgr+ &
463 & bry_cor+ &
464 & bry_str)
465#ifdef MASKING
466 vbar(i,jstr,kout)=vbar(i,jstr,kout)* &
467 &
grid(ng)%vmask(i,jstr)
468#endif
469 END IF
470 END DO
471
472
473
475 DO i=istr,iend
477 vbar(i,jstr,kout)=0.0_r8
478 END IF
479 END DO
480
481#ifdef NESTING
482
483
484
485
491 & (rg.eq.ng).and.(
dxmax(dg).gt.
dxmax(rg)))
THEN
494 DO i=istr,iend
498 cff=0.5_r8*
grid(ng)%om_v(i,jstr)* &
499 & (
grid(ng)%h(i,jstr-1)+zeta(i,jstr-1,kout)+ &
500 &
grid(ng)%h(i,jstr )+zeta(i,jstr ,kout))
502 bry_val=cff1*
refined(cr)%V2d_flux(1,m,tnew)/cff
503# ifdef WEC
504 bry_val=bry_val-
ocean(ng)%vbar_stokes(i,jstr)
505# endif
506# ifdef MASKING
507 bry_val=bry_val*
grid(ng)%vmask(i,jstr)
508# endif
509# ifdef NESTING_DEBUG
511# endif
512 vbar(i,jstr,kout)=bry_val
513 END DO
514 END IF
515 END DO
516#endif
517 END IF
518 END IF
519
520
521
522
523
524 IF (
domain(ng)%Northern_Edge(tile))
THEN
525
526
527
529 DO i=istr,iend+1
530 grad(i,jend )=vbar(i ,jend ,know)- &
531 & vbar(i-1,jend ,know)
532 grad(i,jend+1)=vbar(i ,jend+1,know)- &
533 & vbar(i-1,jend+1,know)
534 END DO
535 DO i=istr,iend
537 dvdt=vbar(i,jend,know)-vbar(i,jend ,kout)
538 dvde=vbar(i,jend,kout)-vbar(i,jend-1,kout)
539
542 obc_out=0.5_r8* &
543 & (
clima(ng)%M2nudgcof(i,jend )+ &
544 &
clima(ng)%M2nudgcof(i,jend+1))
545 obc_in =
obcfac(ng)*obc_out
546 ELSE
549 END IF
550 IF ((dvdt*dvde).lt.0.0_r8) THEN
551 tau=obc_in
552 ELSE
553 tau=obc_out
554 END IF
555#ifdef IMPLICIT_NUDGING
556 IF (tau.gt.0.0_r8) tau=1.0_r8/tau
557#else
558 tau=tau*dt2d
559#endif
560 END IF
561
562 IF ((dvdt*dvde).lt.0.0_r8) dvdt=0.0_r8
563 IF ((dvdt*(grad(i ,jend)+ &
564 & grad(i+1,jend))).gt.0.0_r8) THEN
565 dvdx=grad(i ,jend)
566 ELSE
567 dvdx=grad(i+1,jend)
568 END IF
569 cff=max(dvdx*dvdx+dvde*dvde,eps)
570#ifdef RADIATION_2D
571 cx=min(cff,max(dvdt*dvdx,-cff))
572#else
573 cx=0.0_r8
574#endif
575 ce=dvdt*dvde
576#if defined CELERITY_WRITE && defined FORWARD_WRITE
580#endif
581 vbar(i,jend+1,kout)=(cff*vbar(i,jend+1,know)+ &
582 & ce *vbar(i,jend ,kout)- &
583 & max(cx,0.0_r8)*grad(i ,jend+1)- &
584 & min(cx,0.0_r8)*grad(i+1,jend+1))/ &
585 & (cff+ce)
586
588#ifdef IMPLICIT_NUDGING
589 phi=
dt(ng)/(tau+
dt(ng))
590 vbar(i,jend+1,kout)=(1.0_r8-phi)*vbar(i,jend+1,kout)+ &
592
593#else
594 vbar(i,jend+1,kout)=vbar(i,jend+1,kout)+ &
595 & tau*(
boundary(ng)%vbar_north(i)- &
596 & vbar(i,jend+1,know))
597#endif
598 END IF
599#ifdef MASKING
600 vbar(i,jend+1,kout)=vbar(i,jend+1,kout)* &
601 &
grid(ng)%vmask(i,jend+1)
602#endif
603 END IF
604 END DO
605
606
607
609 DO i=istr,iend
611#if defined SSH_TIDES && !defined UV_TIDES
613 bry_pgr=-
g*(
boundary(ng)%zeta_north(i)- &
614 & zeta(i,jend,know))* &
615 & 0.5_r8*
grid(ng)%pn(i,jend)
616 ELSE
617 bry_pgr=-
g*(zeta(i,jend+1,know)- &
618 & zeta(i,jend ,know))* &
619 & 0.5_r8*(
grid(ng)%pn(i,jend )+ &
620 &
grid(ng)%pn(i,jend+1))
621 END IF
622# ifdef UV_COR
623 bry_cor=-0.125_r8*(ubar(i ,jend ,know)+ &
624 & ubar(i+1,jend ,know)+ &
625 & ubar(i ,jend+1,know)+ &
626 & ubar(i+1,jend+1,know))* &
627 & (
grid(ng)%f(i,jend )+ &
628 &
grid(ng)%f(i,jend+1))
629# else
630 bry_cor=0.0_r8
631# endif
632 cff1=1.0_r8/(0.5_r8*(
grid(ng)%h(i,jend )+ &
633 & zeta(i,jend ,know)+ &
634 &
grid(ng)%h(i,jend+1)+ &
635 & zeta(i,jend+1,know)))
636 bry_str=cff1*(
forces(ng)%svstr(i,jend+1)- &
637 &
forces(ng)%bvstr(i,jend+1))
638 ce=1.0_r8/sqrt(
g*0.5_r8*(
grid(ng)%h(i,jend+1)+ &
639 & zeta(i,jend+1,know)+ &
640 &
grid(ng)%h(i,jend )+ &
641 & zeta(i,jend ,know)))
642# ifdef WEC_VF
643
644 bry_pgr=min(0.5_r8*(
grid(ng)%h(i,jend)+ &
645 & zeta(i,jend,know)), 1.0_r8)*bry_pgr
646 cff=1.0_r8/(0.5_r8*(
grid(ng)%h(i,jend )+ &
647 & zeta(i,jend ,know)+ &
648 &
grid(ng)%h(i,jend+1)+ &
649 & zeta(i,jend+1,know)))
650 cff1=1.5_r8*
pi-
forces(ng)%Dwave(i,jend)- &
651 &
grid(ng)%angler(i,jend)
652 waven=2.0_r8*
pi/max(
forces(ng)%Lwave(i,jend),lwave_min)
653 waveny=waven*sin(cff1)
654 sigma=sqrt(max(
g*waven*tanh(waven/cff),eps))
655 osigma=1.0_r8/sigma
657 &
forces(ng)%Dissip_break(i,jend)*waveny
658 bry_wec=cff1
659# ifdef WEC_ROLLER
660 bry_wec=bry_wec+osigma*
forces(ng)%Dissip_roller(i,jend)* &
661 & waveny
662# endif
663
664# else
665 bry_wec=0.0_r8
666# endif
667 cff2=
grid(ng)%on_v(i,jend+1)*ce
668
669 bry_val=vbar(i,jend,know)+ &
670 & cff2*(bry_pgr+ &
671 & bry_cor+ &
672 & bry_wec+ &
673 & bry_str)
674#else
676#endif
677 cff=1.0_r8/(0.5_r8*(
grid(ng)%h(i,jend )+ &
678 & zeta(i,jend ,know)+ &
679 &
grid(ng)%h(i,jend+1)+ &
680 & zeta(i,jend+1,know)))
682#if defined ATM_PRESS && defined PRESS_COMPENSATE
683 vbar(i,jend+1,kout)=bry_val+ &
684 & ce*(0.5_r8* &
685 & (zeta(i,jend ,know)+ &
686 & zeta(i,jend+1,know)+ &
687 & fac*(
forces(ng)%Pair(i,jend )+ &
688 &
forces(ng)%Pair(i,jend+1)- &
689 & 2.0_r8*oneatm))- &
691#else
692 vbar(i,jend+1,kout)=bry_val+ &
693 & ce*(0.5_r8*(zeta(i,jend ,know)+ &
694 & zeta(i,jend+1,know))- &
696#endif
697#ifdef MASKING
698 vbar(i,jend+1,kout)=vbar(i,jend+1,kout)* &
699 &
grid(ng)%vmask(i,jend+1)
700#endif
701 END IF
702 END DO
703
704
705
707 DO i=istr,iend
709#if defined SSH_TIDES && !defined UV_TIDES
711 bry_pgr=-
g*(
boundary(ng)%zeta_north(i)- &
712 & zeta(i,jend,know))* &
713 & 0.5_r8*
grid(ng)%pn(i,jend)
714 ELSE
715 bry_pgr=-
g*(zeta(i,jend+1,know)- &
716 & zeta(i,jend ,know))* &
717 & 0.5_r8*(
grid(ng)%pn(i,jend )+ &
718 &
grid(ng)%pn(i,jend+1))
719 END IF
720# ifdef UV_COR
721 bry_cor=-0.125_r8*(ubar(i ,jend ,know)+ &
722 & ubar(i+1,jend ,know)+ &
723 & ubar(i ,jend+1,know)+ &
724 & ubar(i+1,jend+1,know))* &
725 & (
grid(ng)%f(i,jend )+ &
726 &
grid(ng)%f(i,jend+1))
727# else
728 bry_cor=0.0_r8
729# endif
730 cff1=1.0_r8/(0.5_r8*(
grid(ng)%h(i,jend )+ &
731 & zeta(i,jend ,know)+ &
732 &
grid(ng)%h(i,jend+1)+ &
733 & zeta(i,jend+1,know)))
734 bry_str=cff1*(
forces(ng)%svstr(i,jend+1)- &
735 &
forces(ng)%bvstr(i,jend+1))
736 ce=1.0_r8/sqrt(
g*0.5_r8*(
grid(ng)%h(i,jend+1)+ &
737 & zeta(i,jend+1,know)+ &
738 &
grid(ng)%h(i,jend )+ &
739 & zeta(i,jend ,know)))
740 cff2=
grid(ng)%on_v(i,jend+1)*ce
741
742 bry_val=vbar(i,jend,know)+ &
743 & cff2*(bry_pgr+ &
744 & bry_cor+ &
745 & bry_str)
746#else
748#endif
749#ifdef WET_DRY
750 cff=0.5_r8*(
grid(ng)%h(i,jend )+ &
751 & zeta(i,jend ,know)+ &
752 &
grid(ng)%h(i,jend+1)+ &
753 & zeta(i,jend+1,know))
754#else
755 cff=0.5_r8*(
grid(ng)%h(i,jend )+ &
756 &
grid(ng)%h(i,jend+1))
757#endif
759 ce=dt2d*cff1*cff*0.5_r8*(
grid(ng)%pn(i,jend )+ &
760 &
grid(ng)%pn(i,jend+1))
761 ze=(0.5_r8+ce)*zeta(i,jend ,know)+ &
762 & (0.5_r8-ce)*zeta(i,jend+1,know)
764 cff2=(1.0_r8-
co/ce)**2
765 cff3=zeta(i,jend,kout)+ &
766 & ce*zeta(i,jend+1,know)- &
767 & (1.0_r8+ce)*zeta(i,jend,know)
768 ze=ze+cff2*cff3
769 END IF
770 vbar(i,jend+1,kout)=0.5_r8* &
771 & ((1.0_r8-ce)*vbar(i,jend+1,know)+ &
772 & ce*vbar(i,jend,know)+ &
773 & bry_val+ &
774 & cff1*(ze-
boundary(ng)%zeta_north(i)))
775#ifdef MASKING
776 vbar(i,jend+1,kout)=vbar(i,jend+1,kout)* &
777 &
grid(ng)%vmask(i,jend+1)
778#endif
779 END IF
780 END DO
781
782
783
785 DO i=istr,iend
787 vbar(i,jend+1,kout)=
boundary(ng)%vbar_north(i)
788#ifdef MASKING
789 vbar(i,jend+1,kout)=vbar(i,jend+1,kout)* &
790 &
grid(ng)%vmask(i,jend+1)
791#endif
792 END IF
793 END DO
794
795
796
798 DO i=istr,iend
800 vbar(i,jend+1,kout)=vbar(i,jend,kout)
801#ifdef MASKING
802 vbar(i,jend+1,kout)=vbar(i,jend+1,kout)* &
803 &
grid(ng)%vmask(i,jend+1)
804#endif
805 END IF
806 END DO
807
808
809
811 DO i=istr,iend
814 bry_pgr=-
g*(
boundary(ng)%zeta_north(i)- &
815 & zeta(i,jend,know))* &
816 & 0.5_r8*
grid(ng)%pn(i,jend)
817 ELSE
818 bry_pgr=-
g*(zeta(i,jend+1,know)- &
819 & zeta(i,jend ,know))* &
820 & 0.5_r8*(
grid(ng)%pn(i,jend )+ &
821 &
grid(ng)%pn(i,jend+1))
822 END IF
823#ifdef UV_COR
824 bry_cor=-0.125_r8*(ubar(i ,jend ,know)+ &
825 & ubar(i+1,jend ,know)+ &
826 & ubar(i ,jend+1,know)+ &
827 & ubar(i+1,jend+1,know))* &
828 & (
grid(ng)%f(i,jend )+ &
829 &
grid(ng)%f(i,jend+1))
830#else
831 bry_cor=0.0_r8
832#endif
833 cff=1.0_r8/(0.5_r8*(
grid(ng)%h(i,jend )+ &
834 & zeta(i,jend ,know)+ &
835 &
grid(ng)%h(i,jend+1)+ &
836 & zeta(i,jend+1,know)))
837 bry_str=cff*(
forces(ng)%svstr(i,jend+1)- &
838 &
forces(ng)%bvstr(i,jend+1))
839 vbar(i,jend+1,kout)=vbar(i,jend+1,know)+ &
840 & dt2d*(bry_pgr+ &
841 & bry_cor+ &
842 & bry_str)
843#ifdef MASKING
844 vbar(i,jend+1,kout)=vbar(i,jend+1,kout)* &
845 &
grid(ng)%vmask(i,jend+1)
846#endif
847 END IF
848 END DO
849
850
851
853 DO i=istr,iend
855 vbar(i,jend+1,kout)=0.0_r8
856 END IF
857 END DO
858
859#ifdef NESTING
860
861
862
863
869 & (rg.eq.ng).and.(
dxmax(dg).gt.
dxmax(rg)))
THEN
872 DO i=istr,iend
876 cff=0.5_r8*
grid(ng)%om_v(i,jend+1)* &
877 & (
grid(ng)%h(i,jend+1)+zeta(i,jend+1,kout)+ &
878 &
grid(ng)%h(i,jend )+zeta(i,jend ,kout))
880 bry_val=cff1*
refined(cr)%V2d_flux(1,m,tnew)/cff
881# ifdef WEC
882 bry_val=bry_val-
ocean(ng)%vbar_stokes(i,jend+1)
883# endif
884# ifdef MASKING
885 bry_val=bry_val*
grid(ng)%vmask(i,jend+1)
886# endif
887# ifdef NESTING_DEBUG
889# endif
890 vbar(i,jend+1,kout)=bry_val
891 END DO
892 END IF
893 END DO
894#endif
895 END IF
896 END IF
897
898
899
900
901
902 IF (
domain(ng)%Western_Edge(tile))
THEN
903
904
905
907 DO j=jstrv-1,jend
908 grad(istr-1,j)=vbar(istr-1,j+1,know)- &
909 & vbar(istr-1,j ,know)
910 grad(istr ,j)=vbar(istr ,j+1,know)- &
911 & vbar(istr ,j ,know)
912 END DO
913 DO j=jstrv,jend
915 dvdt=vbar(istr,j,know)-vbar(istr ,j,kout)
916 dvdx=vbar(istr,j,kout)-vbar(istr+1,j,kout)
917
920 obc_out=0.5_r8* &
921 & (
clima(ng)%M2nudgcof(istr-1,j-1)+ &
922 &
clima(ng)%M2nudgcof(istr-1,j ))
923 obc_in =
obcfac(ng)*obc_out
924 ELSE
927 END IF
928 IF ((dvdt*dvdx).lt.0.0_r8) THEN
929 tau=obc_in
930 ELSE
931 tau=obc_out
932 END IF
933#ifdef IMPLICIT_NUDGING
934 IF (tau.gt.0.0_r8) tau=1.0_r8/tau
935#else
936 tau=tau*dt2d
937#endif
938 END IF
939
940 IF ((dvdt*dvdx).lt.0.0_r8) dvdt=0.0_r8
941 IF ((dvdt*(grad(istr,j-1)+ &
942 & grad(istr,j ))).gt.0.0_r8) THEN
943 dvde=grad(istr,j-1)
944 ELSE
945 dvde=grad(istr,j )
946 END IF
947 cff=max(dvdx*dvdx+dvde*dvde,eps)
948 cx=dvdt*dvdx
949#ifdef RADIATION_2D
950 ce=min(cff,max(dvdt*dvde,-cff))
951#else
952 ce=0.0_r8
953#endif
954#if defined CELERITY_WRITE && defined FORWARD_WRITE
958#endif
959 vbar(istr-1,j,kout)=(cff*vbar(istr-1,j,know)+ &
960 & cx *vbar(istr ,j,kout)- &
961 & max(ce,0.0_r8)*grad(istr-1,j-1)- &
962 & min(ce,0.0_r8)*grad(istr-1,j ))/ &
963 & (cff+cx)
964
966#ifdef IMPLICIT_NUDGING
967 phi=
dt(ng)/(tau+
dt(ng))
968 vbar(istr-1,j,kout)=(1.0_r8-phi)*vbar(istr-1,j,kout)+ &
970#else
971 vbar(istr-1,j,kout)=vbar(istr-1,j,kout)+ &
973 & vbar(istr-1,j,know))
974#endif
975 END IF
976#ifdef MASKING
977 vbar(istr-1,j,kout)=vbar(istr-1,j,kout)* &
978 &
grid(ng)%vmask(istr-1,j)
979#endif
980 END IF
981 END DO
982
983
984
988 DO j=jstrv,jend
990 cff=dt2d*0.5_r8*(
grid(ng)%pm(istr,j-1)+ &
991 &
grid(ng)%pm(istr,j ))
992 cff1=sqrt(
g*0.5_r8*(
grid(ng)%h(istr,j-1)+ &
993 & zeta(istr,j-1,know)+ &
994 &
grid(ng)%h(istr,j )+ &
995 & zeta(istr,j ,know)))
996 cx=cff*cff1
997 cff2=1.0_r8/(1.0_r8+cx)
998 vbar(istr-1,j,kout)=cff2*(vbar(istr-1,j,know)+ &
999 & cx*vbar(istr,j,kout))
1000#ifdef MASKING
1001 vbar(istr-1,j,kout)=vbar(istr-1,j,kout)* &
1002 &
grid(ng)%vmask(istr-1,j)
1003#endif
1004 END IF
1005 END DO
1006
1007
1008
1010 DO j=jstrv,jend
1012 vbar(istr-1,j,kout)=
boundary(ng)%vbar_west(j)
1013#ifdef MASKING
1014 vbar(istr-1,j,kout)=vbar(istr-1,j,kout)* &
1015 &
grid(ng)%vmask(istr-1,j)
1016#endif
1017 END IF
1018 END DO
1019
1020
1021
1023 DO j=jstrv,jend
1025 vbar(istr-1,j,kout)=vbar(istr,j,kout)
1026#ifdef MASKING
1027 vbar(istr-1,j,kout)=vbar(istr-1,j,kout)* &
1028 &
grid(ng)%vmask(istr-1,j)
1029#endif
1030 END IF
1031 END DO
1032
1033
1034
1035
1038 jmin=jstrv
1039 jmax=jend
1040 ELSE
1041 jmin=jstr
1042 jmax=jendr
1043 END IF
1044 DO j=jmin,jmax
1046 vbar(istr-1,j,kout)=
gamma2(ng)*vbar(istr,j,kout)
1047#ifdef MASKING
1048 vbar(istr-1,j,kout)=vbar(istr-1,j,kout)* &
1049 &
grid(ng)%vmask(istr-1,j)
1050#endif
1051 END IF
1052 END DO
1053 END IF
1054 END IF
1055
1056
1057
1058
1059
1060 IF (
domain(ng)%Eastern_Edge(tile))
THEN
1061
1062
1063
1065 DO j=jstrv-1,jend
1066 grad(iend ,j)=vbar(iend ,j+1,know)- &
1067 & vbar(iend ,j ,know)
1068 grad(iend+1,j)=vbar(iend+1,j+1,know)- &
1069 & vbar(iend+1,j ,know)
1070 END DO
1071 DO j=jstrv,jend
1073 dvdt=vbar(iend,j,know)-vbar(iend ,j,kout)
1074 dvdx=vbar(iend,j,kout)-vbar(iend-1,j,kout)
1075
1078 obc_out=0.5_r8* &
1079 & (
clima(ng)%M2nudgcof(iend+1,j-1)+ &
1080 &
clima(ng)%M2nudgcof(iend+1,j ))
1081 obc_in =
obcfac(ng)*obc_out
1082 ELSE
1085 END IF
1086 IF ((dvdt*dvdx).lt.0.0_r8) THEN
1087 tau=obc_in
1088 ELSE
1089 tau=obc_out
1090 END IF
1091 tau=tau*dt2d
1092 END IF
1093
1094 IF ((dvdt*dvdx).lt.0.0_r8) dvdt=0.0_r8
1095 IF ((dvdt*(grad(iend,j-1)+ &
1096 & grad(iend,j ))).gt.0.0_r8) THEN
1097 dvde=grad(iend,j-1)
1098 ELSE
1099 dvde=grad(iend,j )
1100 END IF
1101 cff=max(dvdx*dvdx+dvde*dvde,eps)
1102 cx=dvdt*dvdx
1103#ifdef RADIATION_2D
1104 ce=min(cff,max(dvdt*dvde,-cff))
1105#else
1106 ce=0.0_r8
1107#endif
1108#if defined CELERITY_WRITE && defined FORWARD_WRITE
1112#endif
1113 vbar(iend+1,j,kout)=(cff*vbar(iend+1,j,know)+ &
1114 & cx *vbar(iend ,j,kout)- &
1115 & max(ce,0.0_r8)*grad(iend+1,j-1)- &
1116 & min(ce,0.0_r8)*grad(iend+1,j ))/ &
1117 & (cff+cx)
1118
1120 vbar(iend+1,j,kout)=vbar(iend+1,j,kout)+ &
1121 & tau*(
boundary(ng)%vbar_east(j)- &
1122 & vbar(iend+1,j,know))
1123 END IF
1124#ifdef MASKING
1125 vbar(iend+1,j,kout)=vbar(iend+1,j,kout)* &
1126 &
grid(ng)%vmask(iend+1,j)
1127#endif
1128 END IF
1129 END DO
1130
1131
1132
1136 DO j=jstrv,jend
1138 cff=dt2d*0.5_r8*(
grid(ng)%pm(iend,j-1)+ &
1139 &
grid(ng)%pm(iend,j ))
1140 cff1=sqrt(
g*0.5_r8*(
grid(ng)%h(iend,j-1)+ &
1141 & zeta(iend,j-1,know)+ &
1142 &
grid(ng)%h(iend,j )+ &
1143 & zeta(iend,j ,know)))
1144 cx=cff*cff1
1145 cff2=1.0_r8/(1.0_r8+cx)
1146 vbar(iend+1,j,kout)=cff2*(vbar(iend+1,j,know)+ &
1147 & cx*vbar(iend,j,kout))
1148#ifdef MASKING
1149 vbar(iend+1,j,kout)=vbar(iend+1,j,kout)* &
1150 &
grid(ng)%vmask(iend+1,j)
1151#endif
1152 END IF
1153 END DO
1154
1155
1156
1158 DO j=jstrv,jend
1160 vbar(iend+1,j,kout)=
boundary(ng)%vbar_east(j)
1161#ifdef MASKING
1162 vbar(iend+1,j,kout)=vbar(iend+1,j,kout)* &
1163 &
grid(ng)%vmask(iend+1,j)
1164#endif
1165 END IF
1166 END DO
1167
1168
1169
1171 DO j=jstrv,jend
1173 vbar(iend+1,j,kout)=vbar(iend,j,kout)
1174#ifdef MASKING
1175 vbar(iend+1,j,kout)=vbar(iend+1,j,kout)* &
1176 &
grid(ng)%vmask(iend+1,j)
1177#endif
1178 END IF
1179 END DO
1180
1181
1182
1183
1186 jmin=jstrv
1187 jmax=jend
1188 ELSE
1189 jmin=jstr
1190 jmax=jendr
1191 END IF
1192 DO j=jmin,jmax
1194 vbar(iend+1,j,kout)=
gamma2(ng)*vbar(iend,j,kout)
1195#ifdef MASKING
1196 vbar(iend+1,j,kout)=vbar(iend+1,j,kout)* &
1197 &
grid(ng)%vmask(iend+1,j)
1198#endif
1199 END IF
1200 END DO
1201 END IF
1202 END IF
1203
1204
1205
1206
1207
1209 IF (
domain(ng)%SouthWest_Corner(tile))
THEN
1212 vbar(istr-1,jstr,kout)=0.5_r8*(vbar(istr ,jstr ,kout)+ &
1213 & vbar(istr-1,jstr+1,kout))
1214 END IF
1215 END IF
1216 IF (
domain(ng)%SouthEast_Corner(tile))
THEN
1219 vbar(iend+1,jstr,kout)=0.5_r8*(vbar(iend ,jstr ,kout)+ &
1220 & vbar(iend+1,jstr+1,kout))
1221 END IF
1222 END IF
1223 IF (
domain(ng)%NorthWest_Corner(tile))
THEN
1226 vbar(istr-1,jend+1,kout)=0.5_r8*(vbar(istr-1,jend ,kout)+ &
1227 & vbar(istr ,jend+1,kout))
1228 END IF
1229 END IF
1230 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
1233 vbar(iend+1,jend+1,kout)=0.5_r8*(vbar(iend+1,jend ,kout)+ &
1234 & vbar(iend ,jend+1,kout))
1235 END IF
1236 END IF
1237 END IF
1238
1239#if defined WET_DRY
1240
1241
1242
1243
1244
1246 IF (
domain(ng)%Western_Edge(tile))
THEN
1247 DO j=jstrv,jend
1250 cff1=abs(abs(
grid(ng)%vmask_wet(istr-1,j))-1.0_r8)
1251 cff2=0.5_r8+dsign(0.5_r8,vbar(istr-1,j,kout))* &
1252 &
grid(ng)%vmask_wet(istr-1,j)
1253 cff=0.5_r8*
grid(ng)%vmask_wet(istr-1,j)*cff1+ &
1254 & cff2*(1.0_r8-cff1)
1255 vbar(istr,j,kout)=vbar(istr,j,kout)*cff
1256 END IF
1257 END DO
1258 END IF
1259 IF (
domain(ng)%Eastern_Edge(tile))
THEN
1260 DO j=jstrv,jend
1263 cff1=abs(abs(
grid(ng)%vmask_wet(iend+1,j))-1.0_r8)
1264 cff2=0.5_r8+dsign(0.5_r8,vbar(iend+1,j,kout))* &
1265 &
grid(ng)%vmask_wet(iend+1,j)
1266 cff=0.5_r8*
grid(ng)%vmask_wet(iend+1,j)*cff1+ &
1267 & cff2*(1.0_r8-cff1)
1268 vbar(iend+1,j,kout)=vbar(iend+1,j,kout)*cff
1269 END IF
1270 END DO
1271 END IF
1272 END IF
1273
1275 IF (
domain(ng)%Southern_Edge(tile))
THEN
1276 DO i=istr,iend
1279 cff1=abs(abs(
grid(ng)%vmask_wet(i,jstr))-1.0_r8)
1280 cff2=0.5_r8+dsign(0.5_r8,vbar(i,jstr,kout))* &
1281 &
grid(ng)%vmask_wet(i,jstr)
1282 cff=0.5_r8*
grid(ng)%vmask_wet(i,jstr)*cff1+ &
1283 & cff2*(1.0_r8-cff1)
1284 vbar(i,jstr,kout)=vbar(i,jstr,kout)*cff
1285 END IF
1286 END DO
1287 END IF
1288 IF (
domain(ng)%Northern_Edge(tile))
THEN
1289 DO i=istr,iend
1292 cff1=abs(abs(
grid(ng)%vmask_wet(i,jend+1))-1.0_r8)
1293 cff2=0.5_r8+dsign(0.5_r8,vbar(i,jend+1,kout))* &
1294 &
grid(ng)%vmask_wet(i,jend+1)
1295 cff=0.5_r8*
grid(ng)%vmask_wet(i,jend+1)*cff1+ &
1296 & cff2*(1.0_r8-cff1)
1297 vbar(i,jend+1,kout)=vbar(i,jend+1,kout)*cff
1298 END IF
1299 END DO
1300 END IF
1301 END IF
1302
1304 IF (
domain(ng)%SouthWest_Corner(tile))
THEN
1309 cff1=abs(abs(
grid(ng)%vmask_wet(istr-1,jstr))-1.0_r8)
1310 cff2=0.5_r8+dsign(0.5_r8,vbar(istr-1,jstr,kout))* &
1311 &
grid(ng)%vmask_wet(istr-1,jstr)
1312 cff=0.5_r8*
grid(ng)%vmask_wet(istr-1,jstr)*cff1+ &
1313 & cff2*(1.0_r8-cff1)
1314 vbar(istr-1,jstr,kout)=vbar(istr-1,jstr,kout)*cff
1315 END IF
1316 END IF
1317 IF (
domain(ng)%SouthEast_Corner(tile))
THEN
1322 cff1=abs(abs(
grid(ng)%vmask_wet(iend+1,jstr))-1.0_r8)
1323 cff2=0.5_r8+dsign(0.5_r8,vbar(iend+1,jstr,kout))* &
1324 &
grid(ng)%vmask_wet(iend+1,jstr)
1325 cff=0.5_r8*
grid(ng)%vmask_wet(iend+1,jstr)*cff1+ &
1326 & cff2*(1.0_r8-cff1)
1327 vbar(iend+1,jstr,kout)=vbar(iend+1,jstr,kout)*cff
1328 END IF
1329 END IF
1330 IF (
domain(ng)%NorthWest_Corner(tile))
THEN
1335 cff1=abs(abs(
grid(ng)%vmask_wet(istr-1,jend+1))-1.0_r8)
1336 cff2=0.5_r8+dsign(0.5_r8,vbar(istr-1,jend+1,kout))* &
1337 &
grid(ng)%vmask_wet(istr-1,jend+1)
1338 cff=0.5_r8*
grid(ng)%vmask_wet(istr-1,jend+1)*cff1+ &
1339 & cff2*(1.0_r8-cff1)
1340 vbar(istr-1,jend+1,kout)=vbar(istr-1,jend+1,kout)*cff
1341 END IF
1342 END IF
1343 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
1348 cff1=abs(abs(
grid(ng)%vmask_wet(iend+1,jend+1))-1.0_r8)
1349 cff2=0.5_r8+dsign(0.5_r8,vbar(iend+1,jend+1,kout))* &
1350 &
grid(ng)%vmask_wet(iend+1,jend+1)
1351 cff=0.5_r8*
grid(ng)%vmask_wet(iend+1,jend+1)*cff1+ &
1352 & cff2*(1.0_r8-cff1)
1353 vbar(iend+1,jend+1,kout)=vbar(iend+1,jend+1,kout)*cff
1354 END IF
1355 END IF
1356 END IF
1357#endif
1358
1359 RETURN
1360
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_ngc), dimension(:), allocatable vcontact
type(t_bcp), dimension(:,:), allocatable bry_contact
integer, dimension(:), allocatable rollingindex
type(t_refined), dimension(:), allocatable refined
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
real(r8), dimension(:), allocatable wec_alpha