132
133
137
139# ifdef DISTRIBUTE
141# endif
143
144
145
146 integer, intent(in) :: ng, tile
147 integer, intent(in) :: LBi, UBi, LBj, UBj
148 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
149 integer, intent(in) :: nstp, nnew
150
151# ifdef ASSUMED_SHAPE
152# ifdef MASKING
153 real(r8), intent(in) :: umask(LBi:,LBj:)
154 real(r8), intent(in) :: vmask(LBi:,LBj:)
155# endif
156 real(r8), intent(in) :: Huon(LBi:,LBj:,:)
157 real(r8), intent(in) :: Hvom(LBi:,LBj:,:)
158 real(r8), intent(in) :: Hz(LBi:,LBj:,:)
159 real(r8), intent(in) :: pm(LBi:,LBj:)
160 real(r8), intent(in) :: pn(LBi:,LBj:)
161 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
162 real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
163 real(r8), intent(in) :: ZoBot(LBi:,LBj:)
164 real(r8), intent(in) :: u(LBi:,LBj:,:,:)
165 real(r8), intent(in) :: v(LBi:,LBj:,:,:)
166 real(r8), intent(in) :: W(LBi:,LBj:,0:)
167# ifdef WEC_VF
168 real(r8), intent(in) :: W_stokes(LBi:,LBj:,0:)
169# endif
170 real(r8), intent(in) :: bustr(LBi:,LBj:)
171 real(r8), intent(in) :: bvstr(LBi:,LBj:)
172 real(r8), intent(in) :: sustr(LBi:,LBj:)
173 real(r8), intent(in) :: svstr(LBi:,LBj:)
174# ifdef ZOS_HSIG
175 real(r8), intent(in) :: Hwave(LBi:,LBj:)
176# endif
177# ifdef TKE_WAVEDISS
178 real(r8), intent(in) :: Dissip_break(LBi:,LBj:)
179 real(r8), intent(in) :: Dissip_wcap(LBi:,LBj:)
180# endif
181 real(r8), intent(in) :: bvf(LBi:,LBj:,0:)
182
183 real(r8), intent(inout) :: Akt(LBi:,LBj:,0:,:)
184 real(r8), intent(inout) :: Akv(LBi:,LBj:,0:)
185 real(r8), intent(inout) :: Akk(LBi:,LBj:,0:)
186 real(r8), intent(inout) :: Akp(LBi:,LBj:,0:)
187 real(r8), intent(inout) :: Lscale(LBi:,LBj:,0:)
188 real(r8), intent(inout) :: gls(LBi:,LBj:,0:,:)
189 real(r8), intent(inout) :: tke(LBi:,LBj:,0:,:)
190# else
191# ifdef MASKING
192 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
193 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
194# endif
195 real(r8), intent(in) :: Huon(LBi:UBi,LBj:UBj,N(ng))
196 real(r8), intent(in) :: Hvom(LBi:UBi,LBj:UBj,N(ng))
197 real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
198 real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
199 real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
200 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
201 real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng))
202 real(r8), intent(in) :: ZoBot(LBi:UBi,LBj:UBj)
203 real(r8), intent(in) :: u(LBi:UBi,LBj:UBj,N(ng),2)
204 real(r8), intent(in) :: v(LBi:UBi,LBj:UBj,N(ng),2)
205 real(r8), intent(in) :: W(LBi:UBi,LBj:UBj,0:N(ng))
206# ifdef WEC_VF
207 real(r8), intent(in) :: W_stokes(LBi:UBi,LBj:UBj,0:N(ng))
208# endif
209 real(r8), intent(in) :: bustr(LBi:UBi,LBj:UBj)
210 real(r8), intent(in) :: bvstr(LBi:UBi,LBj:UBj)
211 real(r8), intent(in) :: sustr(LBi:UBi,LBj:UBj)
212 real(r8), intent(in) :: svstr(LBi:UBi,LBj:UBj)
213# if defined ZOS_HSIG
214 real(r8), intent(in) :: Hwave(LBi:UBi,LBj:UBj)
215# endif
216# ifdef TKE_WAVEDISS
217 real(r8), intent(in) :: Dissip_break(LBi:UBi,LBj:UBj)
218 real(r8), intent(in) :: Dissip_wcap(LBi:UBi,LBj:UBj)
219# endif
220 real(r8), intent(in) :: bvf(LBi:UBi,LBj:UBj,0:N(ng))
221
222 real(r8), intent(inout) :: Akt(LBi:UBi,LBj:UBj,0:N(ng),NAT)
223 real(r8), intent(inout) :: Akv(LBi:UBi,LBj:UBj,0:N(ng))
224 real(r8), intent(inout) :: Akk(LBi:UBi,LBj:UBj,0:N(ng))
225 real(r8), intent(inout) :: Akp(LBi:UBi,LBj:UBj,0:N(ng))
226 real(r8), intent(inout) :: Lscale(LBi:UBi,LBj:UBj,0:N(ng))
227 real(r8), intent(inout) :: gls(LBi:UBi,LBj:UBj,0:N(ng),3)
228 real(r8), intent(inout) :: tke(LBi:UBi,LBj:UBj,0:N(ng),3)
229# endif
230
231
232
233 logical :: Lmy25
234 integer :: i, itrc, j, k
235
236 real(r8), parameter :: Gadv = 1.0_r8/3.0_r8
237 real(r8), parameter :: eps = 1.0e-10_r8
238 real(r8) :: Zos_min
239
240 real(r8) :: Gh, Gm, Kprod, Ls_unlmt, Ls_lmt, Pprod, Sh, Sm
241 real(r8) :: cff, cff1, cff2, cff3
242 real(r8) :: cmu_fac1, cmu_fac2, cmu_fac3, cmu_fac4
243 real(r8) :: gls_c3, gls_exp1, gls_fac1, gls_fac2, gls_fac3
244 real(r8) :: gls_fac4, gls_fac5, gls_fac6, ql, sqrt2, strat2
245 real(r8) :: tke_exp1, tke_exp2, tke_exp3, tke_exp4, wall_fac
246 real(r8) :: gls_d, gls_sigp_cb, ogls_sigp, sig_eff
247 real(r8) :: L_sft
248# if defined CRAIG_BANNER || defined TKE_WAVEDISS
249 real(r8) :: cb_wallE
250# endif
251
252 real(r8), dimension(IminS:ImaxS) :: tke_fluxt
253 real(r8), dimension(IminS:ImaxS) :: tke_fluxb
254 real(r8), dimension(IminS:ImaxS) :: gls_fluxt
255 real(r8), dimension(IminS:ImaxS) :: gls_fluxb
256 real(r8), dimension(IminS:ImaxS) :: Zos_eff
257
258 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: BCK
259 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: BCP
260 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: CF
261 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FCK
262 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: FCP
263 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: dU
264 real(r8), dimension(IminS:ImaxS,0:N(ng)) :: dV
265
266 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,0:N(ng)) :: shear2
267 real(r8), dimension(IminS:ImaxS,JminS:JmaxS,0:N(ng)) :: buoy2
268
269 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FEK
270 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FEP
271 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FXK
272 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FXP
273 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Zob_min
274 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: curvK
275 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: curvP
276 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: gradK
277 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: gradP
278
279# include "set_bounds.h"
280
281
282
283
284
285 zos_min=max(
zos(ng),0.0001_r8)
286 DO j=jstr,jend
287 DO i=istr,iend
288 zob_min(i,j)=max(zobot(i,j),0.0001_r8)
289 END DO
290 END DO
291
292 lmy25=.false.
293 IF ((
gls_p(ng).eq.0.0_r8).and. &
294 & (
gls_n(ng).eq.1.0_r8).and. &
295 & (
gls_m(ng).eq.1.0_r8))
THEN
296 lmy25=.true.
297 END IF
298
299# if defined CRAIG_BANNER || defined TKE_WAVEDISS
300 IF (lmy25) THEN
301
302 cb_walle=1.25_r8
303 ELSE
304 cb_walle=1.0_r8
305 END IF
310 & cff1*
gls_n(ng)/3.0_r8*(4.0_r8*
gls_m(ng)+1.0_r8)+ &
311 & cff1**2*
gls_m(ng)/9.0_r8*(2.0_r8+4.0_r8*
gls_m(ng)))
312# else
315# endif
316 ogls_sigp=1.0_r8/gls_sigp_cb
317
318 sqrt2=sqrt(2.0_r8)
321 cmu_fac3=1.0_r8/
gls_cmu0(ng)**2.0_r8
322 cmu_fac4=(1.5_r8*
gls_sigk(ng))**(1.0_r8/3.0_r8)/ &
330 gls_fac6=8.0_r8/
gls_cmu0(ng)**6.0_r8
331
332 gls_exp1=1.0_r8/
gls_n(ng)
335 tke_exp3=0.5_r8+
gls_m(ng)
337
338
339
340
341
342# ifdef RI_SPLINES
343 DO j=jstrm1,jendp1
344 DO i=istrm1,iendp1
345 cf(i,0)=0.0_r8
346 du(i,0)=0.0_r8
347 dv(i,0)=0.0_r8
348 END DO
350 DO i=istrm1,iendp1
351 cff=1.0_r8/(2.0_r8*hz(i,j,k+1)+ &
352 & hz(i,j,k)*(2.0_r8-cf(i,k-1)))
353 cf(i,k)=cff*hz(i,j,k+1)
354 du(i,k)=cff*(3.0_r8*(u(i ,j,k+1,nstp)-u(i, j,k,nstp)+ &
355 & u(i+1,j,k+1,nstp)-u(i+1,j,k,nstp))- &
356 & hz(i,j,k)*du(i,k-1))
357 dv(i,k)=cff*(3.0_r8*(v(i,j ,k+1,nstp)-v(i,j ,k,nstp)+ &
358 & v(i,j+1,k+1,nstp)-v(i,j+1,k,nstp))- &
359 & hz(i,j,k)*dv(i,k-1))
360 END DO
361 END DO
362 DO i=istrm1,iendp1
365 END DO
367 DO i=istrm1,iendp1
368 du(i,k)=du(i,k)-cf(i,k)*du(i,k+1)
369 dv(i,k)=dv(i,k)-cf(i,k)*dv(i,k+1)
370 END DO
371 END DO
373 DO i=istrm1,iendp1
374 shear2(i,j,k)=du(i,k)*du(i,k)+dv(i,k)*dv(i,k)
375 END DO
376 END DO
377 END DO
378# else
380 DO j=jstrm1,jendp1
381 DO i=istrm1,iendp1
382 cff=0.5_r8/(z_r(i,j,k+1)-z_r(i,j,k))
383 shear2(i,j,k)=(cff*(u(i ,j,k+1,nstp)-u(i ,j,k,nstp)+ &
384 & u(i+1,j,k+1,nstp)-u(i+1,j,k,nstp)))**2+ &
385 & (cff*(v(i,j ,k+1,nstp)-v(i,j ,k,nstp)+ &
386 & v(i,j+1,k+1,nstp)-v(i,j+1,k,nstp)))**2
387 END DO
388 END DO
389 END DO
390# endif
391
392
393
395 DO j=jstr-1,jend+1
396 DO i=istr-1,iend+1
397 buoy2(i,j,k)=bvf(i,j,k)
398 END DO
399 END DO
400 END DO
401# ifdef N2S2_HORAVG
402
403
404
405
406
407
409 IF (
domain(ng)%Western_Edge(tile))
THEN
410 DO j=max(1,jstr-1),min(jend+1,
mm(ng))
411 shear2(istr-1,j,k)=shear2(istr,j,k)
412 END DO
413 END IF
414 IF (
domain(ng)%Eastern_Edge(tile))
THEN
415 DO j=max(1,jstr-1),min(jend+1,
mm(ng))
416 shear2(iend+1,j,k)=shear2(iend,j,k)
417 END DO
418 END IF
419 IF (
domain(ng)%Southern_Edge(tile))
THEN
420 DO i=max(1,istr-1),min(iend+1,
lm(ng))
421 shear2(i,jstr-1,k)=shear2(i,jstr,k)
422 END DO
423 END IF
424 IF (
domain(ng)%Northern_Edge(tile))
THEN
425 DO i=max(1,istr-1),min(iend+1,
lm(ng))
426 shear2(i,jend+1,k)=shear2(i,jend,k)
427 END DO
428 END IF
429 IF (
domain(ng)%SouthWest_Corner(tile))
THEN
430 shear2(istr-1,jstr-1,k)=shear2(istr,jstr,k)
431 END IF
432 IF (
domain(ng)%NorthWest_Corner(tile))
THEN
433 shear2(istr-1,jend+1,k)=shear2(istr,jend,k)
434 END IF
435 IF (
domain(ng)%SouthEast_Corner(tile))
THEN
436 shear2(iend+1,jstr-1,k)=shear2(iend,jstr,k)
437 END IF
438 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
439 shear2(iend+1,jend+1,k)=shear2(iend,jend,k)
440 END IF
441
442
443
444 DO j=jstr-1,jend
445 DO i=istr-1,iend
446 buoy2(i,j,0)=0.25_r8*(buoy2(i,j ,k)+buoy2(i+1,j ,k)+ &
447 & buoy2(i,j+1,k)+buoy2(i+1,j+1,k))
448 shear2(i,j,0)=0.25_r8*(shear2(i,j ,k)+shear2(i+1,j ,k)+ &
449 & shear2(i,j+1,k)+shear2(i+1,j+1,k))
450 END DO
451 END DO
452 DO j=jstr,jend
453 DO i=istr,iend
454 buoy2(i,j,k)=0.25_r8*(buoy2(i,j ,0)+buoy2(i-1,j ,0)+ &
455 & buoy2(i,j-1,0)+buoy2(i-1,j-1,0))
456 shear2(i,j,k)=0.25_r8*(shear2(i,j ,0)+shear2(i-1,j ,0)+ &
457 & shear2(i,j-1,0)+shear2(i-1,j-1,0))
458 END DO
459 END DO
460 END DO
461# endif
462
463
464
465
466
467
468
469
470
471
473# ifdef K_C2ADVECTION
474
475
476
477 DO j=jstr,jend
478 DO i=istr,iend+1
479 cff=0.25_r8*(huon(i,j,k)+huon(i,j,k+1))
480 fxk(i,j)=cff*(tke(i,j,k,3)+tke(i-1,j,k,3))
481 fxp(i,j)=cff*(gls(i,j,k,3)+gls(i-1,j,k,3))
482 END DO
483 END DO
484 DO j=jstr,jend+1
485 DO i=istr,iend
486 cff=0.25_r8*(hvom(i,j,k)+hvom(i,j,k+1))
487 fek(i,j)=cff*(tke(i,j,k,3)+tke(i,j-1,k,3))
488 fep(i,j)=cff*(gls(i,j,k,3)+gls(i,j-1,k,3))
489 END DO
490 END DO
491# else
492 DO j=jstr,jend
493 DO i=istrm1,iendp2
494 gradk(i,j)=(tke(i,j,k,3)-tke(i-1,j,k,3))
495# ifdef MASKING
496 gradk(i,j)=gradk(i,j)*umask(i,j)
497# endif
498 gradp(i,j)=(gls(i,j,k,3)-gls(i-1,j,k,3))
499# ifdef MASKING
500 gradp(i,j)=gradp(i,j)*umask(i,j)
501# endif
502 END DO
503 END DO
505 IF (
domain(ng)%Western_Edge(tile))
THEN
506 DO j=jstr,jend
507 gradk(istr-1,j)=gradk(istr,j)
508 gradp(istr-1,j)=gradp(istr,j)
509 END DO
510 END IF
511 END IF
513 IF (
domain(ng)%Eastern_Edge(tile))
THEN
514 DO j=jstr,jend
515 gradk(iend+2,j)=gradk(iend+1,j)
516 gradp(iend+2,j)=gradp(iend+1,j)
517 END DO
518 END IF
519 END IF
520# ifdef K_C4ADVECTION
521
522
523
524 cff1=1.0_r8/6.0_r8
525 DO j=jstr,jend
526 DO i=istr,iend+1
527 cff=0.5_r8*(huon(i,j,k)+huon(i,j,k+1))
528 fxk(i,j)=cff*0.5_r8*(tke(i-1,j,k,3)+tke(i,j,k,3)- &
529 & cff1*(gradk(i+1,j)-gradk(i-1,j)))
530 fxp(i,j)=cff*0.5_r8*(gls(i-1,j,k,3)+gls(i,j,k,3)- &
531 & cff1*(gradp(i+1,j)-gradp(i-1,j)))
532 END DO
533 END DO
534# else
535
536
537
538
539 DO j=jstr,jend
540 DO i=istr-1,iend+1
541 curvk(i,j)=gradk(i+1,j)-gradk(i,j)
542 curvp(i,j)=gradp(i+1,j)-gradp(i,j)
543 END DO
544 END DO
545 DO j=jstr,jend
546 DO i=istr,iend+1
547 cff=0.5_r8*(huon(i,j,k)+huon(i,j,k+1))
548 IF (cff.gt.0.0_r8) THEN
549 cff1=curvk(i-1,j)
550 cff2=curvp(i-1,j)
551 ELSE
552 cff1=curvk(i,j)
553 cff2=curvp(i,j)
554 END IF
555 fxk(i,j)=cff*0.5_r8*(tke(i-1,j,k,3)+tke(i,j,k,3)- &
556 & gadv*cff1)
557 fxp(i,j)=cff*0.5_r8*(gls(i-1,j,k,3)+gls(i,j,k,3)- &
558 & gadv*cff2)
559 END DO
560 END DO
561# endif
562 DO j=jstrm1,jendp2
563 DO i=istr,iend
564 gradk(i,j)=(tke(i,j,k,3)-tke(i,j-1,k,3))
565# ifdef MASKING
566 gradk(i,j)=gradk(i,j)*vmask(i,j)
567# endif
568 gradp(i,j)=(gls(i,j,k,3)-gls(i,j-1,k,3))
569# ifdef MASKING
570 gradp(i,j)=gradp(i,j)*vmask(i,j)
571# endif
572 END DO
573 END DO
575 IF (
domain(ng)%Southern_Edge(tile))
THEN
576 DO i=istr,iend
577 gradk(i,jstr-1)=gradk(i,jstr)
578 gradp(i,jstr-1)=gradp(i,jstr)
579 END DO
580 END IF
581 END IF
583 IF (
domain(ng)%Northern_Edge(tile))
THEN
584 DO i=istr,iend
585 gradk(i,jend+2)=gradk(i,jend+1)
586 gradp(i,jend+2)=gradp(i,jend+1)
587 END DO
588 END IF
589 END IF
590# ifdef K_C4ADVECTION
591 cff1=1.0_r8/6.0_r8
592 DO j=jstr,jend+1
593 DO i=istr,iend
594 cff=0.5_r8*(hvom(i,j,k)+hvom(i,j,k+1))
595 fek(i,j)=cff*0.5_r8*(tke(i,j-1,k,3)+tke(i,j,k,3)- &
596 & cff1*(gradk(i,j+1)-gradk(i,j-1)))
597 fep(i,j)=cff*0.5_r8*(gls(i,j-1,k,3)+gls(i,j,k,3)- &
598 & cff1*(gradp(i,j+1)-gradp(i,j-1)))
599 END DO
600 END DO
601# else
602 DO j=jstr-1,jend+1
603 DO i=istr,iend
604 curvk(i,j)=gradk(i,j+1)-gradk(i,j)
605 curvp(i,j)=gradp(i,j+1)-gradp(i,j)
606 END DO
607 END DO
608 DO j=jstr,jend+1
609 DO i=istr,iend
610 cff=0.5_r8*(hvom(i,j,k)+hvom(i,j,k+1))
611 IF (cff.gt.0.0_r8) THEN
612 cff1=curvk(i,j-1)
613 cff2=curvp(i,j-1)
614 ELSE
615 cff1=curvk(i,j)
616 cff2=curvp(i,j)
617 END IF
618 fek(i,j)=cff*0.5_r8*(tke(i,j-1,k,3)+tke(i,j,k,3)- &
619 & gadv*cff1)
620 fep(i,j)=cff*0.5_r8*(gls(i,j-1,k,3)+gls(i,j,k,3)- &
621 & gadv*cff2)
622 END DO
623 END DO
624# endif
625# endif
626
627
628
629 DO j=jstr,jend
630 DO i=istr,iend
631 cff=
dt(ng)*pm(i,j)*pn(i,j)
632 tke(i,j,k,nnew)=tke(i,j,k,nnew)- &
633 & cff*(fxk(i+1,j)-fxk(i,j)+ &
634 & fek(i,j+1)-fek(i,j))
635 tke(i,j,k,nnew)=max(tke(i,j,k,nnew),
gls_kmin(ng))
636 gls(i,j,k,nnew)=gls(i,j,k,nnew)- &
637 & cff*(fxp(i+1,j)-fxp(i,j)+ &
638 & fep(i,j+1)-fep(i,j))
639 gls(i,j,k,nnew)=max(gls(i,j,k,nnew),
gls_pmin(ng))
640 END DO
641 END DO
642 END DO
643
644
645
646 DO j=jstr,jend
647# ifdef K_C2ADVECTION
649 DO i=istr,iend
650 cff=0.25_r8*(w(i,j,k)+w(i,j,k-1))
651# ifdef WEC_VF
652 cff=cff+0.25_r8*(w_stokes(i,j,k)+w_stokes(i,j,k-1))
653# endif
654 fck(i,k)=cff*(tke(i,j,k,3)+tke(i,j,k-1,3))
655 fcp(i,k)=cff*(gls(i,j,k,3)+gls(i,j,k-1,3))
656 END DO
657 END DO
658# else
659 cff1=7.0_r8/12.0_r8
660 cff2=1.0_r8/12.0_r8
662 DO i=istr,iend
663 cff=0.5*(w(i,j,k)+w(i,j,k-1))
664# ifdef WEC_VF
665 cff=cff+0.5_r8*(w_stokes(i,j,k)+w_stokes(i,j,k-1))
666# endif
667 fck(i,k)=cff*(cff1*(tke(i,j,k-1,3)+ &
668 & tke(i,j,k ,3))- &
669 & cff2*(tke(i,j,k-2,3)+ &
670 & tke(i,j,k+1,3)))
671 fcp(i,k)=cff*(cff1*(gls(i,j,k-1,3)+ &
672 & gls(i,j,k ,3))- &
673 & cff2*(gls(i,j,k-2,3)+ &
674 & gls(i,j,k+1,3)))
675 END DO
676 END DO
677 cff1=1.0_r8/3.0_r8
678 cff2=5.0_r8/6.0_r8
679 cff3=1.0_r8/6.0_r8
680 DO i=istr,iend
681 cff=0.5_r8*(w(i,j,0)+w(i,j,1))
682# ifdef WEC_VF
683 cff=cff+0.5_r8*(w_stokes(i,j,0)+w_stokes(i,j,1))
684# endif
685 fck(i,1)=cff*(cff1*tke(i,j,0,3)+ &
686 & cff2*tke(i,j,1,3)- &
687 & cff3*tke(i,j,2,3))
688 fcp(i,1)=cff*(cff1*gls(i,j,0,3)+ &
689 & cff2*gls(i,j,1,3)- &
690 & cff3*gls(i,j,2,3))
691 cff=0.5_r8*(w(i,j,
n(ng))+w(i,j,
n(ng)-1))
692# ifdef WEC_VF
693 cff=cff+0.5_r8*(w_stokes(i,j,
n(ng))+w_stokes(i,j,
n(ng)-1))
694# endif
695 fck(i,
n(ng))=cff*(cff1*tke(i,j,
n(ng) ,3)+ &
696 & cff2*tke(i,j,
n(ng)-1,3)- &
697 & cff3*tke(i,j,
n(ng)-2,3))
698 fcp(i,
n(ng))=cff*(cff1*gls(i,j,
n(ng) ,3)+ &
699 & cff2*gls(i,j,
n(ng)-1,3)- &
700 & cff3*gls(i,j,
n(ng)-2,3))
701 END DO
702# endif
703
704
705
707 DO i=istr,iend
708 cff=
dt(ng)*pm(i,j)*pn(i,j)
709 tke(i,j,k,nnew)=tke(i,j,k,nnew)- &
710 & cff*(fck(i,k+1)-fck(i,k))
711 tke(i,j,k,nnew)=max(tke(i,j,k,nnew),
gls_kmin(ng))
712 gls(i,j,k,nnew)=gls(i,j,k,nnew)- &
713 & cff*(fcp(i,k+1)-fcp(i,k))
714 gls(i,j,k,nnew)=max(gls(i,j,k,nnew),
gls_pmin(ng))
715 END DO
716 END DO
717
718
719
720
721
722
723
724
726 DO i=istr,iend
728 fck(i,k)=cff*(akk(i,j,k)+akk(i,j,k-1))/hz(i,j,k)
729 fcp(i,k)=cff*(akp(i,j,k)+akp(i,j,k-1))/hz(i,j,k)
730 cf(i,k)=0.0_r8
731 END DO
732 fcp(i,1)=0.0_r8
734 fck(i,1)=0.0_r8
736 END DO
737
738
739
740 DO i=istr,iend
742
743
744
745
746 strat2=buoy2(i,j,k)
747 IF (strat2.gt.0.0_r8) THEN
749 ELSE
751 END IF
752 kprod=shear2(i,j,k)*(akv(i,j,k)-
akv_bak(ng))- &
756
757
758
759
760 cff1=1.0_r8
761 IF (kprod.lt.0.0_r8) THEN
763 cff1=0.0_r8
764 END IF
765 cff2=1.0_r8
766 IF (pprod.lt.0.0_r8) THEN
767 pprod=pprod+gls_c3*strat2*(akt(i,j,k,
itemp)- &
769 cff2=0.0_r8
770 END IF
771
772
773
774 cff=0.5_r8*(hz(i,j,k)+hz(i,j,k+1))
775 tke(i,j,k,nnew)=tke(i,j,k,nnew)+ &
777 gls(i,j,k,nnew)=gls(i,j,k,nnew)+ &
778 &
dt(ng)*cff*pprod*gls(i,j,k,nstp)/ &
780
781
782
783 wall_fac=1.0_r8
784 IF (lmy25) THEN
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
812 & (gls(i,j,k,nstp)**( gls_exp1)*cmu_fac1* &
813 & tke(i,j,k,nstp)**(-tke_exp1)* &
814 & (1.0_r8/ (z_w(i,j,k)-z_w(i,j,0))))**2+ &
816 & (gls(i,j,k,nstp)**( gls_exp1)*cmu_fac1* &
817 & tke(i,j,k,nstp)**(-tke_exp1)* &
818 & (1.0_r8/ (z_w(i,j,
n(ng))-z_w(i,j,k))))**2
819 END IF
820
821 bck(i,k)=cff*(1.0_r8+
dt(ng)* &
822 & gls(i,j,k,nstp)**(-gls_exp1)*cmu_fac2* &
823 & tke(i,j,k,nstp)**( tke_exp2)+ &
824 &
dt(ng)*(1.0_r8-cff1)*strat2* &
826 & tke(i,j,k,nstp))- &
827 & fck(i,k)-fck(i,k+1)
828 bcp(i,k)=cff*(1.0_r8+
dt(ng)*
gls_c2(ng)*wall_fac* &
829 & gls(i,j,k,nstp)**(-gls_exp1)*cmu_fac2* &
830 & tke(i,j,k,nstp)**( tke_exp2)+ &
831 &
dt(ng)*(1.0_r8-cff2)*gls_c3*strat2* &
833 & tke(i,j,k,nstp))- &
834 & fcp(i,k)-fcp(i,k+1)
835 END DO
836 END DO
837
838
839
840
841
842
843
844
845
846 DO i=istr,iend
847# if defined CRAIG_BANNER
848 tke(i,j,
n(ng),nnew)=max(cmu_fac4*0.5_r8* &
849 & sqrt((sustr(i,j)+sustr(i+1,j))**2+ &
850 & (svstr(i,j)+svstr(i,j+1))**2)* &
853# elif defined TKE_WAVEDISS
854 tke(i,j,
n(ng),nnew)=max(cmu_fac4*(
sz_alpha(ng)* &
855 & (dissip_break(i,j)+ &
856 & dissip_wcap(i,j)))** &
858# else
859 tke(i,j,
n(ng),nnew)=max(cmu_fac3*0.5_r8* &
860 & sqrt((sustr(i,j)+sustr(i+1,j))**2+ &
861 & (svstr(i,j)+svstr(i,j+1))**2), &
863# endif
864 tke(i,j,0,nnew)=max(cmu_fac3*0.5_r8* &
865 & sqrt((bustr(i,j)+bustr(i+1,j))**2+ &
866 & (bvstr(i,j)+bvstr(i,j+1))**2), &
868# if defined CHARNOK
870 & sqrt((sustr(i,j)+sustr(i+1,j))**2+ &
871 & (svstr(i,j)+svstr(i,j+1))**2), &
872 & zos_min)
873# elif defined ZOS_HSIG
875# else
876 zos_eff(i)=zos_min
877# endif
879 & tke(i,j,
n(ng),nnew)**
gls_m(ng)* &
880 & (l_sft*zos_eff(i))**
gls_n(ng), &
883 gls(i,j,0,nnew)=max(cff*tke(i,j,0,nnew)**(
gls_m(ng)), &
885 END DO
886
887
888
889 DO i=istr,iend
890# if defined CRAIG_BANNER
892 & (0.50_r8* &
893 & sqrt((sustr(i,j)+sustr(i+1,j))**2+ &
894 & (svstr(i,j)+svstr(i,j+1))**2))**1.5_r8
895# elif defined TKE_WAVEDISS
896 tke_fluxt(i)=
dt(ng)*
sz_alpha(ng)*(dissip_break(i,j)+ &
897 & dissip_wcap(i,j))
898# else
899 tke_fluxt(i)=0.0_r8
900# endif
901 tke_fluxb(i)=0.0_r8
902
903 cff=1.0_r8/bck(i,
n(ng)-1)
904 cf(i,
n(ng)-1)=cff*fck(i,
n(ng)-1)
905 tke(i,j,
n(ng)-1,nnew)=cff*(tke(i,j,
n(ng)-1,nnew)+tke_fluxt(i))
906 END DO
907 DO i=istr,iend
909 cff=1.0_r8/(bck(i,k)-cf(i,k+1)*fck(i,k+1))
910 cf(i,k)=cff*fck(i,k)
911 tke(i,j,k,nnew)=cff*(tke(i,j,k,nnew)- &
912 & fck(i,k+1)*tke(i,j,k+1,nnew))
913 END DO
914 tke(i,j,1,nnew)=tke(i,j,1,nnew)-cff*tke_fluxb(i)
915 END DO
917 DO i=istr,iend
918 tke(i,j,k,nnew)=tke(i,j,k,nnew)-cf(i,k)*tke(i,j,k-1,nnew)
919 END DO
920 END DO
921
922
923
924 DO i=istr,iend
925 cff=0.5_r8*(tke(i,j,
n(ng),nnew)+tke(i,j,
n(ng)-1,nnew))
926 gls_fluxt(i)=
dt(ng)*gls_fac3*cff**
gls_m(ng)* &
927 & l_sft**(
gls_n(ng))* &
928 & (zos_eff(i)+0.5_r8*hz(i,j,
n(ng)))** &
929 & (
gls_n(ng)-1.0_r8)* &
930 & 0.5_r8*(akp(i,j,
n(ng))+akp(i,j,
n(ng)-1))
931# ifdef CRAIG_BANNER
932 gls_fluxt(i)=gls_fluxt(i)-
dt(ng)* &
934 & cff**(
gls_m(ng)-1.0_r8)* &
935 & ((zos_eff(i)+0.5_r8*hz(i,j,
n(ng)))*l_sft)** &
938 & (0.5_r8* &
939 & sqrt((sustr(i,j)+sustr(i+1,j))**2+ &
940 & (svstr(i,j)+svstr(i,j+1))**2))**1.5_r8
941# elif defined TKE_WAVEDISS
942 gls_fluxt(i)=gls_fluxt(i)-
dt(ng)* &
944 & cff**(
gls_m(ng)-1.0_r8)* &
945 & ((zos_eff(i)+0.5_r8*hz(i,j,
n(ng)))*l_sft)** &
948 & (dissip_break(i,j)+dissip_wcap(i,j))
949# endif
950 cff=0.5_r8*(tke(i,j,0,nnew)+tke(i,j,1,nnew))
951 gls_fluxb(i)=
dt(ng)*gls_fac2*(cff**
gls_m(ng))* &
952 & (0.5_r8*hz(i,j,1)+zob_min(i,j))** &
953 & (
gls_n(ng)-1.0_r8)* &
954 & 0.5_r8*(akp(i,j,0)+akp(i,j,1))
955
956 cff=1.0_r8/bcp(i,
n(ng)-1)
957 cf(i,
n(ng)-1)=cff*fcp(i,
n(ng)-1)
958 gls(i,j,
n(ng)-1,nnew)=cff*(gls(i,j,
n(ng)-1,nnew)-gls_fluxt(i))
959 END DO
960 DO i=istr,iend
962 cff=1.0_r8/(bcp(i,k)-cf(i,k+1)*fcp(i,k+1))
963 cf(i,k)=cff*fcp(i,k)
964 gls(i,j,k,nnew)=cff*(gls(i,j,k,nnew)- &
965 & fcp(i,k+1)*gls(i,j,k+1,nnew))
966 END DO
967 gls(i,j,1,nnew)=gls(i,j,1,nnew)-cff*gls_fluxb(i)
968
969 END DO
971 DO i=istr,iend
972 gls(i,j,k,nnew)=gls(i,j,k,nnew)-cf(i,k)*gls(i,j,k-1,nnew)
973
974 END DO
975 END DO
976
977
978
979
980
981 DO i=istr,iend
983
984
985
986 tke(i,j,k,nnew)=max(tke(i,j,k,nnew),
gls_kmin(ng))
987 gls(i,j,k,nnew)=max(gls(i,j,k,nnew),
gls_pmin(ng))
988 IF (
gls_n(ng).ge.0.0_r8)
THEN
989 gls(i,j,k,nnew)=min(gls(i,j,k,nnew),gls_fac5* &
990 & tke(i,j,k,nnew)**(tke_exp4)* &
991 & (sqrt(max(0.0_r8, &
992 & buoy2(i,j,k)))+eps)** &
994 ELSE
995 gls(i,j,k,nnew)=max(gls(i,j,k,nnew),gls_fac5* &
996 & tke(i,j,k,nnew)**(tke_exp4)* &
997 & (sqrt(max(0.0_r8, &
998 & buoy2(i,j,k)))+eps)** &
1000 END IF
1001 ls_unlmt=max(eps, &
1002 & gls(i,j,k,nnew)**( gls_exp1)*cmu_fac1* &
1003 & tke(i,j,k,nnew)**(-tke_exp1))
1004 IF (buoy2(i,j,k).gt.0.0_r8) THEN
1005 ls_lmt=min(ls_unlmt, &
1006 & sqrt(0.56_r8*tke(i,j,k,nnew)/ &
1007 & (max(0.0_r8,buoy2(i,j,k))+eps)))
1008 ELSE
1009 ls_lmt=ls_unlmt
1010 END IF
1011
1012
1013
1015 & tke(i,j,k,nnew)**
gls_m(ng)* &
1017
1018
1019
1020
1021 gh=min(
gls_gh0,-buoy2(i,j,k)*ls_lmt*ls_lmt/ &
1022 & (2.0_r8*tke(i,j,k,nnew)))
1026# if defined CANUTO_A || defined CANUTO_B
1027
1028
1029
1032 gm=min(gm,shear2(i,j,k)*ls_lmt*ls_lmt/ &
1033 & (2.0_r8*tke(i,j,k,nnew)))
1034
1035
1036
1037
1040 &
gls_b5*gls_fac6**2*gm*gm
1043 sm=max(sm,0.0_r8)
1044 sh=max(sh,0.0_r8)
1045
1046
1047
1050# elif defined KANTHA_CLAYSON
1054# else /* Galperin */
1058# endif
1059
1060
1061
1062
1063
1064 ql=sqrt2*0.5_r8*(ls_lmt*sqrt(tke(i,j,k,nnew))+ &
1065 & lscale(i,j,k)*sqrt(tke(i,j,k,nstp)))
1068 akt(i,j,k,itrc)=
akt_bak(itrc,ng)+sh*ql
1069 END DO
1070
1071
1072
1073
1076
1077# if defined CRAIG_BANNER || defined TKE_WAVEDISS
1078
1079
1080
1081
1082 pprod=
gls_c1(ng)*shear2(i,j,k)*akv(i,j,k)
1083 cff=cmu_fac2*tke(i,j,k,nnew)**(1.5_r8+tke_exp1)* &
1084 & gls(i,j,k,nnew)**(-1.0_r8/
gls_n(ng))
1085 cff2=min(pprod/cff, 1.0_r8)
1086 sig_eff=cff2*
gls_sigp(ng)+(1.0_r8-cff2)*gls_sigp_cb
1087 akp(i,j,k)=
akp_bak(ng)+sm*ql/sig_eff
1088# else
1089 akp(i,j,k)=
akp_bak(ng)+sm*ql*ogls_sigp
1090# endif
1091
1092
1093
1094 lscale(i,j,k)=ls_lmt
1095 END DO
1096
1097
1098
1099
1101 & sqrt(tke(i,j,
n(ng),nnew))
1103 & sqrt(tke(i,j,0,nnew))
1104
1107 akp(i,j,
n(ng))=
akp_bak(ng)+akv(i,j,
n(ng))*ogls_sigp
1109
1111 akt(i,j,
n(ng),itrc)=
akt_bak(itrc,ng)
1112 akt(i,j,0,itrc)=
akt_bak(itrc,ng)
1113 END DO
1114
1115# if defined LIMIT_VDIFF || defined LIMIT_VVISC
1116
1117
1118
1119
1120
1121
1123# ifdef LIMIT_VDIFF
1125 akt(i,j,k,itrc)=min(
akt_limit(itrc,ng), akt(i,j,k,itrc))
1126 END DO
1127# endif
1128# ifdef LIMIT_VVISC
1129 akv(i,j,k)=min(
akv_limit(ng), akv(i,j,k))
1130# endif
1131 END DO
1132# endif
1133 END DO
1134 END DO
1135
1136
1137
1138
1139
1141 IF (
domain(ng)%Western_Edge(tile))
THEN
1142 DO j=jstr,jend
1144 akt(istr-1,j,k,itrc)=akt(istr,j,k,itrc)
1145 END DO
1146 akv(istr-1,j,k)=akv(istr,j,k)
1147 END DO
1148 END IF
1149 IF (
domain(ng)%Eastern_Edge(tile))
THEN
1150 DO j=jstr,jend
1152 akt(iend+1,j,k,itrc)=akt(iend,j,k,itrc)
1153 END DO
1154 akv(iend+1,j,k)=akv(iend,j,k)
1155 END DO
1156
1157 END IF
1158 IF (
domain(ng)%Southern_Edge(tile))
THEN
1159 DO i=istr,iend
1161 akt(i,jstr-1,k,itrc)=akt(i,jstr,k,itrc)
1162 END DO
1163 akv(i,jstr-1,k)=akv(i,jstr,k)
1164 END DO
1165 END IF
1166 IF (
domain(ng)%Northern_Edge(tile))
THEN
1167 DO i=istr,iend
1169 akt(i,jend+1,k,itrc)=akt(i,jend,k,itrc)
1170 END DO
1171 akv(i,jend+1,k)=akv(i,jend,k)
1172 END DO
1173 END IF
1174 IF (
domain(ng)%SouthWest_Corner(tile))
THEN
1176 akt(istr-1,jstr-1,k,itrc)=0.5_r8* &
1177 & (akt(istr ,jstr-1,k,itrc)+ &
1178 & akt(istr-1,jstr ,k,itrc))
1179 END DO
1180 akv(istr-1,jstr-1,k)=0.5_r8* &
1181 & (akv(istr ,jstr-1,k)+ &
1182 & akv(istr-1,jstr ,k))
1183 END IF
1184 IF (
domain(ng)%SouthEast_Corner(tile))
THEN
1186 akt(iend+1,jstr-1,k,itrc)=0.5_r8* &
1187 & (akt(iend ,jstr-1,k,itrc)+ &
1188 & akt(iend+1,jstr ,k,itrc))
1189 END DO
1190 akv(iend+1,jstr-1,k)=0.5_r8* &
1191 & (akv(iend ,jstr-1,k)+ &
1192 & akv(iend+1,jstr ,k))
1193 END IF
1194 IF (
domain(ng)%NorthWest_Corner(tile))
THEN
1196 akt(istr-1,jend+1,k,itrc)=0.5_r8* &
1197 & (akt(istr ,jend+1,k,itrc)+ &
1198 & akt(istr-1,jend ,k,itrc))
1199 END DO
1200 akv(istr-1,jend+1,k)=0.5_r8* &
1201 & (akv(istr ,jend+1,k)+ &
1202 & akv(istr-1,jend ,k))
1203 END IF
1204 IF (
domain(ng)%NorthEast_Corner(tile))
THEN
1206 akt(iend+1,jend+1,k,itrc)=0.5_r8* &
1207 & (akt(iend ,jend+1,k,itrc)+ &
1208 & akt(iend+1,jend ,k,itrc))
1209 END DO
1210 akv(iend+1,jend+1,k)=0.5_r8* &
1211 & (akv(iend ,jend+1,k)+ &
1212 & akv(iend+1,jend ,k))
1213 END IF
1214 END DO
1215
1217 & lbi, ubi, lbj, ubj,
n(ng), &
1218 & imins, imaxs, jmins, jmaxs, &
1219 & nnew, nstp, &
1220 & gls, tke)
1221
1224 & lbi, ubi, lbj, ubj, 0,
n(ng), &
1225 & tke(:,:,:,nnew))
1227 & lbi, ubi, lbj, ubj, 0,
n(ng), &
1228 & gls(:,:,:,nnew))
1230 & lbi, ubi, lbj, ubj, 0,
n(ng), &
1231 & akv)
1234 & lbi, ubi, lbj, ubj, 0,
n(ng), &
1235 & akt(:,:,:,itrc))
1236 END DO
1237 END IF
1238
1239# ifdef DISTRIBUTE
1241 & lbi, ubi, lbj, ubj, 0,
n(ng), &
1244 & tke(:,:,:,nnew), &
1245 & gls(:,:,:,nnew), &
1246 & akv)
1248 & lbi, ubi, lbj, ubj, 0,
n(ng), 1,
nat, &
1251 & akt)
1252# endif
1253
1254 RETURN
subroutine exchange_w3d_tile(ng, tile, lbi, ubi, lbj, ubj, lbk, ubk, a)
integer, dimension(:), allocatable n
type(t_domain), dimension(:), allocatable domain
integer, dimension(:), allocatable lm
integer, dimension(:), allocatable mm
real(r8), dimension(:), allocatable gls_p
real(r8), dimension(:), allocatable gls_n
real(r8), dimension(:), allocatable gls_m
real(r8), dimension(:), allocatable zos_hsig_alpha
real(r8), parameter gls_ghcri
real(dp), dimension(:), allocatable dt
real(r8), dimension(:), allocatable sz_alpha
logical, dimension(:), allocatable ewperiodic
logical, dimension(:), allocatable nsperiodic
real(r8), dimension(:), allocatable gls_c2
real(r8), dimension(:), allocatable charnok_alpha
real(r8), dimension(:), allocatable zos
real(r8), dimension(:), allocatable akk_bak
real(r8), dimension(:), allocatable gls_cmu0
real(r8), parameter gls_ghmin
real(r8), dimension(:), allocatable gls_sigk
real(r8), dimension(:), allocatable crgban_cw
real(r8), dimension(:,:), allocatable akt_bak
real(r8), dimension(:), allocatable akv_bak
real(r8), dimension(:), allocatable gls_pmin
real(r8), parameter gls_e2
logical, dimension(:,:), allocatable compositegrid
integer, parameter isouth
real(r8), dimension(:), allocatable akp_bak
real(r8), dimension(:,:), allocatable akt_limit
real(r8), dimension(:), allocatable gls_c1
real(r8), dimension(:), allocatable akv_limit
real(r8), parameter gls_gh0
real(r8), dimension(:), allocatable gls_sigp
real(r8), dimension(:), allocatable gls_kmin
integer, parameter inorth
real(r8), dimension(:), allocatable gls_c3m
real(r8), dimension(:), allocatable gls_c3p
subroutine mp_exchange4d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, lbk, ubk, lbt, ubt, nghost, ew_periodic, ns_periodic, a, b, c)
subroutine mp_exchange3d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, lbk, ubk, nghost, ew_periodic, ns_periodic, a, b, c, d)
subroutine, public tkebc_tile(ng, tile, lbi, ubi, lbj, ubj, ubk, imins, imaxs, jmins, jmaxs, nout, nstp, gls, tke)