55
56
63
64
65
66 integer, intent(in) :: ng, tile
67 integer, intent(in) :: LBi, UBi, LBj, UBj, UBk
68 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
69 integer, intent(in) :: nstp, nout
70
71# ifdef ASSUMED_SHAPE
72 real(r8), intent(inout) :: u(LBi:,LBj:,:,:)
73# else
74 real(r8), intent(inout) :: u(LBi:UBi,LBj:UBj,UBk,2)
75# endif
76
77
78
79 integer :: Imin, Imax
80 integer :: i, j, k
81
82 real(r8), parameter :: eps = 1.0e-20_r8
83
84 real(r8) :: Ce, Cx, cff, dUde, dUdt, dUdx
85 real(r8) :: obc_in, obc_out, phi, tau
86
87 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: grad
88
89# include "set_bounds.h"
90
91
92
93
94
95 IF (
domain(ng)%Western_Edge(tile))
THEN
96
97
98
101 DO j=jstr,jend+1
102 grad(istr ,j)=u(istr ,j ,k,nstp)- &
103 & u(istr ,j-1,k,nstp)
104 grad(istr+1,j)=u(istr+1,j ,k,nstp)- &
105 & u(istr+1,j-1,k,nstp)
106 END DO
107 DO j=jstr,jend
109 dudt=u(istr+1,j,k,nstp)-u(istr+1,j,k,nout)
110 dudx=u(istr+1,j,k,nout)-u(istr+2,j,k,nout)
111
114 obc_out=0.5_r8* &
115 & (
clima(ng)%M3nudgcof(istr-1,j,k)+ &
116 &
clima(ng)%M3nudgcof(istr ,j,k))
117 obc_in =
obcfac(ng)*obc_out
118 ELSE
121 END IF
122 IF ((dudt*dudx).lt.0.0_r8) THEN
123 tau=obc_in
124 ELSE
125 tau=obc_out
126 END IF
127# ifdef IMPLICIT_NUDGING
128 IF (tau.gt.0.0_r8) tau=1.0_r8/tau
129# else
131# endif
132 END IF
133
134 IF ((dudt*dudx).lt.0.0_r8) dudt=0.0_r8
135 IF ((dudt*(grad(istr+1,j )+ &
136 & grad(istr+1,j+1))).gt.0.0_r8) THEN
137 dude=grad(istr+1,j )
138 ELSE
139 dude=grad(istr+1,j+1)
140 END IF
141 cff=max(dudx*dudx+dude*dude,eps)
142 cx=dudt*dudx
143# ifdef RADIATION_2D
144 ce=min(cff,max(dudt*dude,-cff))
145# else
146 ce=0.0_r8
147# endif
148# if defined CELERITY_WRITE && defined FORWARD_WRITE
152# endif
153 u(istr,j,k,nout)=(cff*u(istr ,j,k,nstp)+ &
154 & cx *u(istr+1,j,k,nout)- &
155 & max(ce,0.0_r8)*grad(istr,j )- &
156 & min(ce,0.0_r8)*grad(istr,j+1))/ &
157 & (cff+cx)
158
160# ifdef IMPLICIT_NUDGING
161 phi=
dt(ng)/(tau+
dt(ng))
162 u(istr,j,k,nout)=(1.0_r8-phi)*u(istr,j,k,nout)+ &
164# else
165 u(istr,j,k,nout)=u(istr,j,k,nout)+ &
167 & u(istr,j,k,nstp))
168# endif
169 END IF
170# ifdef MASKING
171 u(istr,j,k,nout)=u(istr,j,k,nout)* &
172 &
grid(ng)%umask(istr,j)
173# endif
174# ifdef WET_DRY
175 u(istr,j,k,nout)=u(istr,j,k,nout)* &
176 &
grid(ng)%umask_wet(istr,j)
177# endif
178 END IF
179 END DO
180 END DO
181
182
183
186 DO j=jstr,jend
188 u(istr,j,k,nout)=
boundary(ng)%u_west(j,k)
189# ifdef MASKING
190 u(istr,j,k,nout)=u(istr,j,k,nout)* &
191 &
grid(ng)%umask(istr,j)
192# endif
193# ifdef WET_DRY
194 u(istr,j,k,nout)=u(istr,j,k,nout)* &
195 &
grid(ng)%umask_wet(istr,j)
196# endif
197 END IF
198 END DO
199 END DO
200
201
202
205 DO j=jstr,jend
207 u(istr,j,k,nout)=u(istr+1,j,k,nout)
208# ifdef MASKING
209 u(istr,j,k,nout)=u(istr,j,k,nout)* &
210 &
grid(ng)%umask(istr,j)
211# endif
212# ifdef WET_DRY
213 u(istr,j,k,nout)=u(istr,j,k,nout)* &
214 &
grid(ng)%umask_wet(istr,j)
215# endif
216 END IF
217 END DO
218 END DO
219
220
221
224 DO j=jstr,jend
226 u(istr,j,k,nout)=0.0_r8
227 END IF
228 END DO
229 END DO
230 END IF
231 END IF
232
233
234
235
236
237 IF (
domain(ng)%Eastern_Edge(tile))
THEN
238
239
240
243 DO j=jstr,jend+1
244 grad(iend ,j)=u(iend ,j ,k,nstp)- &
245 & u(iend ,j-1,k,nstp)
246 grad(iend+1,j)=u(iend+1,j ,k,nstp)- &
247 & u(iend+1,j-1,k,nstp)
248 END DO
249 DO j=jstr,jend
251 dudt=u(iend,j,k,nstp)-u(iend ,j,k,nout)
252 dudx=u(iend,j,k,nout)-u(iend-1,j,k,nout)
253
256 obc_out=0.5_r8* &
257 & (
clima(ng)%M3nudgcof(iend ,j,k)+ &
258 &
clima(ng)%M3nudgcof(iend+1,j,k))
259 obc_in =
obcfac(ng)*obc_out
260 ELSE
263 END IF
264 IF ((dudt*dudx).lt.0.0_r8) THEN
265 tau=obc_in
266 ELSE
267 tau=obc_out
268 END IF
269# ifdef IMPLICIT_NUDGING
270 IF (tau.gt.0.0_r8) tau=1.0_r8/tau
271# else
273# endif
274 END IF
275
276 IF ((dudt*dudx).lt.0.0_r8) dudt=0.0_r8
277 IF ((dudt*(grad(iend,j )+ &
278 & grad(iend,j+1))).gt.0.0_r8) THEN
279 dude=grad(iend,j )
280 ELSE
281 dude=grad(iend,j+1)
282 END IF
283 cff=max(dudx*dudx+dude*dude,eps)
284 cx=dudt*dudx
285# ifdef RADIATION_2D
286 ce=min(cff,max(dudt*dude,-cff))
287# else
288 ce=0.0_r8
289# endif
290# if defined CELERITY_WRITE && defined FORWARD_WRITE
294# endif
295 u(iend+1,j,k,nout)=(cff*u(iend+1,j,k,nstp)+ &
296 & cx *u(iend ,j,k,nout)- &
297 & max(ce,0.0_r8)*grad(iend+1,j )- &
298 & min(ce,0.0_r8)*grad(iend+1,j+1))/ &
299 & (cff+cx)
300
302# ifdef IMPLICIT_NUDGING
303 phi=
dt(ng)/(tau+
dt(ng))
304 u(iend+1,j,k,nout)=(1.0_r8-phi)*u(iend+1,j,k,nout)+ &
306# else
307 u(iend+1,j,k,nout)=u(iend+1,j,k,nout)+ &
309 & u(iend+1,j,k,nstp))
310# endif
311 END IF
312# ifdef MASKING
313 u(iend+1,j,k,nout)=u(iend+1,j,k,nout)* &
314 &
grid(ng)%umask(iend+1,j)
315# endif
316# ifdef WET_DRY
317 u(iend+1,j,k,nout)=u(iend+1,j,k,nout)* &
318 &
grid(ng)%umask_wet(iend+1,j)
319# endif
320 END IF
321 END DO
322 END DO
323
324
325
328 DO j=jstr,jend
330 u(iend+1,j,k,nout)=
boundary(ng)%u_east(j,k)
331# ifdef MASKING
332 u(iend+1,j,k,nout)=u(iend+1,j,k,nout)* &
333 &
grid(ng)%umask(iend+1,j)
334# endif
335# ifdef WET_DRY
336 u(iend+1,j,k,nout)=u(iend+1,j,k,nout)* &
337 &
grid(ng)%umask_wet(iend+1,j)
338# endif
339 END IF
340 END DO
341 END DO
342
343
344
347 DO j=jstr,jend
349 u(iend+1,j,k,nout)=u(iend,j,k,nout)
350# ifdef MASKING
351 u(iend+1,j,k,nout)=u(iend+1,j,k,nout)* &
352 &
grid(ng)%umask(iend+1,j)
353# endif
354# ifdef WET_DRY
355 u(iend+1,j,k,nout)=u(iend+1,j,k,nout)* &
356 &
grid(ng)%umask_wet(iend+1,j)
357# endif
358 END IF
359 END DO
360 END DO
361
362
363
366 DO j=jstr,jend
368 u(iend+1,j,k,nout)=0.0_r8
369 END IF
370 END DO
371 END DO
372 END IF
373 END IF
374
375
376
377
378
379 IF (
domain(ng)%Southern_Edge(tile))
THEN
380
381
382
385 DO i=istru-1,iend
386 grad(i,jstr-1)=u(i+1,jstr-1,k,nstp)- &
387 & u(i ,jstr-1,k,nstp)
388 grad(i,jstr )=u(i+1,jstr ,k,nstp)- &
389 & u(i ,jstr ,k,nstp)
390 END DO
391 DO i=istru,iend
393 dudt=u(i,jstr,k,nstp)-u(i,jstr ,k,nout)
394 dude=u(i,jstr,k,nout)-u(i,jstr+1,k,nout)
395
398 obc_out=0.5_r8* &
399 & (
clima(ng)%M3nudgcof(i-1,jstr-1,k)+ &
400 &
clima(ng)%M3nudgcof(i ,jstr-1,k))
401 obc_in =
obcfac(ng)*obc_out
402 ELSE
405 END IF
406 IF ((dudt*dude).lt.0.0_r8) THEN
407 tau=obc_in
408 ELSE
409 tau=obc_out
410 END IF
411# ifdef IMPLICIT_NUDGING
412 IF (tau.gt.0.0_r8) tau=1.0_r8/tau
413# else
415# endif
416 END IF
417
418 IF ((dudt*dude).lt.0.0_r8) dudt=0.0_r8
419 IF ((dudt*(grad(i-1,jstr)+ &
420 & grad(i ,jstr))).gt.0.0_r8) THEN
421 dudx=grad(i-1,jstr)
422 ELSE
423 dudx=grad(i ,jstr)
424 END IF
425 cff=max(dudx*dudx+dude*dude,eps)
426# ifdef RADIATION_2D
427 cx=min(cff,max(dudt*dudx,-cff))
428# else
429 cx=0.0_r8
430# endif
431 ce=dudt*dude
432# if defined CELERITY_WRITE && defined FORWARD_WRITE
436# endif
437 u(i,jstr-1,k,nout)=(cff*u(i,jstr-1,k,nstp)+ &
438 & ce *u(i,jstr ,k,nout)- &
439 & max(cx,0.0_r8)*grad(i-1,jstr-1)- &
440 & min(cx,0.0_r8)*grad(i ,jstr-1))/ &
441 & (cff+ce)
442
444# ifdef IMPLICIT_NUDGING
445 phi=
dt(ng)/(tau+
dt(ng))
446 u(i,jstr-1,k,nout)=(1.0_r8-phi)*u(i,jstr-1,k,nout)+ &
448# else
449 u(i,jstr-1,k,nout)=u(i,jstr-1,k,nout)+ &
451 & u(i,jstr-1,k,nstp))
452# endif
453 END IF
454# ifdef MASKING
455 u(i,jstr-1,k,nout)=u(i,jstr-1,k,nout)* &
456 &
grid(ng)%umask(i,jstr-1)
457# endif
458# ifdef WET_DRY
459 u(i,jstr-1,k,nout)=u(i,jstr-1,k,nout)* &
460 &
grid(ng)%umask_wet(i,jstr-1)
461# endif
462 END IF
463 END DO
464 END DO
465
466
467
470 DO i=istru,iend
472 u(i,jstr-1,k,nout)=
boundary(ng)%u_south(i,k)
473# ifdef MASKING
474 u(i,jstr-1,k,nout)=u(i,jstr-1,k,nout)* &
475 &
grid(ng)%umask(i,jstr-1)
476# endif
477# ifdef WET_DRY
478 u(i,jstr-1,k,nout)=u(i,jstr-1,k,nout)* &
479 &
grid(ng)%umask_wet(i,jstr-1)
480# endif
481 END IF
482 END DO
483 END DO
484
485
486
489 DO i=istru,iend
491 u(i,jstr-1,k,nout)=u(i,jstr,k,nout)
492# ifdef MASKING
493 u(i,jstr-1,k,nout)=u(i,jstr-1,k,nout)* &
494 &
grid(ng)%umask(i,jstr-1)
495# endif
496# ifdef WET_MASK
497 u(i,jstr-1,k,nout)=u(i,jstr-1,k,nout)* &
498 &
grid(ng)%umask_wet(i,jstr-1)
499# endif
500 END IF
501 END DO
502 END DO
503
504
505
506
509 imin=istru
510 imax=iend
511 ELSE
512 imin=istr
513 imax=iendr
514 END IF
516 DO i=imin,imax
518 u(i,jstr-1,k,nout)=
gamma2(ng)*u(i,jstr,k,nout)
519# ifdef MASKING
520 u(i,jstr-1,k,nout)=u(i,jstr-1,k,nout)* &
521 &
grid(ng)%umask(i,jstr-1)
522# endif
523# ifdef WET_DRY
524 u(i,jstr-1,k,nout)=u(i,jstr-1,k,nout)* &
525 &
grid(ng)%umask_wet(i,jstr-1)
526# endif
527 END IF
528 END DO
529 END DO
530 END IF
531 END IF
532
533
534
535
536
537 IF (
domain(ng)%Northern_Edge(tile))
THEN
538
539
540
543 DO i=istru-1,iend
544 grad(i,jend )=u(i+1,jend ,k,nstp)- &
545 & u(i ,jend ,k,nstp)
546 grad(i,jend+1)=u(i+1,jend+1,k,nstp)- &
547 & u(i ,jend+1,k,nstp)
548 END DO
549 DO i=istru,iend
551 dudt=u(i,jend,k,nstp)-u(i,jend ,k,nout)
552 dude=u(i,jend,k,nout)-u(i,jend-1,k,nout)
553
556 obc_out=0.5_r8* &
557 & (
clima(ng)%M3nudgcof(i-1,jend+1,k)+ &
558 &
clima(ng)%M3nudgcof(i ,jend+1,k))
559 obc_in =
obcfac(ng)*obc_out
560 ELSE
563 END IF
564 IF ((dudt*dude).lt.0.0_r8) THEN
565 tau=obc_in
566 ELSE
567 tau=obc_out
568 END IF
569# ifdef IMPLICIT_NUDGING
570 IF (tau.gt.0.0_r8) tau=1.0_r8/tau
571# else
573# endif
574 END IF
575
576 IF ((dudt*dude).lt.0.0_r8) dudt=0.0_r8
577 IF ((dudt*(grad(i-1,jend)+ &
578 & grad(i ,jend))).gt.0.0_r8) THEN
579 dudx=grad(i-1,jend)
580 ELSE
581 dudx=grad(i ,jend)
582 END IF
583 cff=max(dudx*dudx+dude*dude,eps)
584# ifdef RADIATION_2D
585 cx=min(cff,max(dudt*dudx,-cff))
586# else
587 cx=0.0_r8
588# endif
589 ce=dudt*dude
590# if defined CELERITY_WRITE && defined FORWARD_WRITE
594# endif
595 u(i,jend+1,k,nout)=(cff*u(i,jend+1,k,nstp)+ &
596 & ce *u(i,jend ,k,nout)- &
597 & max(cx,0.0_r8)*grad(i-1,jend+1)- &
598 & min(cx,0.0_r8)*grad(i ,jend+1))/ &
599 & (cff+ce)
600
602# ifdef IMPLICIT_NUDGING
603 phi=
dt(ng)/(tau+
dt(ng))
604 u(i,jend+1,k,nout)=(1.0_r8-phi)*u(i,jend+1,k,nout)+ &
606# else
607 u(i,jend+1,k,nout)=u(i,jend+1,k,nout)+ &
609 & u(i,jend+1,k,nstp))
610# endif
611 END IF
612# ifdef MASKING
613 u(i,jend+1,k,nout)=u(i,jend+1,k,nout)* &
614 &
grid(ng)%umask(i,jend+1)
615# endif
616# ifdef WET_DRY
617 u(i,jend+1,k,nout)=u(i,jend+1,k,nout)* &
618 &
grid(ng)%umask_wet(i,jend+1)
619# endif
620 END IF
621 END DO
622 END DO
623
624
625
628 DO i=istru,iend
630 u(i,jend+1,k,nout)=
boundary(ng)%u_north(i,k)
631# ifdef MASKING
632 u(i,jend+1,k,nout)=u(i,jend+1,k,nout)* &
633 &
grid(ng)%umask(i,jend+1)
634# endif
635# ifdef WET_DRY
636 u(i,jend+1,k,nout)=u(i,jend+1,k,nout)* &
637 &
grid(ng)%umask_wet(i,jend+1)
638# endif
639 END IF
640 END DO
641 END DO
642
643
644
647 DO i=istru,iend
649 u(i,jend+1,k,nout)=u(i,jend,k,nout)
650# ifdef MASKING
651 u(i,jend+1,k,nout)=u(i,jend+1,k,nout)* &
652 &
grid(ng)%umask(i,jend+1)
653# endif
654# ifdef WET_DRY
655 u(i,jend+1,k,nout)=u(i,jend+1,k,nout)* &
656 &
grid(ng)%umask_wet(i,jend+1)
657# endif
658 END IF
659 END DO
660 END DO
661
662
663
664
667 imin=istru
668 imax=iend
669 ELSE
670 imin=istr
671 imax=iendr
672 END IF
674 DO i=imin,imax
676 u(i,jend+1,k,nout)=
gamma2(ng)*u(i,jend,k,nout)
677# ifdef MASKING
678 u(i,jend+1,k,nout)=u(i,jend+1,k,nout)* &
679 &
grid(ng)%umask(i,jend+1)
680# endif
681# ifdef WET_DRY
682 u(i,jend+1,k,nout)=u(i,jend+1,k,nout)* &
683 &
grid(ng)%umask_wet(i,jend+1)
684# endif
685 END IF
686 END DO
687 END DO
688 END IF
689 END IF
690
691
692
693
694
696 IF (
domain(ng)%SouthWest_Corner(tile))
THEN
700 u(istr,jstr-1,k,nout)=0.5_r8*(u(istr+1,jstr-1,k,nout)+ &
701 & u(istr ,jstr ,k,nout))
702 END DO
703 END IF
704 END IF
705 IF (
domain(ng)%SouthEast_Corner(tile))
THEN
709 u(iend+1,jstr-1,k,nout)=0.5_r8*(u(iend ,jstr-1,k,nout)+ &
710 & u(iend+1,jstr ,k,nout))
711 END DO
712 END IF
713 END IF
714 IF (
domain(ng)%NorthWest_Corner(tile))
THEN
718 u(istr,jend+1,k,nout)=0.5_r8*(u(istr ,jend ,k,nout)+ &
719 & u(istr+1,jend+1,k,nout))
720 END DO
721 END IF
722 END IF
723 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
727 u(iend+1,jend+1,k,nout)=0.5_r8*(u(iend+1,jend ,k,nout)+ &
728 & u(iend ,jend+1,k,nout))
729 END DO
730 END IF
731 END IF
732 END IF
733
734 RETURN
type(t_boundary), dimension(:), allocatable boundary
type(t_apply), dimension(:), allocatable lbc_apply
type(t_clima), dimension(:), allocatable clima
type(t_grid), dimension(:), allocatable grid
type(t_lbc), dimension(:,:,:), allocatable lbc
type(t_domain), dimension(:), allocatable domain
real(dp), dimension(:,:), allocatable m3obc_out
real(dp), dimension(:), allocatable dt
logical, dimension(:), allocatable ewperiodic
logical, dimension(:), allocatable nsperiodic
real(dp), dimension(:), allocatable obcfac
real(r8), dimension(:), allocatable gamma2
logical, dimension(:), allocatable lnudgem3clm
integer, parameter isouth
integer, parameter inorth
real(dp), dimension(:,:), allocatable m3obc_in