3#if defined SOLVE3D && (defined MY25_MIXING || defined GLS_MIXING)
26 SUBROUTINE tkebc (ng, tile, nout)
35 integer,
intent(in) :: ng, tile, nout
42 & lbi, ubi, lbj, ubj,
n(ng), &
43 & imins, imaxs, jmins, jmaxs, &
52 & LBi, UBi, LBj, UBj, UBk, &
53 & IminS, ImaxS, JminS, JmaxS, &
65 integer,
intent(in) :: ng, tile
66 integer,
intent(in) :: LBi, UBi, LBj, UBj, UBk
67 integer,
intent(in) :: IminS, ImaxS, JminS, JmaxS
68 integer,
intent(in) :: nout, nstp
71 real(r8),
intent(inout) :: gls(LBi:,LBj:,0:,:)
72 real(r8),
intent(inout) :: tke(LBi:,LBj:,0:,:)
74 real(r8),
intent(inout) :: gls(LBi:UBi,LBj:UBj,0:UBk,3)
75 real(r8),
intent(inout) :: tke(LBi:UBi,LBj:UBj,0:UBk,3)
82 real(r8),
parameter :: eps = 1.0e-20_r8
84 real(r8) :: Ce, Cx, cff, dKde, dKdt, dKdx
86 real(r8),
dimension(IminS:ImaxS,JminS:JmaxS) :: grad
87 real(r8),
dimension(IminS:ImaxS,JminS:JmaxS) :: gradL
89# include "set_bounds.h"
95 IF (
domain(ng)%Western_Edge(tile))
THEN
102 grad(istr-1,j)=tke(istr-1,j ,k,nstp)- &
103 & tke(istr-1,j-1,k,nstp)
105 grad(istr-1,j)=grad(istr-1,j)* &
106 &
grid(ng)%vmask(istr-1,j)
108 grad(istr ,j)=tke(istr ,j ,k,nstp)- &
109 & tke(istr ,j-1,k,nstp)
111 grad(istr ,j)=grad(istr ,j)* &
112 &
grid(ng)%vmask(istr ,j)
114 gradl(istr-1,j)=gls(istr-1,j ,k,nstp)- &
115 & gls(istr-1,j-1,k,nstp)
117 gradl(istr-1,j)=gradl(istr-1,j)* &
118 &
grid(ng)%vmask(istr-1,j)
120 gradl(istr ,j)=gls(istr ,j ,k,nstp)- &
121 & gls(istr ,j-1,k,nstp)
123 gradl(istr ,j)=gradl(istr ,j)* &
124 &
grid(ng)%vmask(istr ,j)
128 IF (lbc_apply(ng)%west(j))
THEN
129 dkdt=tke(istr,j,k,nstp)-tke(istr ,j,k,nout)
130 dkdx=tke(istr,j,k,nstp)-tke(istr+1,j,k,nstp)
131 IF ((dkdt*dkdx).lt.0.0_r8) dkdt=0.0_r8
132 IF ((dkdt*(grad(istr,j )+ &
133 & grad(istr,j+1))).gt.0.0_r8)
THEN
138 cff=dkdt/max(dkdx*dkdx+dkde*dkde,eps)
139 cx=min(1.0_r8,cff*dkdx)
141 ce=min(1.0_r8,max(cff*dkde,-1.0_r8))
145 tke(istr-1,j,k,nout)=(1.0_r8-cx)*tke(istr-1,j,k,nstp)+ &
146 & cx*tke(istr,j,k,nstp)- &
147 & max(ce,0.0_r8)*grad(istr-1,j )- &
148 & min(ce,0.0_r8)*grad(istr-1,j+1)
150 tke(istr-1,j,k,nout)=tke(istr-1,j,k,nout)* &
151 &
grid(ng)%rmask(istr-1,j)
153 dkdt=gls(istr,j,k,nstp)-gls(istr ,j,k,nout)
154 dkdx=gls(istr,j,k,nstp)-gls(istr+1,j,k,nstp)
155 IF ((dkdt*dkdx).lt.0.0_r8) dkdt=0.0_r8
156 IF ((dkdt*(gradl(istr,j )+ &
157 & gradl(istr,j+1))).gt.0.0_r8)
THEN
162 cff=dkdt/max(dkdx*dkdx+dkde*dkde,eps)
163 cx=min(1.0_r8,cff*dkdx)
165 ce=min(1.0_r8,max(cff*dkde,-1.0_r8))
169 gls(istr-1,j,k,nout)=(1.0_r8-cx)*gls(istr-1,j,k,nstp)+ &
170 & cx*gls(istr,j,k,nstp)- &
171 & max(ce,0.0_r8)*gradl(istr-1,j )- &
172 & min(ce,0.0_r8)*gradl(istr-1,j+1)
174 gls(istr-1,j,k,nout)=gls(istr-1,j,k,nout)* &
175 &
grid(ng)%rmask(istr-1,j)
186 IF (lbc_apply(ng)%west(j))
THEN
187 tke(istr-1,j,k,nout)=tke(istr,j,k,nout)
189 tke(istr-1,j,k,nout)=tke(istr-1,j,k,nout)* &
190 &
grid(ng)%rmask(istr-1,j)
192 gls(istr-1,j,k,nout)=gls(istr,j,k,nout)
194 gls(istr-1,j,k,nout)=gls(istr-1,j,k,nout)* &
195 &
grid(ng)%rmask(istr-1,j)
206 IF (lbc_apply(ng)%west(j))
THEN
207 tke(istr-1,j,k,nout)=tke(istr,j,k,nout)
209 tke(istr-1,j,k,nout)=tke(istr-1,j,k,nout)* &
210 &
grid(ng)%rmask(istr-1,j)
212 gls(istr-1,j,k,nout)=gls(istr,j,k,nout)
214 gls(istr-1,j,k,nout)=gls(istr-1,j,k,nout)* &
215 &
grid(ng)%rmask(istr-1,j)
227 IF (
domain(ng)%Eastern_Edge(tile))
THEN
234 grad(iend ,j)=tke(iend ,j ,k,nstp)- &
235 & tke(iend ,j-1,k,nstp)
237 grad(iend ,j)=grad(iend ,j)* &
238 &
grid(ng)%vmask(iend ,j)
240 grad(iend+1,j)=tke(iend+1,j ,k,nstp)- &
241 & tke(iend+1,j-1,k,nstp)
243 grad(iend+1,j)=grad(iend+1,j)* &
244 &
grid(ng)%vmask(iend+1,j)
246 gradl(iend ,j)=gls(iend ,j ,k,nstp)- &
247 & gls(iend ,j-1,k,nstp)
249 gradl(iend ,j)=gradl(iend ,j)* &
250 &
grid(ng)%vmask(iend ,j)
252 gradl(iend+1,j)=gls(iend+1,j ,k,nstp)- &
253 & gls(iend+1,j-1,k,nstp)
255 gradl(iend+1,j)=gradl(iend+1,j)* &
256 &
grid(ng)%vmask(iend+1,j)
260 IF (lbc_apply(ng)%east(j))
THEN
261 dkdt=tke(iend,j,k,nstp)-tke(iend ,j,k,nout)
262 dkdx=tke(iend,j,k,nstp)-tke(iend-1,j,k,nstp)
263 IF ((dkdt*dkdx).lt.0.0_r8) dkdt=0.0_r8
264 IF ((dkdt*(grad(iend,j )+ &
265 & grad(iend,j+1))).gt.0.0_r8)
THEN
270 cff=dkdt/max(dkdx*dkdx+dkde*dkde,eps)
271 cx=min(1.0_r8,cff*dkdx)
273 ce=min(1.0_r8,max(cff*dkde,-1.0_r8))
277 tke(iend+1,j,k,nout)=(1.0_r8-cx)*tke(iend+1,j,k,nstp)+ &
278 & cx*tke(iend,j,k,nstp)- &
279 & max(ce,0.0_r8)*grad(iend+1,j )- &
280 & min(ce,0.0_r8)*grad(iend+1,j+1)
282 tke(iend+1,j,k,nout)=tke(iend+1,j,k,nout)* &
283 &
grid(ng)%rmask(iend+1,j)
285 dkdt=gls(iend,j,k,nstp)-gls(iend ,j,k,nout)
286 dkdx=gls(iend,j,k,nstp)-gls(iend-1,j,k,nstp)
287 IF ((dkdt*dkdx).lt.0.0_r8) dkdt=0.0_r8
288 IF ((dkdt*(gradl(iend,j )+ &
289 & gradl(iend,j+1))).gt.0.0_r8)
THEN
294 cff=dkdt/max(dkdx*dkdx+dkde*dkde,eps)
295 cx=min(1.0_r8,cff*dkdx)
297 ce=min(1.0_r8,max(cff*dkde,-1.0_r8))
301 gls(iend+1,j,k,nout)=(1.0_r8-cx)*gls(iend+1,j,k,nstp)+ &
302 & cx*gls(iend,j,k,nstp)- &
303 & max(ce,0.0_r8)*gradl(iend+1,j )- &
304 & min(ce,0.0_r8)*gradl(iend+1,j+1)
306 gls(iend+1,j,k,nout)=gls(iend+1,j,k,nout)* &
307 &
grid(ng)%rmask(iend+1,j)
318 IF (lbc_apply(ng)%east(j))
THEN
319 tke(iend+1,j,k,nout)=tke(iend,j,k,nout)
321 tke(iend+1,j,k,nout)=tke(iend+1,j,k,nout)* &
322 &
grid(ng)%rmask(iend+1,j)
324 gls(iend+1,j,k,nout)=gls(iend,j,k,nout)
326 gls(iend+1,j,k,nout)=gls(iend+1,j,k,nout)* &
327 &
grid(ng)%rmask(iend+1,j)
338 IF (lbc_apply(ng)%east(j))
THEN
339 tke(iend+1,j,k,nout)=tke(iend,j,k,nout)
341 tke(iend+1,j,k,nout)=tke(iend+1,j,k,nout)* &
342 &
grid(ng)%rmask(iend+1,j)
344 gls(iend+1,j,k,nout)=gls(iend,j,k,nout)
346 gls(iend+1,j,k,nout)=gls(iend+1,j,k,nout)* &
347 &
grid(ng)%rmask(iend+1,j)
359 IF (
domain(ng)%Southern_Edge(tile))
THEN
366 grad(i,jstr )=tke(i ,jstr ,k,nstp)- &
367 & tke(i-1,jstr ,k,nstp)
369 grad(i,jstr )=grad(i,jstr )*
grid(ng)%umask(i,jstr )
371 grad(i,jstr-1)=tke(i ,jstr-1,k,nstp)- &
372 & tke(i-1,jstr-1,k,nstp)
374 grad(i,jstr-1)=grad(i,jstr-1)*
grid(ng)%umask(i,jstr-1)
376 gradl(i,jstr )=gls(i ,jstr ,k,nstp)- &
377 & gls(i-1,jstr ,k,nstp)
379 gradl(i,jstr )=gradl(i,jstr )*
grid(ng)%umask(i,jstr )
381 gradl(i,jstr-1)=gls(i ,jstr-1,k,nstp)- &
382 & gls(i-1,jstr-1,k,nstp)
384 gradl(i,jstr-1)=gradl(i,jstr-1)*
grid(ng)%umask(i,jstr-1)
388 IF (lbc_apply(ng)%south(i))
THEN
389 dkdt=tke(i,jstr,k,nstp)-tke(i,jstr ,k,nout)
390 dkde=tke(i,jstr,k,nstp)-tke(i,jstr+1,k,nstp)
391 IF ((dkdt*dkde).lt.0.0_r8) dkdt=0.0_r8
392 IF ((dkdt*(grad(i ,jstr)+ &
393 & grad(i+1,jstr))).gt.0.0_r8)
THEN
398 cff=dkdt/max(dkdx*dkdx+dkde*dkde,eps)
400 cx=min(1.0_r8,max(cff*dkdx,-1.0_r8))
404 ce=min(1.0_r8,cff*dkde)
405 tke(i,jstr-1,k,nout)=(1.0_r8-ce)*tke(i,jstr-1,k,nstp)+ &
406 & ce*tke(i,jstr,k,nstp)- &
407 & max(cx,0.0_r8)*grad(i ,jstr-1)- &
408 & min(cx,0.0_r8)*grad(i+1,jstr-1)
410 tke(i,jstr-1,k,nout)=tke(i,jstr-1,k,nout)* &
411 &
grid(ng)%rmask(i,jstr-1)
413 dkdt=gls(i,jstr,k,nstp)-gls(i,jstr ,k,nout)
414 dkde=gls(i,jstr,k,nstp)-gls(i,jstr+1,k,nstp)
415 IF ((dkdt*dkde).lt.0.0_r8) dkdt=0.0_r8
416 IF ((dkdt*(gradl(i ,jstr)+ &
417 & gradl(i+1,jstr))).gt.0.0_r8)
THEN
422 cff=dkdt/max(dkdx*dkdx+dkde*dkde,eps)
424 cx=min(1.0_r8,max(cff*dkdx,-1.0_r8))
428 ce=min(1.0_r8,cff*dkde)
429 gls(i,jstr-1,k,nout)=(1.0_r8-ce)*gls(i,jstr-1,k,nstp)+ &
430 & ce*gls(i,jstr,k,nstp)- &
431 & max(cx,0.0_r8)*gradl(i ,jstr-1)- &
432 & min(cx,0.0_r8)*gradl(i+1,jstr-1)
434 gls(i,jstr-1,k,nout)=gls(i,jstr-1,k,nout)* &
435 &
grid(ng)%rmask(i,jstr-1)
446 IF (lbc_apply(ng)%south(i))
THEN
447 tke(i,jstr-1,k,nout)=tke(i,jstr,k,nout)
449 tke(i,jstr-1,k,nout)=tke(i,jstr-1,k,nout)* &
450 &
grid(ng)%rmask(i,jstr-1)
452 gls(i,jstr-1,k,nout)=gls(i,jstr,k,nout)
454 gls(i,jstr-1,k,nout)=gls(i,jstr-1,k,nout)* &
455 &
grid(ng)%rmask(i,jstr-1)
466 IF (lbc_apply(ng)%south(i))
THEN
467 tke(i,jstr-1,k,nout)=tke(i,jstr,k,nout)
469 tke(i,jstr-1,k,nout)=tke(i,jstr-1,k,nout)* &
470 &
grid(ng)%rmask(i,jstr-1)
472 gls(i,jstr-1,k,nout)=gls(i,jstr,k,nout)
474 gls(i,jstr-1,k,nout)=gls(i,jstr-1,k,nout)* &
475 &
grid(ng)%rmask(i,jstr-1)
487 IF (
domain(ng)%Northern_Edge(tile))
THEN
494 grad(i,jend )=tke(i ,jend ,k,nstp)- &
495 & tke(i-1,jend ,k,nstp)
497 grad(i,jend )=grad(i,jend )* &
498 &
grid(ng)%umask(i,jend )
500 grad(i,jend+1)=tke(i ,jend+1,k,nstp)- &
501 & tke(i-1,jend+1,k,nstp)
503 grad(i,jend+1)=grad(i,jend+1)* &
504 &
grid(ng)%umask(i,jend+1)
506 gradl(i,jend )=gls(i ,jend ,k,nstp)- &
507 & gls(i-1,jend ,k,nstp)
509 gradl(i,jend )=gradl(i,jend )* &
510 &
grid(ng)%umask(i,jend )
512 gradl(i,jend+1)=gls(i ,jend+1,k,nstp)- &
513 & gls(i-1,jend+1,k,nstp)
515 gradl(i,jend+1)=gradl(i,jend+1)* &
516 &
grid(ng)%umask(i,jend+1)
520 IF (lbc_apply(ng)%north(i))
THEN
521 dkdt=tke(i,jend,k,nstp)-tke(i,jend ,k,nout)
522 dkde=tke(i,jend,k,nstp)-tke(i,jend-1,k,nstp)
523 IF ((dkdt*dkde).lt.0.0_r8) dkdt=0.0_r8
524 IF ((dkdt*(grad(i ,jend)+ &
525 & grad(i+1,jend))).gt.0.0_r8)
THEN
530 cff=dkdt/max(dkdx*dkdx+dkde*dkde,eps)
532 cx=min(1.0_r8,max(cff*dkdx,-1.0_r8))
536 ce=min(1.0_r8,cff*dkde)
537 tke(i,jend+1,k,nout)=(1.0_r8-ce)*tke(i,jend+1,k,nstp)+ &
538 & ce*tke(i,jend,k,nstp)- &
539 & max(cx,0.0_r8)*grad(i ,jend+1)- &
540 & min(cx,0.0_r8)*grad(i+1,jend+1)
542 tke(i,jend+1,k,nout)=tke(i,jend+1,k,nout)* &
543 &
grid(ng)%rmask(i,jend+1)
545 dkdt=gls(i,jend,k,nstp)-gls(i,jend ,k,nout)
546 dkde=gls(i,jend,k,nstp)-gls(i,jend-1,k,nstp)
547 IF ((dkdt*dkde).lt.0.0_r8) dkdt=0.0_r8
548 IF ((dkdt*(gradl(i ,jend)+ &
549 & gradl(i+1,jend))).gt.0.0_r8)
THEN
554 cff=dkdt/max(dkdx*dkdx+dkde*dkde,eps)
556 cx=min(1.0_r8,max(cff*dkdx,-1.0_r8))
560 ce=min(1.0_r8,cff*dkde)
561 gls(i,jend+1,k,nout)=(1.0_r8-ce)*gls(i,jend+1,k,nstp)+ &
562 & ce*gls(i,jend,k,nstp)- &
563 & max(cx,0.0_r8)*gradl(i ,jend+1)- &
564 & min(cx,0.0_r8)*gradl(i+1,jend+1)
566 gls(i,jend+1,k,nout)=gls(i,jend+1,k,nout)* &
567 &
grid(ng)%rmask(i,jend+1)
578 IF (lbc_apply(ng)%north(i))
THEN
579 tke(i,jend+1,k,nout)=tke(i,jend,k,nout)
581 tke(i,jend+1,k,nout)=tke(i,jend+1,k,nout)* &
582 &
grid(ng)%rmask(i,jend+1)
584 gls(i,jend+1,k,nout)=gls(i,jend,k,nout)
586 gls(i,jend+1,k,nout)=gls(i,jend+1,k,nout)* &
587 &
grid(ng)%rmask(i,jend+1)
598 IF (lbc_apply(ng)%north(i))
THEN
599 tke(i,jend+1,k,nout)=tke(i,jend,k,nout)
601 tke(i,jend+1,k,nout)=tke(i,jend+1,k,nout)* &
602 &
grid(ng)%rmask(i,jend+1)
604 gls(i,jend+1,k,nout)=gls(i,jend,k,nout)
606 gls(i,jend+1,k,nout)=gls(i,jend+1,k,nout)* &
607 &
grid(ng)%rmask(i,jend+1)
620 IF (
domain(ng)%SouthWest_Corner(tile))
THEN
621 IF (lbc_apply(ng)%south(istr-1).and. &
622 & lbc_apply(ng)%west (jstr-1))
THEN
624 tke(istr-1,jstr-1,k,nout)=0.5_r8* &
625 & (tke(istr ,jstr-1,k,nout)+ &
626 & tke(istr-1,jstr ,k,nout))
627 gls(istr-1,jstr-1,k,nout)=0.5_r8* &
628 & (gls(istr ,jstr-1,k,nout)+ &
629 & gls(istr-1,jstr ,k,nout))
633 IF (
domain(ng)%SouthEast_Corner(tile))
THEN
634 IF (lbc_apply(ng)%south(iend+1).and. &
635 & lbc_apply(ng)%east (jstr-1))
THEN
637 tke(iend+1,jstr-1,k,nout)=0.5_r8* &
638 & (tke(iend ,jstr-1,k,nout)+ &
639 & tke(iend+1,jstr ,k,nout))
640 gls(iend+1,jstr-1,k,nout)=0.5_r8* &
641 & (gls(iend ,jstr-1,k,nout)+ &
642 & gls(iend+1,jstr ,k,nout))
646 IF (
domain(ng)%NorthWest_Corner(tile))
THEN
647 IF (lbc_apply(ng)%north(istr-1).and. &
648 & lbc_apply(ng)%west (jend+1))
THEN
650 tke(istr-1,jend+1,k,nout)=0.5_r8* &
651 & (tke(istr ,jend+1,k,nout)+ &
652 & tke(istr-1,jend ,k,nout))
653 gls(istr-1,jend+1,k,nout)=0.5_r8* &
654 & (gls(istr ,jend+1,k,nout)+ &
655 & gls(istr-1,jend ,k,nout))
659 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
660 IF (lbc_apply(ng)%north(iend+1).and. &
661 & lbc_apply(ng)%east (jend+1))
THEN
663 tke(iend+1,jend+1,k,nout)=0.5_r8* &
664 & (tke(iend ,jend+1,k,nout)+ &
665 & tke(iend+1,jend ,k,nout))
666 gls(iend+1,jend+1,k,nout)=0.5_r8* &
667 & (gls(iend ,jend+1,k,nout)+ &
668 & gls(iend+1,jend ,k,nout))
type(t_grid), dimension(:), allocatable grid
type(t_mixing), dimension(:), allocatable mixing
integer, dimension(:), allocatable n
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
integer, dimension(:), allocatable nstp
subroutine tkebc(ng, tile, nout)
subroutine, public tkebc_tile(ng, tile, lbi, ubi, lbj, ubj, ubk, imins, imaxs, jmins, jmaxs, nout, nstp, gls, tke)