27 |
C GGL90TKE :: sub-grid turbulent kinetic energy (m^2/s^2) |
C GGL90TKE :: sub-grid turbulent kinetic energy (m^2/s^2) |
28 |
C GGL90viscAz :: GGL90 eddy viscosity coefficient (m^2/s) |
C GGL90viscAz :: GGL90 eddy viscosity coefficient (m^2/s) |
29 |
C GGL90diffKzT :: GGL90 diffusion coefficient for temperature (m^2/s) |
C GGL90diffKzT :: GGL90 diffusion coefficient for temperature (m^2/s) |
|
C |
|
30 |
C \ev |
C \ev |
31 |
|
|
32 |
C !USES: ============================================================ |
C !USES: ============================================================ |
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 |
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 |
136 |
KappaE(I,J,K) = 0. _d 0 |
KappaE(I,J,K) = 0. _d 0 |
137 |
TKEPrandtlNumber(I,J,K) = 1. _d 0 |
TKEPrandtlNumber(I,J,K) = 1. _d 0 |
138 |
GGL90mixingLength(I,J,K) = GGL90mixingLengthMin |
GGL90mixingLength(I,J,K) = GGL90mixingLengthMin |
|
rMixingLength(I,J,K) = 0. _d 0 |
|
139 |
ENDDO |
ENDDO |
140 |
ENDDO |
ENDDO |
141 |
ENDDO |
ENDDO |
144 |
rhoK(I,J) = 0. _d 0 |
rhoK(I,J) = 0. _d 0 |
145 |
rhoKm1(I,J) = 0. _d 0 |
rhoKm1(I,J) = 0. _d 0 |
146 |
totalDepth(I,J) = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj) |
totalDepth(I,J) = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj) |
147 |
|
rMixingLength(i,j,1) = 0. _d 0 |
148 |
mxLength_Dn(I,J,1) = GGL90mixingLengthMin |
mxLength_Dn(I,J,1) = GGL90mixingLengthMin |
149 |
|
SQRTTKE(i,j,1) = SQRT( GGL90TKE(i,j,1,bi,bj) ) |
150 |
ENDDO |
ENDDO |
151 |
ENDDO |
ENDDO |
152 |
|
|
168 |
DO J=jMin,jMax |
DO J=jMin,jMax |
169 |
DO I=iMin,iMax |
DO I=iMin,iMax |
170 |
SQRTTKE(i,j,k)=SQRT( GGL90TKE(I,J,K,bi,bj) ) |
SQRTTKE(i,j,k)=SQRT( GGL90TKE(I,J,K,bi,bj) ) |
171 |
C |
|
172 |
C buoyancy frequency |
C buoyancy frequency |
|
C |
|
173 |
Nsquare(i,j,k) = - gravity*recip_rhoConst*recip_drC(K) |
Nsquare(i,j,k) = - gravity*recip_rhoConst*recip_drC(K) |
174 |
& * ( rhoKm1(I,J) - rhoK(I,J) )*maskC(I,J,K,bi,bj) |
& * ( rhoKm1(I,J) - rhoK(I,J) )*maskC(I,J,K,bi,bj) |
175 |
cC vertical shear term (dU/dz)^2+(dV/dz)^2 |
cC vertical shear term (dU/dz)^2+(dV/dz)^2 |
303 |
DO I=iMin,iMax |
DO I=iMin,iMax |
304 |
GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K), |
GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K), |
305 |
& mxLength_Dn(I,J,K)) |
& mxLength_Dn(I,J,K)) |
306 |
tmpmlx = SQRT(GGL90mixingLength(I,J,K) * mxLength_Dn(I,J,K)) |
tmpmlx = SQRT( GGL90mixingLength(I,J,K)*mxLength_Dn(I,J,K) ) |
307 |
tmpmlx = MAX( tmpmlx, GGL90mixingLengthMin) |
tmpmlx = MAX( tmpmlx, GGL90mixingLengthMin) |
308 |
rMixingLength(I,J,K) = 1. _d 0 / tmpmlx |
rMixingLength(I,J,K) = 1. _d 0 / tmpmlx |
309 |
ENDDO |
ENDDO |
414 |
ENDDO |
ENDDO |
415 |
ENDDO |
ENDDO |
416 |
C Compute divergence of fluxes |
C Compute divergence of fluxes |
417 |
|
C- jmc: concerned about missing a deltaT since gTKE is already the future TKE. |
418 |
DO j=1-Oly,sNy+Oly-1 |
DO j=1-Oly,sNy+Oly-1 |
419 |
DO i=1-Olx,sNx+Olx-1 |
DO i=1-Olx,sNx+Olx-1 |
420 |
gTKE(i,j,k)=gTKE(i,j,k) |
gTKE(i,j,k)=gTKE(i,j,k) |
442 |
ENDDO |
ENDDO |
443 |
ENDDO |
ENDDO |
444 |
DO k=2,Nr |
DO k=2,Nr |
445 |
|
C- jmc: concerned that a(k=2) should always be zero |
446 |
|
C and would be better without recip_hFacC |
447 |
km1=max(2,k-1) |
km1=max(2,k-1) |
448 |
DO j=jMin,jMax |
DO j=jMin,jMax |
449 |
DO i=iMin,iMax |
DO i=iMin,iMax |
462 |
ENDDO |
ENDDO |
463 |
ENDDO |
ENDDO |
464 |
DO k=2,Nr |
DO k=2,Nr |
465 |
|
C- jmc: concerned that c(k) from k=kLow to k=Nr should always be zero |
466 |
DO j=jMin,jMax |
DO j=jMin,jMax |
467 |
DO i=iMin,iMax |
DO i=iMin,iMax |
468 |
|
C- jmc: this is dangerous since klowC could be zero if land column |
469 |
|
C and would be better without recip_hFacC |
470 |
kp1=min(klowC(i,j,bi,bj),k+1) |
kp1=min(klowC(i,j,bi,bj),k+1) |
471 |
c(i,j,k) = -deltaTggl90 |
c(i,j,k) = -deltaTggl90 |
472 |
c & *recip_drF( k )*recip_hFacI(i,j,k,bi,bj) |
c & *recip_drF( k )*recip_hFacI(i,j,k,bi,bj) |
488 |
ENDDO |
ENDDO |
489 |
C end set up matrix |
C end set up matrix |
490 |
|
|
|
C |
|
491 |
C Apply boundary condition |
C Apply boundary condition |
492 |
C |
C- jmc: concerned about conservation when a or c are changed after computing b |
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 |
503 |
|
C- jmc: would be much better to update the provisional TKE (i.e. gTKE) at k=2 |
504 |
gTKE(I,J,kSurf) = MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare) |
gTKE(I,J,kSurf) = MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare) |
505 |
& *maskC(I,J,kSurf,bi,bj) |
& *maskC(I,J,kSurf,bi,bj) |
506 |
C Dirichlet bottom boundary condition for TKE = GGL90TKEbottom |
C Dirichlet bottom boundary condition for TKE = GGL90TKEbottom |
510 |
c(I,J,kBottom) = 0. _d 0 |
c(I,J,kBottom) = 0. _d 0 |
511 |
ENDDO |
ENDDO |
512 |
ENDDO |
ENDDO |
513 |
C |
|
514 |
C solve tri-diagonal system, and store solution on gTKE (previously rhs) |
C solve tri-diagonal system, and store solution on gTKE (previously rhs) |
515 |
C |
c CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax, |
516 |
|
c I a, b, c, |
517 |
|
c U gTKE, |
518 |
|
c O errCode, |
519 |
|
c I bi, bj, myThid ) |
520 |
CALL GGL90_SOLVE( bi, bj, iMin, iMax, jMin, jMax, |
CALL GGL90_SOLVE( bi, bj, iMin, iMax, jMin, jMax, |
521 |
I a, b, c, |
I a, b, c, |
522 |
U gTKE, |
U gTKE, |
523 |
I myThid ) |
I myThid ) |
524 |
C |
|
525 |
C now update TKE |
C now update TKE |
|
C |
|
526 |
DO K=1,Nr |
DO K=1,Nr |
527 |
DO J=jMin,jMax |
DO J=jMin,jMax |
528 |
DO I=iMin,iMax |
DO I=iMin,iMax |