--- MITgcm/pkg/ggl90/ggl90_calc.F 2009/01/30 02:23:56 1.11 +++ MITgcm/pkg/ggl90/ggl90_calc.F 2009/11/14 14:27:56 1.15 @@ -1,4 +1,4 @@ -C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/ggl90/ggl90_calc.F,v 1.11 2009/01/30 02:23:56 dfer Exp $ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/ggl90/ggl90_calc.F,v 1.15 2009/11/14 14:27:56 dfer Exp $ C $Name: $ #include "GGL90_OPTIONS.h" @@ -11,29 +11,26 @@ I bi, bj, myTime, myThid ) C !DESCRIPTION: \bv -C /==========================================================\ +C *==========================================================* C | SUBROUTINE GGL90_CALC | C | o Compute all GGL90 fields defined in GGL90.h | -C |==========================================================| +C *==========================================================* C | Equation numbers refer to | C | Gaspar et al. (1990), JGR 95 (C9), pp 16,179 | C | Some parts of the implementation follow Blanke and | C | Delecuse (1993), JPO, and OPA code, in particular the | C | computation of the | C | mixing length = max(min(lk,depth),lkmin) | -C \==========================================================/ - IMPLICIT NONE -C -C-------------------------------------------------------------------- +C *==========================================================* C global parameters updated by ggl90_calc -C GGL90TKE - sub-grid turbulent kinetic energy (m^2/s^2) -C GGL90viscAz - GGL90 eddy viscosity coefficient (m^2/s) -C GGL90diffKzT - GGL90 diffusion coefficient for temperature (m^2/s) -C +C GGL90TKE :: sub-grid turbulent kinetic energy (m^2/s^2) +C GGL90viscAz :: GGL90 eddy viscosity coefficient (m^2/s) +C GGL90diffKzT :: GGL90 diffusion coefficient for temperature (m^2/s) C \ev C !USES: ============================================================ + IMPLICIT NONE #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" @@ -43,18 +40,19 @@ #include "GRID.h" C !INPUT PARAMETERS: =================================================== -c Routine arguments -c bi, bj - array indices on which to apply calculations -c myTime - Current time in simulation - +C Routine arguments +C bi, bj :: array indices on which to apply calculations +C myTime :: Current time in simulation +C myThid :: My Thread Id number INTEGER bi, bj - INTEGER myThid _RL myTime + INTEGER myThid +CEOP #ifdef ALLOW_GGL90 C !LOCAL VARIABLES: ==================================================== -c Local constants +C Local constants C iMin, iMax, jMin, jMax, I, J - array computation indices C K, Kp1, km1, kSurf, kBottom - vertical loop indices C ab15, ab05 - weights for implicit timestepping @@ -86,22 +84,24 @@ _RL RiNumber _RL TKEdissipation _RL tempU, tempV, prTemp - _RL MaxLength + _RL MaxLength, tmpmlx _RL TKEPrandtlNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL GGL90mixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) - _RL rMixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) + _RL rMixingLength (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) + _RL mxLength_Dn (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL KappaE (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL rhoK (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL rhoKm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL totalDepth (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL gTKE (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) -C tri-diagonal matrix +C- tri-diagonal matrix _RL a(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL c(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) + INTEGER errCode #ifdef ALLOW_GGL90_HORIZDIFF -C xA, yA - area of lateral faces -C dfx, dfy - diffusive flux across lateral faces +C- xA, yA - area of lateral faces +C- dfx, dfy - diffusive flux across lateral faces _RL xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy) @@ -113,7 +113,6 @@ p8=0.125 _d 0 p16=0.0625 _d 0 #endif -CEOP iMin = 2-OLx iMax = sNx+OLx-1 jMin = 2-OLy @@ -121,7 +120,7 @@ C set separate time step (should be deltaTtracer) deltaTggl90 = dTtracerLev(1) -C + kSurf = 1 C implicit timestepping weights for dissipation ab15 = 1.5 _d 0 @@ -135,17 +134,19 @@ DO I=1-Olx,sNx+Olx gTKE(I,J,K) = 0. _d 0 KappaE(I,J,K) = 0. _d 0 - TKEPrandtlNumber(I,J,K) = 0. _d 0 + TKEPrandtlNumber(I,J,K) = 1. _d 0 GGL90mixingLength(I,J,K) = GGL90mixingLengthMin - rMixingLength(I,J,K) = 0. _d 0 ENDDO ENDDO ENDDO DO J=1-Oly,sNy+Oly DO I=1-Olx,sNx+Olx - rhoK (I,J) = 0. _d 0 - rhoKm1 (I,J) = 0. _d 0 - totalDepth(I,J) = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj) + rhoK(I,J) = 0. _d 0 + rhoKm1(I,J) = 0. _d 0 + totalDepth(I,J) = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj) + rMixingLength(i,j,1) = 0. _d 0 + mxLength_Dn(I,J,1) = GGL90mixingLengthMin + SQRTTKE(i,j,1) = SQRT( GGL90TKE(i,j,1,bi,bj) ) ENDDO ENDDO @@ -167,9 +168,8 @@ DO J=jMin,jMax DO I=iMin,iMax SQRTTKE(i,j,k)=SQRT( GGL90TKE(I,J,K,bi,bj) ) -C + C buoyancy frequency -C Nsquare(i,j,k) = - gravity*recip_rhoConst*recip_drC(K) & * ( rhoKm1(I,J) - rhoK(I,J) )*maskC(I,J,K,bi,bj) cC vertical shear term (dU/dz)^2+(dV/dz)^2 @@ -192,29 +192,54 @@ ENDDO ENDDO -C- Impose upper bound for mixing length (total depth) +C- Impose upper and lower bound for mixing length IF ( mxlMaxFlag .EQ. 0 ) THEN +C- DO k=2,Nr DO J=jMin,jMax DO I=iMin,iMax MaxLength=totalDepth(I,J) GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K), - & MaxLength) + & MaxLength) ENDDO ENDDO ENDDO + + DO k=2,Nr + DO J=jMin,jMax + DO I=iMin,iMax + GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K), + & GGL90mixingLengthMin) + rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K) + ENDDO + ENDDO + ENDDO + ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN +C- DO k=2,Nr DO J=jMin,jMax DO I=iMin,iMax MaxLength=MIN(Ro_surf(I,J,bi,bj)-rF(k),rF(k)-R_low(I,J,bi,bj)) c MaxLength=MAX(MaxLength,20. _d 0) GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K), - & MaxLength) + & MaxLength) + ENDDO + ENDDO + ENDDO + + DO k=2,Nr + DO J=jMin,jMax + DO I=iMin,iMax + GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K), + & GGL90mixingLengthMin) + rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K) ENDDO ENDDO ENDDO + ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN +C- DO k=2,Nr DO J=jMin,jMax DO I=iMin,iMax @@ -234,23 +259,71 @@ DO I=iMin,iMax GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K), & GGL90mixingLength(I,J,K+1)+drF(k)) - ENDDO - ENDDO + ENDDO + ENDDO ENDDO - ELSE - STOP 'GGL90_CALC: Wrong mxlMaxFlag (mixing lenght limit)' - ENDIF -C- Impose minimum mixing length (to avoid division by zero) - DO k=2,Nr + DO k=2,Nr + DO J=jMin,jMax + DO I=iMin,iMax + GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K), + & GGL90mixingLengthMin) + rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K) + ENDDO + ENDDO + ENDDO + + ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN +C- + DO k=2,Nr + DO J=jMin,jMax + DO I=iMin,iMax + mxLength_Dn(I,J,K) = MIN(GGL90mixingLength(I,J,K), + & mxLength_Dn(I,J,K-1)+drF(k-1)) + ENDDO + ENDDO + ENDDO DO J=jMin,jMax DO I=iMin,iMax - GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K), - & GGL90mixingLengthMin) - rMixingLength(I,J,K) = 1. _d 0 /GGL90mixingLength(I,J,K) + GGL90mixingLength(I,J,Nr) = MIN(GGL90mixingLength(I,J,Nr), + & GGL90mixingLengthMin+drF(Nr)) + ENDDO + ENDDO + DO k=Nr-1,2,-1 + DO J=jMin,jMax + DO I=iMin,iMax + GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K), + & GGL90mixingLength(I,J,K+1)+drF(k)) + ENDDO + ENDDO + ENDDO + + DO k=2,Nr + DO J=jMin,jMax + DO I=iMin,iMax + GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K), + & mxLength_Dn(I,J,K)) + tmpmlx = SQRT( GGL90mixingLength(I,J,K)*mxLength_Dn(I,J,K) ) + tmpmlx = MAX( tmpmlx, GGL90mixingLengthMin) + rMixingLength(I,J,K) = 1. _d 0 / tmpmlx + ENDDO ENDDO ENDDO - ENDDO + + ELSE + STOP 'GGL90_CALC: Wrong mxlMaxFlag (mixing length limit)' + ENDIF + +C- Impose minimum mixing length (to avoid division by zero) +c DO k=2,Nr +c DO J=jMin,jMax +c DO I=iMin,iMax +c GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K), +c & GGL90mixingLengthMin) +c rMixingLength(I,J,K) = 1. _d 0 /GGL90mixingLength(I,J,K) +c ENDDO +c ENDDO +c ENDDO DO k=2,Nr Km1 = K-1 @@ -264,18 +337,19 @@ & -( vVel(I,J,K ,bi,bj)+vVel(I,J+1,K ,bi,bj)) ) & *recip_drC(K) verticalShear = tempU*tempU + tempV*tempV - RiNumber = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps) + RiNumber = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps) C compute Prandtl number (always greater than 0) prTemp = 1. _d 0 IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber TKEPrandtlNumber(I,J,K) = MIN(10. _d 0,prTemp) +c TKEPrandtlNumber(I,J,K) = 1. _d 0 C viscosity and diffusivity KappaM = GGL90ck*GGL90mixingLength(I,J,K)*SQRTTKE(i,j,k) KappaH = KappaM/TKEPrandtlNumber(I,J,K) C Set a minium (= background) and maximum value - KappaM = MAX(KappaM,viscAr) + KappaM = MAX(KappaM,viscArNr(k)) KappaH = MAX(KappaH,diffKrNrT(k)) KappaM = MIN(KappaM,GGL90viscMax) KappaH = MIN(KappaH,GGL90diffMax) @@ -346,7 +420,7 @@ & -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj) & *( (dfx(i+1,j)-dfx(i,j)) & +(dfy(i,j+1)-dfy(i,j)) - & ) + & )*deltaTggl90 ENDDO ENDDO C end of k-loop @@ -367,14 +441,16 @@ ENDDO ENDDO DO k=2,Nr - km1=max(2,k-1) + km1=MAX(2,k-1) DO j=jMin,jMax DO i=iMin,iMax +C- We keep recip_hFacC in the diffusive flux calculation, +C- but no hFacC in TKE volume control +C- No need for maskC(k-1) with recip_hFacC(k-1) a(i,j,k) = -deltaTggl90 -c & *recip_drF(km1)*recip_hFacI(i,j,k,bi,bj) & *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj) & *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1)) - & *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) ENDDO ENDDO ENDDO @@ -387,30 +463,33 @@ DO k=2,Nr DO j=jMin,jMax DO i=iMin,iMax - kp1=min(klowC(i,j,bi,bj),k+1) + kp1=MAX(1,MIN(klowC(i,j,bi,bj),k+1)) +C- We keep recip_hFacC in the diffusive flux calculation, +C- but no hFacC in TKE volume control +C- No need for maskC(k) with recip_hFacC(k) c(i,j,k) = -deltaTggl90 -c & *recip_drF( k )*recip_hFacI(i,j,k,bi,bj) & *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj) & *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1)) - & *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) ENDDO ENDDO ENDDO C-- Center diagonal DO k=1,Nr + km1 = MAX(k-1,1) DO j=jMin,jMax DO i=iMin,iMax b(i,j,k) = 1. _d 0 - c(i,j,k) - a(i,j,k) - & + ab15*deltaTggl90*GGL90ceps*SQRT(GGL90TKE(I,J,K,bi,bj)) - & *rMixingLength(I,J,K)*maskC(i,j,k,bi,bj) + & + ab15*deltaTggl90*GGL90ceps*SQRTTKE(I,J,K) + & * rMixingLength(I,J,K) + & * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj) ENDDO ENDDO ENDDO C end set up matrix -C C Apply boundary condition -C + kp1 = MIN(Nr,kSurf+1) DO J=jMin,jMax DO I=iMin,iMax C estimate friction velocity uStar from surface forcing @@ -423,6 +502,9 @@ C Dirichlet surface boundary condition for TKE gTKE(I,J,kSurf) = MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare) & *maskC(I,J,kSurf,bi,bj) + gTKE(i,j,kp1) = gTKE(i,j,kp1) + & - a(i,j,kp1)*gTKE(i,j,kSurf) + a(i,j,kp1) = 0. _d 0 C Dirichlet bottom boundary condition for TKE = GGL90TKEbottom kBottom = MAX(kLowC(I,J,bi,bj),1) gTKE(I,J,kBottom) = gTKE(I,J,kBottom) @@ -430,16 +512,19 @@ c(I,J,kBottom) = 0. _d 0 ENDDO ENDDO -C + C solve tri-diagonal system, and store solution on gTKE (previously rhs) -C - CALL GGL90_SOLVE( bi, bj, iMin, iMax, jMin, jMax, - I a, b, c, - U gTKE, - I myThid ) -C + CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax, + I a, b, c, + U gTKE, + O errCode, + I bi, bj, myThid ) +c CALL GGL90_SOLVE( bi, bj, iMin, iMax, jMin, jMax, +c I a, b, c, +c U gTKE, +c I myThid ) + C now update TKE -C DO K=1,Nr DO J=jMin,jMax DO I=iMin,iMax @@ -505,4 +590,3 @@ RETURN END -