56
57
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) :: nout, nstp
70
71# ifdef ASSUMED_SHAPE
72 real(r8), intent(inout) :: gls(LBi:,LBj:,0:,:)
73 real(r8), intent(inout) :: tke(LBi:,LBj:,0:,:)
74# else
75 real(r8), intent(inout) :: gls(LBi:UBi,LBj:UBj,0:UBk,3)
76 real(r8), intent(inout) :: tke(LBi:UBi,LBj:UBj,0:UBk,3)
77# endif
78
79
80
81 integer :: i, j, k
82
83 real(r8), parameter :: eps = 1.0e-20_r8
84
85 real(r8) :: Ce, Cx, cff, dKde, dKdt, dKdx
86
87 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: grad
88 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: gradL
89
90# include "set_bounds.h"
91
92
93
94
95
96 IF (
domain(ng)%Western_Edge(tile))
THEN
97
98
99
102 DO j=jstr,jend+1
103 grad(istr-1,j)=tke(istr-1,j ,k,nstp)- &
104 & tke(istr-1,j-1,k,nstp)
105# ifdef MASKING
106 grad(istr-1,j)=grad(istr-1,j)* &
107 &
grid(ng)%vmask(istr-1,j)
108# endif
109 grad(istr ,j)=tke(istr ,j ,k,nstp)- &
110 & tke(istr ,j-1,k,nstp)
111# ifdef MASKING
112 grad(istr ,j)=grad(istr ,j)* &
113 &
grid(ng)%vmask(istr ,j)
114# endif
115 gradl(istr-1,j)=gls(istr-1,j ,k,nstp)- &
116 & gls(istr-1,j-1,k,nstp)
117# ifdef MASKING
118 gradl(istr-1,j)=gradl(istr-1,j)* &
119 &
grid(ng)%vmask(istr-1,j)
120# endif
121 gradl(istr ,j)=gls(istr ,j ,k,nstp)- &
122 & gls(istr ,j-1,k,nstp)
123# ifdef MASKING
124 gradl(istr ,j)=gradl(istr ,j)* &
125 &
grid(ng)%vmask(istr ,j)
126# endif
127 END DO
128 DO j=jstr,jend
130 dkdt=tke(istr,j,k,nstp)-tke(istr ,j,k,nout)
131 dkdx=tke(istr,j,k,nout)-tke(istr+1,j,k,nout)
132 IF ((dkdt*dkdx).lt.0.0_r8) dkdt=0.0_r8
133 IF ((dkdt*(grad(istr,j )+ &
134 & grad(istr,j+1))).gt.0.0_r8) THEN
135 dkde=grad(istr,j )
136 ELSE
137 dkde=grad(istr,j+1)
138 END IF
139 cff=max(dkdx*dkdx+dkde*dkde,eps)
140 cx=dkdt*dkdx
141# ifdef RADIATION_2D
142 ce=min(cff,max(dkdt*dkde,-cff))
143# else
144 ce=0.0_r8
145# endif
146 tke(istr-1,j,k,nout)=(cff*tke(istr-1,j,k,nstp)+ &
147 & cx *tke(istr ,j,k,nout)- &
148 & max(ce,0.0_r8)* &
149 & grad(istr-1,j )- &
150 & min(ce,0.0_r8)* &
151 & grad(istr-1,j+1))/ &
152 & (cff+cx)
153# ifdef MASKING
154 tke(istr-1,j,k,nout)=tke(istr-1,j,k,nout)* &
155 &
grid(ng)%rmask(istr-1,j)
156# endif
157 dkdt=gls(istr,j,k,nstp)-gls(istr ,j,k,nout)
158 dkdx=gls(istr,j,k,nout)-gls(istr+1,j,k,nout)
159 IF ((dkdt*dkdx).lt.0.0_r8) dkdt=0.0_r8
160 IF ((dkdt*(gradl(istr,j )+ &
161 & gradl(istr,j+1))).gt.0.0_r8) THEN
162 dkde=gradl(istr,j )
163 ELSE
164 dkde=gradl(istr,j+1)
165 END IF
166 cff=max(dkdx*dkdx+dkde*dkde,eps)
167 cx=dkdt*dkdx
168# ifdef RADIATION_2D
169 ce=min(cff,max(dkdt*dkde,-cff))
170# else
171 ce=0.0_r8
172# endif
173 gls(istr-1,j,k,nout)=(cff*gls(istr-1,j,k,nstp)+ &
174 & cx *gls(istr ,j,k,nout)- &
175 & max(ce,0.0_r8)* &
176 & gradl(istr-1,j )- &
177 & min(ce,0.0_r8)* &
178 & gradl(istr-1,j+1))/ &
179 & (cff+cx)
180# ifdef MASKING
181 gls(istr-1,j,k,nout)=gls(istr-1,j,k,nout)* &
182 &
grid(ng)%rmask(istr-1,j)
183# endif
184 END IF
185 END DO
186 END DO
187
188
189
192 DO j=jstr,jend
194 tke(istr-1,j,k,nout)=tke(istr,j,k,nout)
195# ifdef MASKING
196 tke(istr-1,j,k,nout)=tke(istr-1,j,k,nout)* &
197 &
grid(ng)%rmask(istr-1,j)
198# endif
199 gls(istr-1,j,k,nout)=gls(istr,j,k,nout)
200# ifdef MASKING
201 gls(istr-1,j,k,nout)=gls(istr-1,j,k,nout)* &
202 &
grid(ng)%rmask(istr-1,j)
203# endif
204 END IF
205 END DO
206 END DO
207
208
209
212 DO j=jstr,jend
214 tke(istr-1,j,k,nout)=tke(istr,j,k,nout)
215# ifdef MASKING
216 tke(istr-1,j,k,nout)=tke(istr-1,j,k,nout)* &
217 &
grid(ng)%rmask(istr-1,j)
218# endif
219 gls(istr-1,j,k,nout)=gls(istr,j,k,nout)
220# ifdef MASKING
221 gls(istr-1,j,k,nout)=gls(istr-1,j,k,nout)* &
222 &
grid(ng)%rmask(istr-1,j)
223# endif
224 END IF
225 END DO
226 END DO
227 END IF
228 END IF
229
230
231
232
233
234 IF (
domain(ng)%Eastern_Edge(tile))
THEN
235
236
237
240 DO j=jstr,jend+1
241 grad(iend ,j)=tke(iend ,j ,k,nstp)- &
242 & tke(iend ,j-1,k,nstp)
243# ifdef MASKING
244 grad(iend ,j)=grad(iend ,j)* &
245 &
grid(ng)%vmask(iend ,j)
246# endif
247 grad(iend+1,j)=tke(iend+1,j ,k,nstp)- &
248 & tke(iend+1,j-1,k,nstp)
249# ifdef MASKING
250 grad(iend+1,j)=grad(iend+1,j)* &
251 &
grid(ng)%vmask(iend+1,j)
252# endif
253 gradl(iend ,j)=gls(iend ,j ,k,nstp)- &
254 & gls(iend ,j-1,k,nstp)
255# ifdef MASKING
256 gradl(iend ,j)=gradl(iend ,j)* &
257 &
grid(ng)%vmask(iend ,j)
258# endif
259 gradl(iend+1,j)=gls(iend+1,j ,k,nstp)- &
260 & gls(iend+1,j-1,k,nstp)
261# ifdef MASKING
262 gradl(iend+1,j)=gradl(iend+1,j)* &
263 &
grid(ng)%vmask(iend+1,j)
264# endif
265 END DO
266 DO j=jstr,jend
268 dkdt=tke(iend,j,k,nstp)-tke(iend ,j,k,nout)
269 dkdx=tke(iend,j,k,nout)-tke(iend-1,j,k,nout)
270 IF ((dkdt*dkdx).lt.0.0_r8) dkdt=0.0_r8
271 IF ((dkdt*(grad(iend,j )+ &
272 & grad(iend,j+1))).gt.0.0_r8) THEN
273 dkde=grad(iend,j )
274 ELSE
275 dkde=grad(iend,j+1)
276 END IF
277 cff=max(dkdx*dkdx+dkde*dkde,eps)
278 cx=dkdt*dkdx
279# ifdef RADIATION_2D
280 ce=min(cff,max(dkdt*dkde,-cff))
281# else
282 ce=0.0_r8
283# endif
284 tke(iend+1,j,k,nout)=(cff*tke(iend+1,j,k,nstp)+ &
285 & cx *tke(iend ,j,k,nout)- &
286 & max(ce,0.0_r8)* &
287 & grad(iend+1,j )- &
288 & min(ce,0.0_r8)* &
289 & grad(iend+1,j+1))/ &
290 & (cff+cx)
291# ifdef MASKING
292 tke(iend+1,j,k,nout)=tke(iend+1,j,k,nout)* &
293 &
grid(ng)%rmask(iend+1,j)
294# endif
295 dkdt=gls(iend,j,k,nstp)-gls(iend ,j,k,nout)
296 dkdx=gls(iend,j,k,nout)-gls(iend-1,j,k,nout)
297 IF ((dkdt*dkdx).lt.0.0_r8) dkdt=0.0_r8
298 IF ((dkdt*(gradl(iend,j )+ &
299 & gradl(iend,j+1))).gt.0.0_r8) THEN
300 dkde=gradl(iend,j )
301 ELSE
302 dkde=gradl(iend,j+1)
303 END IF
304 cff=max(dkdx*dkdx+dkde*dkde,eps)
305 cx=dkdt*dkdx
306# ifdef RADIATION_2D
307 ce=min(cff,max(dkdt*dkde,-cff))
308# else
309 ce=0.0_r8
310# endif
311 gls(iend+1,j,k,nout)=(cff*gls(iend+1,j,k,nstp)+ &
312 & cx *gls(iend ,j,k,nout)- &
313 & max(ce,0.0_r8)* &
314 & gradl(iend+1,j )- &
315 & min(ce,0.0_r8)* &
316 & gradl(iend+1,j+1))/ &
317 & (cff+cx)
318# ifdef MASKING
319 gls(iend+1,j,k,nout)=gls(iend+1,j,k,nout)* &
320 &
grid(ng)%rmask(iend+1,j)
321# endif
322 END IF
323 END DO
324 END DO
325
326
327
330 DO j=jstr,jend
332 tke(iend+1,j,k,nout)=tke(iend,j,k,nout)
333# ifdef MASKING
334 tke(iend+1,j,k,nout)=tke(iend+1,j,k,nout)* &
335 &
grid(ng)%rmask(iend+1,j)
336# endif
337 gls(iend+1,j,k,nout)=gls(iend,j,k,nout)
338# ifdef MASKING
339 gls(iend+1,j,k,nout)=gls(iend+1,j,k,nout)* &
340 &
grid(ng)%rmask(iend+1,j)
341# endif
342 END IF
343 END DO
344 END DO
345
346
347
350 DO j=jstr,jend
352 tke(iend+1,j,k,nout)=tke(iend,j,k,nout)
353# ifdef MASKING
354 tke(iend+1,j,k,nout)=tke(iend+1,j,k,nout)* &
355 &
grid(ng)%rmask(iend+1,j)
356# endif
357 gls(iend+1,j,k,nout)=gls(iend,j,k,nout)
358# ifdef MASKING
359 gls(iend+1,j,k,nout)=gls(iend+1,j,k,nout)* &
360 &
grid(ng)%rmask(iend+1,j)
361# endif
362 END IF
363 END DO
364 END DO
365 END IF
366 END IF
367
368
369
370
371
372 IF (
domain(ng)%Southern_Edge(tile))
THEN
373
374
375
378 DO i=istr,iend+1
379 grad(i,jstr )=tke(i ,jstr ,k,nstp)- &
380 & tke(i-1,jstr ,k,nstp)
381# ifdef MASKING
382 grad(i,jstr )=grad(i,jstr )*
grid(ng)%umask(i,jstr )
383# endif
384 grad(i,jstr-1)=tke(i ,jstr-1,k,nstp)- &
385 & tke(i-1,jstr-1,k,nstp)
386# ifdef MASKING
387 grad(i,jstr-1)=grad(i,jstr-1)*
grid(ng)%umask(i,jstr-1)
388# endif
389 gradl(i,jstr )=gls(i ,jstr ,k,nstp)- &
390 & gls(i-1,jstr ,k,nstp)
391# ifdef MASKING
392 gradl(i,jstr )=gradl(i,jstr )*
grid(ng)%umask(i,jstr )
393# endif
394 gradl(i,jstr-1)=gls(i ,jstr-1,k,nstp)- &
395 & gls(i-1,jstr-1,k,nstp)
396# ifdef MASKING
397 gradl(i,jstr-1)=gradl(i,jstr-1)*
grid(ng)%umask(i,jstr-1)
398# endif
399 END DO
400 DO i=istr,iend
402 dkdt=tke(i,jstr,k,nstp)-tke(i,jstr ,k,nout)
403 dkde=tke(i,jstr,k,nout)-tke(i,jstr+1,k,nout)
404 IF ((dkdt*dkde).lt.0.0_r8) dkdt=0.0_r8
405 IF ((dkdt*(grad(i ,jstr)+ &
406 & grad(i+1,jstr))).gt.0.0_r8) THEN
407 dkdx=grad(i ,jstr)
408 ELSE
409 dkdx=grad(i+1,jstr)
410 END IF
411 cff=max(dkdx*dkdx+dkde*dkde, eps)
412# ifdef RADIATION_2D
413 cx=min(cff,max(dkdt*dkdx,-cff))
414# else
415 cx=0.0_r8
416# endif
417 ce=dkdt*dkde
418 tke(i,jstr-1,k,nout)=(cff*tke(i,jstr-1,k,nstp)+ &
419 & ce *tke(i,jstr ,k,nout)- &
420 & max(cx,0.0_r8)* &
421 & grad(i ,jstr-1)- &
422 & min(cx,0.0_r8)* &
423 & grad(i+1,jstr-1))/ &
424 & (cff+ce)
425# ifdef MASKING
426 tke(i,jstr-1,k,nout)=tke(i,jstr-1,k,nout)* &
427 &
grid(ng)%rmask(i,jstr-1)
428# endif
429 dkdt=gls(i,jstr,k,nstp)-gls(i,jstr ,k,nout)
430 dkde=gls(i,jstr,k,nout)-gls(i,jstr+1,k,nout)
431 IF ((dkdt*dkde).lt.0.0_r8) dkdt=0.0_r8
432 IF ((dkdt*(gradl(i ,jstr)+ &
433 & gradl(i+1,jstr))).gt.0.0_r8) THEN
434 dkdx=gradl(i ,jstr)
435 ELSE
436 dkdx=gradl(i+1,jstr)
437 END IF
438 cff=max(dkdx*dkdx+dkde*dkde,eps)
439# ifdef RADIATION_2D
440 cx=min(cff,max(dkdt*dkdx,-cff))
441# else
442 cx=0.0_r8
443# endif
444 ce=dkdt*dkde
445 gls(i,jstr-1,k,nout)=(cff*gls(i,jstr-1,k,nstp)+ &
446 & ce *gls(i,jstr ,k,nout)- &
447 & max(cx,0.0_r8)* &
448 & gradl(i ,jstr-1)- &
449 & min(cx,0.0_r8)* &
450 & gradl(i+1,jstr-1))/ &
451 & (cff+ce)
452# ifdef MASKING
453 gls(i,jstr-1,k,nout)=gls(i,jstr-1,k,nout)* &
454 &
grid(ng)%rmask(i,jstr-1)
455# endif
456 END IF
457 END DO
458 END DO
459
460
461
464 DO i=istr,iend
466 tke(i,jstr-1,k,nout)=tke(i,jstr,k,nout)
467# ifdef MASKING
468 tke(i,jstr-1,k,nout)=tke(i,jstr-1,k,nout)* &
469 &
grid(ng)%rmask(i,jstr-1)
470# endif
471 gls(i,jstr-1,k,nout)=gls(i,jstr,k,nout)
472# ifdef MASKING
473 gls(i,jstr-1,k,nout)=gls(i,jstr-1,k,nout)* &
474 &
grid(ng)%rmask(i,jstr-1)
475# endif
476 END IF
477 END DO
478 END DO
479
480
481
484 DO i=istr,iend
486 tke(i,jstr-1,k,nout)=tke(i,jstr,k,nout)
487# ifdef MASKING
488 tke(i,jstr-1,k,nout)=tke(i,jstr-1,k,nout)* &
489 &
grid(ng)%rmask(i,jstr-1)
490# endif
491 gls(i,jstr-1,k,nout)=gls(i,jstr,k,nout)
492# ifdef MASKING
493 gls(i,jstr-1,k,nout)=gls(i,jstr-1,k,nout)* &
494 &
grid(ng)%rmask(i,jstr-1)
495# endif
496 END IF
497 END DO
498 END DO
499 END IF
500 END IF
501
502
503
504
505
506 IF (
domain(ng)%Northern_Edge(tile))
THEN
507
508
509
512 DO i=istr,iend+1
513 grad(i,jend )=tke(i ,jend ,k,nstp)- &
514 & tke(i-1,jend ,k,nstp)
515# ifdef MASKING
516 grad(i,jend )=grad(i,jend )* &
517 &
grid(ng)%umask(i,jend )
518# endif
519 grad(i,jend+1)=tke(i ,jend+1,k,nstp)- &
520 & tke(i-1,jend+1,k,nstp)
521# ifdef MASKING
522 grad(i,jend+1)=grad(i,jend+1)* &
523 &
grid(ng)%umask(i,jend+1)
524# endif
525 gradl(i,jend )=gls(i ,jend ,k,nstp)- &
526 & gls(i-1,jend ,k,nstp)
527# ifdef MASKING
528 gradl(i,jend )=gradl(i,jend )* &
529 &
grid(ng)%umask(i,jend )
530# endif
531 gradl(i,jend+1)=gls(i ,jend+1,k,nstp)- &
532 & gls(i-1,jend+1,k,nstp)
533# ifdef MASKING
534 gradl(i,jend+1)=gradl(i,jend+1)* &
535 &
grid(ng)%umask(i,jend+1)
536# endif
537 END DO
538 DO i=istr,iend
540 dkdt=tke(i,jend,k,nstp)-tke(i,jend ,k,nout)
541 dkde=tke(i,jend,k,nout)-tke(i,jend-1,k,nout)
542 IF ((dkdt*dkde).lt.0.0_r8) dkdt=0.0_r8
543 IF ((dkdt*(grad(i ,jend)+ &
544 & grad(i+1,jend))).gt.0.0_r8) THEN
545 dkdx=grad(i ,jend)
546 ELSE
547 dkdx=grad(i+1,jend)
548 END IF
549 cff=max(dkdx*dkdx+dkde*dkde,eps)
550# ifdef RADIATION_2D
551 cx=min(cff,max(dkdt*dkdx,-cff))
552# else
553 cx=0.0_r8
554# endif
555 ce=dkdt*dkde
556 tke(i,jend+1,k,nout)=(cff*tke(i,jend+1,k,nstp)+ &
557 & ce *tke(i,jend ,k,nout)- &
558 & max(cx,0.0_r8)* &
559 & grad(i ,jend+1)- &
560 & min(cx,0.0_r8)* &
561 & grad(i+1,jend+1))/ &
562 & (cff+ce)
563# ifdef MASKING
564 tke(i,jend+1,k,nout)=tke(i,jend+1,k,nout)* &
565 &
grid(ng)%rmask(i,jend+1)
566# endif
567 dkdt=gls(i,jend,k,nstp)-gls(i,jend ,k,nout)
568 dkde=gls(i,jend,k,nout)-gls(i,jend-1,k,nout)
569 IF ((dkdt*dkde).lt.0.0_r8) dkdt=0.0_r8
570 IF ((dkdt*(gradl(i ,jend)+ &
571 & gradl(i+1,jend))).gt.0.0_r8) THEN
572 dkdx=gradl(i ,jend)
573 ELSE
574 dkdx=gradl(i+1,jend)
575 END IF
576 cff=max(dkdx*dkdx+dkde*dkde,eps)
577# ifdef RADIATION_2D
578 cx=min(cff,max(dkdt*dkdx,-cff))
579# else
580 cx=0.0_r8
581# endif
582 ce=dkdt*dkde
583 gls(i,jend+1,k,nout)=(cff*gls(i,jend+1,k,nstp)+ &
584 & ce *gls(i,jend ,k,nout)- &
585 & max(cx,0.0_r8)* &
586 & gradl(i ,jend+1)- &
587 & min(cx,0.0_r8)* &
588 & gradl(i+1,jend+1))/ &
589 & (cff+ce)
590# ifdef MASKING
591 gls(i,jend+1,k,nout)=gls(i,jend+1,k,nout)* &
592 &
grid(ng)%rmask(i,jend+1)
593# endif
594 END IF
595 END DO
596 END DO
597
598
599
602 DO i=istr,iend
604 tke(i,jend+1,k,nout)=tke(i,jend,k,nout)
605# ifdef MASKING
606 tke(i,jend+1,k,nout)=tke(i,jend+1,k,nout)* &
607 &
grid(ng)%rmask(i,jend+1)
608# endif
609 gls(i,jend+1,k,nout)=gls(i,jend,k,nout)
610# ifdef MASKING
611 gls(i,jend+1,k,nout)=gls(i,jend+1,k,nout)* &
612 &
grid(ng)%rmask(i,jend+1)
613# endif
614 END IF
615 END DO
616 END DO
617
618
619
622 DO i=istr,iend
624 tke(i,jend+1,k,nout)=tke(i,jend,k,nout)
625# ifdef MASKING
626 tke(i,jend+1,k,nout)=tke(i,jend+1,k,nout)* &
627 &
grid(ng)%rmask(i,jend+1)
628# endif
629 gls(i,jend+1,k,nout)=gls(i,jend,k,nout)
630# ifdef MASKING
631 gls(i,jend+1,k,nout)=gls(i,jend+1,k,nout)* &
632 &
grid(ng)%rmask(i,jend+1)
633# endif
634 END IF
635 END DO
636 END DO
637 END IF
638 END IF
639
640
641
642
643
645 IF (
domain(ng)%SouthWest_Corner(tile))
THEN
649 tke(istr-1,jstr-1,k,nout)=0.5_r8* &
650 & (tke(istr ,jstr-1,k,nout)+ &
651 & tke(istr-1,jstr ,k,nout))
652 gls(istr-1,jstr-1,k,nout)=0.5_r8* &
653 & (gls(istr ,jstr-1,k,nout)+ &
654 & gls(istr-1,jstr ,k,nout))
655 END DO
656 END IF
657 END IF
658 IF (
domain(ng)%SouthEast_Corner(tile))
THEN
662 tke(iend+1,jstr-1,k,nout)=0.5_r8* &
663 & (tke(iend ,jstr-1,k,nout)+ &
664 & tke(iend+1,jstr ,k,nout))
665 gls(iend+1,jstr-1,k,nout)=0.5_r8* &
666 & (gls(iend ,jstr-1,k,nout)+ &
667 & gls(iend+1,jstr ,k,nout))
668 END DO
669 END IF
670 END IF
671 IF (
domain(ng)%NorthWest_Corner(tile))
THEN
675 tke(istr-1,jend+1,k,nout)=0.5_r8* &
676 & (tke(istr ,jend+1,k,nout)+ &
677 & tke(istr-1,jend ,k,nout))
678 gls(istr-1,jend+1,k,nout)=0.5_r8* &
679 & (gls(istr ,jend+1,k,nout)+ &
680 & gls(istr-1,jend ,k,nout))
681 END DO
682 END IF
683 END IF
684 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
688 tke(iend+1,jend+1,k,nout)=0.5_r8* &
689 & (tke(iend ,jend+1,k,nout)+ &
690 & tke(iend+1,jend ,k,nout))
691 gls(iend+1,jend+1,k,nout)=0.5_r8* &
692 & (gls(iend ,jend+1,k,nout)+ &
693 & gls(iend+1,jend ,k,nout))
694 END DO
695 END IF
696 END IF
697 END IF
698
699 RETURN
type(t_apply), dimension(:), allocatable lbc_apply
type(t_grid), dimension(:), allocatable grid
type(t_lbc), dimension(:,:,:), allocatable lbc
type(t_domain), dimension(:), allocatable domain
logical, dimension(:), allocatable ewperiodic
logical, dimension(:), allocatable nsperiodic
integer, parameter isouth
integer, parameter inorth