98 |
_RL a(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
_RL a(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
99 |
_RL b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
_RL b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
100 |
_RL c(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
_RL c(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) |
101 |
c INTEGER errCode |
INTEGER errCode |
102 |
#ifdef ALLOW_GGL90_HORIZDIFF |
#ifdef ALLOW_GGL90_HORIZDIFF |
103 |
C- xA, yA - area of lateral faces |
C- xA, yA - area of lateral faces |
104 |
C- dfx, dfy - diffusive flux across lateral faces |
C- dfx, dfy - diffusive flux across lateral faces |
414 |
ENDDO |
ENDDO |
415 |
ENDDO |
ENDDO |
416 |
C Compute divergence of fluxes |
C Compute divergence of fluxes |
|
C- jmc: concerned about missing a deltaT since gTKE is already the future TKE. |
|
417 |
DO j=1-Oly,sNy+Oly-1 |
DO j=1-Oly,sNy+Oly-1 |
418 |
DO i=1-Olx,sNx+Olx-1 |
DO i=1-Olx,sNx+Olx-1 |
419 |
gTKE(i,j,k)=gTKE(i,j,k) |
gTKE(i,j,k)=gTKE(i,j,k) |
420 |
& -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj) |
& -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj) |
421 |
& *( (dfx(i+1,j)-dfx(i,j)) |
& *( (dfx(i+1,j)-dfx(i,j)) |
422 |
& +(dfy(i,j+1)-dfy(i,j)) |
& +(dfy(i,j+1)-dfy(i,j)) |
423 |
& ) |
& )*deltaTggl90 |
424 |
ENDDO |
ENDDO |
425 |
ENDDO |
ENDDO |
426 |
C end of k-loop |
C end of k-loop |
441 |
ENDDO |
ENDDO |
442 |
ENDDO |
ENDDO |
443 |
DO k=2,Nr |
DO k=2,Nr |
444 |
C- jmc: concerned that a(k=2) should always be zero |
km1=MAX(2,k-1) |
|
C and would be better without recip_hFacC |
|
|
km1=max(2,k-1) |
|
445 |
DO j=jMin,jMax |
DO j=jMin,jMax |
446 |
DO i=iMin,iMax |
DO i=iMin,iMax |
447 |
|
C- We keep recip_hFacC in the diffusive flux calculation, |
448 |
|
C- but no hFacC in TKE volume control |
449 |
|
C- No need for maskC(k-1) with recip_hFacC(k-1) |
450 |
a(i,j,k) = -deltaTggl90 |
a(i,j,k) = -deltaTggl90 |
|
c & *recip_drF(km1)*recip_hFacI(i,j,k,bi,bj) |
|
451 |
& *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj) |
& *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj) |
452 |
& *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1)) |
& *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1)) |
453 |
& *recip_drC(k)*maskC(i,j,k,bi,bj)*maskC(i,j,k-1,bi,bj) |
& *recip_drC(k)*maskC(i,j,k,bi,bj) |
454 |
ENDDO |
ENDDO |
455 |
ENDDO |
ENDDO |
456 |
ENDDO |
ENDDO |
461 |
ENDDO |
ENDDO |
462 |
ENDDO |
ENDDO |
463 |
DO k=2,Nr |
DO k=2,Nr |
|
C- jmc: concerned that c(k) from k=kLow to k=Nr should always be zero |
|
464 |
DO j=jMin,jMax |
DO j=jMin,jMax |
465 |
DO i=iMin,iMax |
DO i=iMin,iMax |
466 |
C- jmc: this is dangerous since klowC could be zero if land column |
kp1=MAX(1,MIN(klowC(i,j,bi,bj),k+1)) |
467 |
C and would be better without recip_hFacC |
C- We keep recip_hFacC in the diffusive flux calculation, |
468 |
kp1=min(klowC(i,j,bi,bj),k+1) |
C- but no hFacC in TKE volume control |
469 |
|
C- No need for maskC(k) with recip_hFacC(k) |
470 |
c(i,j,k) = -deltaTggl90 |
c(i,j,k) = -deltaTggl90 |
|
c & *recip_drF( k )*recip_hFacI(i,j,k,bi,bj) |
|
471 |
& *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj) |
& *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj) |
472 |
& *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1)) |
& *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1)) |
473 |
& *recip_drC(k)*maskC(i,j,k,bi,bj)*maskC(i,j,k-1,bi,bj) |
& *recip_drC(k)*maskC(i,j,k-1,bi,bj) |
474 |
ENDDO |
ENDDO |
475 |
ENDDO |
ENDDO |
476 |
ENDDO |
ENDDO |
477 |
C-- Center diagonal |
C-- Center diagonal |
478 |
DO k=1,Nr |
DO k=1,Nr |
479 |
|
km1 = MAX(k-1,1) |
480 |
DO j=jMin,jMax |
DO j=jMin,jMax |
481 |
DO i=iMin,iMax |
DO i=iMin,iMax |
482 |
b(i,j,k) = 1. _d 0 - c(i,j,k) - a(i,j,k) |
b(i,j,k) = 1. _d 0 - c(i,j,k) - a(i,j,k) |
483 |
& + ab15*deltaTggl90*GGL90ceps*SQRTTKE(I,J,K) |
& + ab15*deltaTggl90*GGL90ceps*SQRTTKE(I,J,K) |
484 |
& *rMixingLength(I,J,K)*maskC(i,j,k,bi,bj) |
& * rMixingLength(I,J,K) |
485 |
|
& * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj) |
486 |
ENDDO |
ENDDO |
487 |
ENDDO |
ENDDO |
488 |
ENDDO |
ENDDO |
489 |
C end set up matrix |
C end set up matrix |
490 |
|
|
491 |
C Apply boundary condition |
C Apply boundary condition |
492 |
C- jmc: concerned about conservation when a or c are changed after computing b |
kp1 = MIN(Nr,kSurf+1) |
493 |
DO J=jMin,jMax |
DO J=jMin,jMax |
494 |
DO I=iMin,iMax |
DO I=iMin,iMax |
495 |
C estimate friction velocity uStar from surface forcing |
C estimate friction velocity uStar from surface forcing |
500 |
& + surfaceForcingV(I, J+1,bi,bj) ) )**2 |
& + surfaceForcingV(I, J+1,bi,bj) ) )**2 |
501 |
& ) |
& ) |
502 |
C Dirichlet surface boundary condition for TKE |
C Dirichlet surface boundary condition for TKE |
|
C- jmc: would be much better to update the provisional TKE (i.e. gTKE) at k=2 |
|
503 |
gTKE(I,J,kSurf) = MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare) |
gTKE(I,J,kSurf) = MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare) |
504 |
& *maskC(I,J,kSurf,bi,bj) |
& *maskC(I,J,kSurf,bi,bj) |
505 |
|
gTKE(i,j,kp1) = gTKE(i,j,kp1) |
506 |
|
& - a(i,j,kp1)*gTKE(i,j,kSurf) |
507 |
|
a(i,j,kp1) = 0. _d 0 |
508 |
C Dirichlet bottom boundary condition for TKE = GGL90TKEbottom |
C Dirichlet bottom boundary condition for TKE = GGL90TKEbottom |
509 |
kBottom = MAX(kLowC(I,J,bi,bj),1) |
kBottom = MAX(kLowC(I,J,bi,bj),1) |
510 |
gTKE(I,J,kBottom) = gTKE(I,J,kBottom) |
gTKE(I,J,kBottom) = gTKE(I,J,kBottom) |
514 |
ENDDO |
ENDDO |
515 |
|
|
516 |
C solve tri-diagonal system, and store solution on gTKE (previously rhs) |
C solve tri-diagonal system, and store solution on gTKE (previously rhs) |
517 |
c CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax, |
CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax, |
518 |
c I a, b, c, |
I a, b, c, |
519 |
c U gTKE, |
U gTKE, |
520 |
c O errCode, |
O errCode, |
521 |
c I bi, bj, myThid ) |
I bi, bj, myThid ) |
522 |
CALL GGL90_SOLVE( bi, bj, iMin, iMax, jMin, jMax, |
c CALL GGL90_SOLVE( bi, bj, iMin, iMax, jMin, jMax, |
523 |
I a, b, c, |
c I a, b, c, |
524 |
U gTKE, |
c U gTKE, |
525 |
I myThid ) |
c I myThid ) |
526 |
|
|
527 |
C now update TKE |
C now update TKE |
528 |
DO K=1,Nr |
DO K=1,Nr |