/[MITgcm]/MITgcm/pkg/ggl90/ggl90_calc.F
ViewVC logotype

Diff of /MITgcm/pkg/ggl90/ggl90_calc.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.19 by jmc, Thu Aug 19 23:52:37 2010 UTC revision 1.34 by jmc, Thu Jan 14 17:50:25 2016 UTC
# Line 8  C !ROUTINE: GGL90_CALC Line 8  C !ROUTINE: GGL90_CALC
8    
9  C !INTERFACE: ======================================================  C !INTERFACE: ======================================================
10        SUBROUTINE GGL90_CALC(        SUBROUTINE GGL90_CALC(
11       I     bi, bj, myTime, myIter, myThid )       I                 bi, bj, sigmaR, myTime, myIter, myThid )
12    
13  C !DESCRIPTION: \bv  C !DESCRIPTION: \bv
14  C     *==========================================================*  C     *==========================================================*
# Line 35  C !USES: =============================== Line 35  C !USES: ===============================
35  #include "EEPARAMS.h"  #include "EEPARAMS.h"
36  #include "PARAMS.h"  #include "PARAMS.h"
37  #include "DYNVARS.h"  #include "DYNVARS.h"
 #include "GGL90.h"  
38  #include "FFIELDS.h"  #include "FFIELDS.h"
39  #include "GRID.h"  #include "GRID.h"
40    #include "GGL90.h"
41    
42  C !INPUT PARAMETERS: ===================================================  C !INPUT PARAMETERS: ===================================================
43  C Routine arguments  C Routine arguments
44  C     bi, bj :: array indices on which to apply calculations  C     bi, bj :: Current tile indices
45    C     sigmaR :: Vertical gradient of iso-neutral density
46  C     myTime :: Current time in simulation  C     myTime :: Current time in simulation
47  C     myIter :: Current time-step number  C     myIter :: Current time-step number
48  C     myThid :: My Thread Id number  C     myThid :: My Thread Id number
49        INTEGER bi, bj        INTEGER bi, bj
50          _RL     sigmaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
51        _RL     myTime        _RL     myTime
52        INTEGER myIter        INTEGER myIter
53        INTEGER myThid        INTEGER myThid
 CEOP  
54    
55  #ifdef ALLOW_GGL90  #ifdef ALLOW_GGL90
56    
# Line 58  C Local constants Line 59  C Local constants
59  C     iMin,iMax,jMin,jMax :: index boundaries of computation domain  C     iMin,iMax,jMin,jMax :: index boundaries of computation domain
60  C     i, j, k, kp1,km1 :: array computation indices  C     i, j, k, kp1,km1 :: array computation indices
61  C     kSurf, kBottom   :: vertical indices of domain boundaries  C     kSurf, kBottom   :: vertical indices of domain boundaries
62    C     hFac/hFacI       :: fractional thickness of W-cell
63  C     explDissFac      :: explicit Dissipation Factor (in [0-1])  C     explDissFac      :: explicit Dissipation Factor (in [0-1])
64  C     implDissFac      :: implicit Dissipation Factor (in [0-1])  C     implDissFac      :: implicit Dissipation Factor (in [0-1])
65  C     uStarSquare      :: square of friction velocity  C     uStarSquare      :: square of friction velocity
# Line 72  C     GGL90mixingLength:: mixing length Line 74  C     GGL90mixingLength:: mixing length
74  C         rMixingLength:: inverse of mixing length  C         rMixingLength:: inverse of mixing length
75  C     totalDepth       :: thickness of water column (inverse of recip_Rcol)  C     totalDepth       :: thickness of water column (inverse of recip_Rcol)
76  C     TKEPrandtlNumber :: here, an empirical function of the Richardson number  C     TKEPrandtlNumber :: here, an empirical function of the Richardson number
 C     rhoK, rhoKm1     :: density at layer k and km1 (relative to k)  
77        INTEGER iMin ,iMax ,jMin ,jMax        INTEGER iMin ,iMax ,jMin ,jMax
78        INTEGER i, j, k, kp1, km1, kSurf, kBottom        INTEGER i, j, k, kp1, km1, kSurf, kBottom
79        _RL     explDissFac, implDissFac        _RL     explDissFac, implDissFac
80        _RL     uStarSquare        _RL     uStarSquare
81        _RL     verticalShear        _RL     verticalShear(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
82        _RL     KappaM, KappaH        _RL     KappaM(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
83          _RL     KappaH
84  c     _RL     Nsquare  c     _RL     Nsquare
85        _RL     Nsquare(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     Nsquare(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
86        _RL     deltaTggl90        _RL     deltaTggl90
87  c     _RL     SQRTTKE  c     _RL     SQRTTKE
88        _RL     SQRTTKE(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     SQRTTKE(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
89        _RL     RiNumber        _RL     RiNumber
90    #ifdef ALLOW_GGL90_IDEMIX
91          _RL     IDEMIX_RiNumber
92    #endif
93        _RL     TKEdissipation        _RL     TKEdissipation
94        _RL     tempU, tempV, prTemp        _RL     tempU, tempUp, tempV, tempVp, prTemp
95        _RL     MaxLength, tmpmlx, tmpVisc        _RL     MaxLength, tmpmlx, tmpVisc
96        _RL     TKEPrandtlNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     TKEPrandtlNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
97        _RL     GGL90mixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     GGL90mixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
98        _RL     rMixingLength    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     rMixingLength    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
99        _RL     mxLength_Dn      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     mxLength_Dn      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
100        _RL     KappaE           (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)  
101        _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102        _RL     GGL90visctmp     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     GGL90visctmp     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
103    #ifdef ALLOW_DIAGNOSTICS
104          _RL     surf_flx_tke     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105    #endif /* ALLOW_DIAGNOSTICS */
106    C     hFac(I)  :: fractional thickness of W-cell
107          _RL       hFac
108    #ifdef ALLOW_GGL90_IDEMIX
109          _RL       hFacI(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
110    #endif /* ALLOW_GGL90_IDEMIX */
111          _RL recip_hFacI(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
112  C-    tri-diagonal matrix  C-    tri-diagonal matrix
113        _RL     a(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     a3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
114        _RL     b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     b3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
115        _RL     c(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     c3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
116        INTEGER errCode        INTEGER errCode
117  #ifdef ALLOW_GGL90_HORIZDIFF  #ifdef ALLOW_GGL90_HORIZDIFF
118  C     xA, yA   :: area of lateral faces  C     xA, yA   :: area of lateral faces
# Line 114  C     gTKE     :: right hand side of dif Line 126  C     gTKE     :: right hand side of dif
126  #endif /* ALLOW_GGL90_HORIZDIFF */  #endif /* ALLOW_GGL90_HORIZDIFF */
127  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
128        _RL p4, p8, p16        _RL p4, p8, p16
       p4=0.25 _d 0  
       p8=0.125 _d 0  
       p16=0.0625 _d 0  
129  #endif  #endif
130        iMin = 2-OLx  CEOP
131        iMax = sNx+OLx-1  
132        jMin = 2-OLy        PARAMETER( iMin = 2-OLx, iMax = sNx+OLx-1 )
133        jMax = sNy+OLy-1        PARAMETER( jMin = 2-OLy, jMax = sNy+OLy-1 )
134    #ifdef ALLOW_GGL90_SMOOTH
135          p4  = 0.25   _d 0
136          p8  = 0.125  _d 0
137          p16 = 0.0625 _d 0
138    #endif
139    
140  C     set separate time step (should be deltaTtracer)  C     set separate time step (should be deltaTtracer)
141        deltaTggl90 = dTtracerLev(1)        deltaTggl90 = dTtracerLev(1)
# Line 131  C     explicit/implicit timestepping wei Line 145  C     explicit/implicit timestepping wei
145        explDissFac = 0. _d 0        explDissFac = 0. _d 0
146        implDissFac = 1. _d 0 - explDissFac        implDissFac = 1. _d 0 - explDissFac
147    
148    C     For nonlinear free surface and especially with r*-coordinates, the
149    C     hFacs change every timestep, so we need to update them here in the
150    C     case of using IDEMIX.
151           DO K=1,Nr
152            Km1 = MAX(K-1,1)
153            DO j=1-OLy,sNy+OLy
154             DO i=1-OLx,sNx+OLx
155              hFac =
156         &         MIN(.5 _d 0,_hFacC(i,j,km1,bi,bj) ) +
157         &         MIN(.5 _d 0,_hFacC(i,j,k  ,bi,bj) )
158              recip_hFacI(I,J,K)=0. _d 0
159              IF ( hFac .NE. 0. _d 0 )
160         &         recip_hFacI(I,J,K)=1. _d 0/hFac
161    #ifdef ALLOW_GGL90_IDEMIX
162              hFacI(i,j,k) = hFac
163    #endif /* ALLOW_GGL90_IDEMIX */
164             ENDDO
165            ENDDO
166           ENDDO
167    
168  C     Initialize local fields  C     Initialize local fields
169        DO k = 1, Nr        DO k = 1, Nr
170         DO j=1-Oly,sNy+Oly         DO j=1-OLy,sNy+OLy
171          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
172             rMixingLength(i,j,k)     = 0. _d 0
173             mxLength_Dn(i,j,k)       = 0. _d 0
174             GGL90visctmp(i,j,k)      = 0. _d 0
175           KappaE(i,j,k)            = 0. _d 0           KappaE(i,j,k)            = 0. _d 0
176           TKEPrandtlNumber(i,j,k)  = 1. _d 0           TKEPrandtlNumber(i,j,k)  = 1. _d 0
177           GGL90mixingLength(i,j,k) = GGL90mixingLengthMin           GGL90mixingLength(i,j,k) = GGL90mixingLengthMin
178           GGL90visctmp(i,j,k)      = 0. _d 0           GGL90visctmp(i,j,k)      = 0. _d 0
179    #ifndef SOLVE_DIAGONAL_LOWMEMORY
180             a3d(i,j,k) = 0. _d 0
181             b3d(i,j,k) = 1. _d 0
182             c3d(i,j,k) = 0. _d 0
183    #endif
184             Nsquare(i,j,k) = 0. _d 0
185             SQRTTKE(i,j,k) = 0. _d 0
186          ENDDO          ENDDO
187         ENDDO         ENDDO
188        ENDDO        ENDDO
189        DO j=1-Oly,sNy+Oly        DO j=1-OLy,sNy+OLy
190         DO i=1-Olx,sNx+Olx         DO i=1-OLx,sNx+OLx
191          rhoK(i,j)          = 0. _d 0          KappaM(i,j)        = 0. _d 0
192          rhoKm1(i,j)        = 0. _d 0          verticalShear(i,j) = 0. _d 0
193          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)
194          rMixingLength(i,j,1) = 0. _d 0          rMixingLength(i,j,1) = 0. _d 0
195          mxLength_Dn(i,j,1) = GGL90mixingLengthMin          mxLength_Dn(i,j,1) = GGL90mixingLengthMin
196          SQRTTKE(i,j,1) = SQRT( GGL90TKE(i,j,1,bi,bj) )          SQRTTKE(i,j,1) = SQRT( GGL90TKE(i,j,1,bi,bj) )
197    #ifdef ALLOW_GGL90_HORIZDIFF
198            xA(i,j)  = 0. _d 0
199            yA(i,j)  = 0. _d 0
200            dfx(i,j) = 0. _d 0
201            dfy(i,j) = 0. _d 0
202            gTKE(i,j) = 0. _d 0
203    #endif /* ALLOW_GGL90_HORIZDIFF */
204         ENDDO         ENDDO
205        ENDDO        ENDDO
206    
207  C     start k-loop  #ifdef ALLOW_GGL90_IDEMIX
208          IF ( useIDEMIX) CALL GGL90_IDEMIX(
209         &   bi, bj, hFacI, recip_hFacI, sigmaR, myTime, myIter, myThid )
210    #endif /* ALLOW_GGL90_IDEMIX */
211    
212        DO k = 2, Nr        DO k = 2, Nr
        km1 = k-1  
 c      kp1 = MIN(Nr,k+1)  
        CALL FIND_RHO_2D(  
      I      iMin, iMax, jMin, jMax, k,  
      I      theta(1-OLx,1-OLy,km1,bi,bj), salt(1-OLx,1-OLy,km1,bi,bj),  
      O      rhoKm1,  
      I      km1, bi, bj, myThid )  
   
        CALL FIND_RHO_2D(  
      I      iMin, iMax, jMin, jMax, k,  
      I      theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),  
      O      rhoK,  
      I      k, bi, bj, myThid )  
213         DO j=jMin,jMax         DO j=jMin,jMax
214          DO i=iMin,iMax          DO i=iMin,iMax
215           SQRTTKE(i,j,k)=SQRT( GGL90TKE(i,j,k,bi,bj) )           SQRTTKE(i,j,k)=SQRT( GGL90TKE(i,j,k,bi,bj) )
216    
217  C     buoyancy frequency  C     buoyancy frequency
218           Nsquare(i,j,k) = - gravity*recip_rhoConst*recip_drC(k)           Nsquare(i,j,k) = gravity*gravitySign*recip_rhoConst
219       &        * ( rhoKm1(i,j) - rhoK(i,j) )*maskC(i,j,k,bi,bj)       &                  * sigmaR(i,j,k)
220  cC     vertical shear term (dU/dz)^2+(dV/dz)^2  C     vertical shear term (dU/dz)^2+(dV/dz)^2 is computed later
221  c         tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)  C     to save some memory
 c     &                 -( uVel(i,j,k  ,bi,bj)+uVel(i+1,j,k  ,bi,bj)) )  
 c     &        *recip_drC(k)  
 c         tempV= .5 _d 0*( vVel(i,j,km1,bi,bj)+vVel(i,j+1,km1,bi,bj)  
 c     &                 -( vVel(i,j,k  ,bi,bj)+vVel(i,j+1,k  ,bi,bj)) )  
 c     &        *recip_drC(k)  
 c         verticalShear = tempU*tempU + tempV*tempV  
 c         RiNumber   = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)  
 cC     compute Prandtl number (always greater than 0)  
 c         prTemp = 1. _d 0  
 c         IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber  
 c         TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)  
222  C     mixing length  C     mixing length
223           GGL90mixingLength(i,j,k) = SQRTTWO *           GGL90mixingLength(i,j,k) = SQRTTWO *
224       &        SQRTTKE(i,j,k)/SQRT( MAX(Nsquare(i,j,k),GGL90eps) )       &        SQRTTKE(i,j,k)/SQRT( MAX(Nsquare(i,j,k),GGL90eps) )
# Line 195  C     mixing length Line 226  C     mixing length
226         ENDDO         ENDDO
227        ENDDO        ENDDO
228    
229  C- Impose upper and lower bound for mixing length  C- ensure mixing between first and second level
230          IF (mxlSurfFlag) THEN
231           DO j=jMin,jMax
232            DO i=iMin,iMax
233             GGL90mixingLength(i,j,2)=drF(1)
234            ENDDO
235           ENDDO
236          ENDIF
237    
238    C--   Impose upper and lower bound for mixing length
239    C--   Impose minimum mixing length to avoid division by zero
240        IF ( mxlMaxFlag .EQ. 0 ) THEN        IF ( mxlMaxFlag .EQ. 0 ) THEN
241  C-  
242         DO k=2,Nr         DO k=2,Nr
243          DO j=jMin,jMax          DO j=jMin,jMax
244           DO i=iMin,iMax           DO i=iMin,iMax
# Line 219  C- Line 260  C-
260         ENDDO         ENDDO
261    
262        ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN
263  C-  
264         DO k=2,Nr         DO k=2,Nr
265          DO j=jMin,jMax          DO j=jMin,jMax
266           DO i=iMin,iMax           DO i=iMin,iMax
# Line 242  c         MaxLength=MAX(MaxLength,20. _d Line 283  c         MaxLength=MAX(MaxLength,20. _d
283         ENDDO         ENDDO
284    
285        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN
286  C-  
 cgf ensure mixing between first and second level  
 c      DO j=jMin,jMax  
 c        DO i=iMin,iMax  
 c         GGL90mixingLength(i,j,2)=drF(1)  
 c        ENDDO  
 c      ENDDO  
 cgf  
287         DO k=2,Nr         DO k=2,Nr
288          DO j=jMin,jMax          DO j=jMin,jMax
289           DO i=iMin,iMax           DO i=iMin,iMax
# Line 284  cgf Line 318  cgf
318         ENDDO         ENDDO
319    
320        ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN
321  C-  
322         DO k=2,Nr         DO k=2,Nr
323          DO j=jMin,jMax          DO j=jMin,jMax
324           DO i=iMin,iMax           DO i=iMin,iMax
# Line 324  C- Line 358  C-
358         STOP 'GGL90_CALC: Wrong mxlMaxFlag (mixing length limit)'         STOP 'GGL90_CALC: Wrong mxlMaxFlag (mixing length limit)'
359        ENDIF        ENDIF
360    
361  C- Impose minimum mixing length (to avoid division by zero)  C     start "proper" k-loop (the code above was moved out and up to
362  c      DO k=2,Nr  C     implemement various mixing parameters efficiently)
 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  
   
363        DO k=2,Nr        DO k=2,Nr
364         km1 = k-1         km1 = k-1
365    
366  #ifdef ALLOW_GGL90_HORIZDIFF  #ifdef ALLOW_GGL90_HORIZDIFF
367         IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN        IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
368  C     horizontal diffusion of TKE (requires an exchange in  C     horizontal diffusion of TKE (requires an exchange in
369  C      do_fields_blocking_exchanges)  C      do_fields_blocking_exchanges)
370  C     common factors  C     common factors
371          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
372           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
373            xA(i,j) = _dyG(i,j,bi,bj)            xA(i,j) = _dyG(i,j,bi,bj)*drC(k)*
374       &         *drF(k)*_hFacW(i,j,k,bi,bj)       &                 (min(.5 _d 0,_hFacW(i,j,k-1,bi,bj) ) +
375            yA(i,j) = _dxG(i,j,bi,bj)       &                  min(.5 _d 0,_hFacW(i,j,k  ,bi,bj) ) )
376       &         *drF(k)*_hFacS(i,j,k,bi,bj)            yA(i,j) = _dxG(i,j,bi,bj)*drC(k)*
377         &                 (min(.5 _d 0,_hFacS(i,j,k-1,bi,bj) ) +
378         &                  min(.5 _d 0,_hFacS(i,j,k  ,bi,bj) ) )
379           ENDDO           ENDDO
380          ENDDO          ENDDO
381  C     Compute diffusive fluxes  C     Compute diffusive fluxes
382  C     ... across x-faces  C     ... across x-faces
383          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
384           dfx(1-Olx,j)=0. _d 0           dfx(1-OLx,j)=0. _d 0
385           DO i=1-Olx+1,sNx+Olx           DO i=1-OLx+1,sNx+OLx
386            dfx(i,j) = -GGL90diffTKEh*xA(i,j)            dfx(i,j) = -GGL90diffTKEh*xA(i,j)
387       &      *_recip_dxC(i,j,bi,bj)       &      *_recip_dxC(i,j,bi,bj)
388       &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i-1,j,k,bi,bj))       &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i-1,j,k,bi,bj))
389    #ifdef ISOTROPIC_COS_SCALING
390       &      *CosFacU(j,bi,bj)       &      *CosFacU(j,bi,bj)
391    #endif /* ISOTROPIC_COS_SCALING */
392           ENDDO           ENDDO
393          ENDDO          ENDDO
394  C     ... across y-faces  C     ... across y-faces
395          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
396           dfy(i,1-Oly)=0. _d 0           dfy(i,1-OLy)=0. _d 0
397          ENDDO          ENDDO
398          DO j=1-Oly+1,sNy+Oly          DO j=1-OLy+1,sNy+OLy
399           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
400            dfy(i,j) = -GGL90diffTKEh*yA(i,j)            dfy(i,j) = -GGL90diffTKEh*yA(i,j)
401       &      *_recip_dyC(i,j,bi,bj)       &      *_recip_dyC(i,j,bi,bj)
402       &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i,j-1,k,bi,bj))       &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i,j-1,k,bi,bj))
# Line 377  C     ... across y-faces Line 406  C     ... across y-faces
406           ENDDO           ENDDO
407          ENDDO          ENDDO
408  C     Compute divergence of fluxes  C     Compute divergence of fluxes
409          DO j=1-Oly,sNy+Oly-1          DO j=1-OLy,sNy+OLy-1
410           DO i=1-Olx,sNx+Olx-1           DO i=1-OLx,sNx+OLx-1
411            gTKE(i,j) =            gTKE(i,j) = -recip_drC(k)*recip_rA(i,j,bi,bj)
412       &    -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj)       &         *recip_hFacI(i,j,k)
413       &         *( (dfx(i+1,j)-dfx(i,j))       &         *((dfx(i+1,j)-dfx(i,j))
414       &           +(dfy(i,j+1)-dfy(i,j))       &         + (dfy(i,j+1)-dfy(i,j)) )
      &          )  
415           ENDDO           ENDDO
416          ENDDO          ENDDO
417  C      end if GGL90diffTKEh .eq. 0.  C     end if GGL90diffTKEh .eq. 0.
418         ENDIF         ENDIF
419  #endif /* ALLOW_GGL90_HORIZDIFF */  #endif /* ALLOW_GGL90_HORIZDIFF */
420    
421    C     viscosity and diffusivity
422         DO j=jMin,jMax         DO j=jMin,jMax
423          DO i=iMin,iMax          DO i=iMin,iMax
424  C     vertical shear term (dU/dz)^2+(dV/dz)^2           KappaM(i,j) = GGL90ck*GGL90mixingLength(i,j,k)*SQRTTKE(i,j,k)
425           tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)           GGL90visctmp(i,j,k) = MAX(KappaM(i,j),diffKrNrS(k))
426       &                 -( uVel(i,j,k  ,bi,bj)+uVel(i+1,j,k  ,bi,bj)) )       &                            * maskC(i,j,k,bi,bj)
427       &        *recip_drC(k)  C        note: storing GGL90visctmp like this, and using it later to compute
428           tempV= .5 _d 0*( vVel(i,j,km1,bi,bj)+vVel(i,j+1,km1,bi,bj)  C              GGL9rdiffKr etc. is robust in case of smoothing (e.g. see OPA)
429       &                 -( vVel(i,j,k  ,bi,bj)+vVel(i,j+1,k  ,bi,bj)) )           KappaM(i,j) = MAX(KappaM(i,j),viscArNr(k)) * maskC(i,j,k,bi,bj)
430       &        *recip_drC(k)          ENDDO
431           verticalShear = tempU*tempU + tempV*tempV         ENDDO
432           RiNumber = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)  
433    C     compute vertical shear (dU/dz)^2+(dV/dz)^2
434           IF ( calcMeanVertShear ) THEN
435    C     by averaging (@ grid-cell center) the 4 vertical shear compon @ U,V pos.
436            DO j=jMin,jMax
437             DO i=iMin,iMax
438              tempU  = ( uVel( i ,j,km1,bi,bj) - uVel( i ,j,k,bi,bj) )
439              tempUp = ( uVel(i+1,j,km1,bi,bj) - uVel(i+1,j,k,bi,bj) )
440              tempV  = ( vVel(i, j ,km1,bi,bj) - vVel(i, j ,k,bi,bj) )
441              tempVp = ( vVel(i,j+1,km1,bi,bj) - vVel(i,j+1,k,bi,bj) )
442              verticalShear(i,j) = (
443         &                 ( tempU*tempU + tempUp*tempUp )*halfRL
444         &               + ( tempV*tempV + tempVp*tempVp )*halfRL
445         &                         )*recip_drC(k)*recip_drC(k)
446             ENDDO
447            ENDDO
448           ELSE
449    C     from the averaged flow at grid-cell center (2 compon x 2 pos.)
450            DO j=jMin,jMax
451             DO i=iMin,iMax
452              tempU = ( ( uVel(i,j,km1,bi,bj) + uVel(i+1,j,km1,bi,bj) )
453         &             -( uVel(i,j,k  ,bi,bj) + uVel(i+1,j,k  ,bi,bj) )
454         &            )*halfRL*recip_drC(k)
455              tempV = ( ( vVel(i,j,km1,bi,bj) + vVel(i,j+1,km1,bi,bj) )
456         &             -( vVel(i,j,k  ,bi,bj) + vVel(i,j+1,k  ,bi,bj) )
457         &            )*halfRL*recip_drC(k)
458              verticalShear(i,j) = tempU*tempU + tempV*tempV
459             ENDDO
460            ENDDO
461           ENDIF
462    
463  C     compute Prandtl number (always greater than 0)  C     compute Prandtl number (always greater than 0)
464           prTemp = 1. _d 0  #ifdef ALLOW_GGL90_IDEMIX
465           IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber         IF ( useIDEMIX ) THEN
466           TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)          DO j=jMin,jMax
467  c         TKEPrandtlNumber(i,j,k) = 1. _d 0           DO i=iMin,iMax
468    C     account for partical cell factor in vertical shear:
469              verticalShear(i,j) = verticalShear(i,j)
470         &                       * recip_hFacI(i,j,k)*recip_hFacI(i,j,k)
471              RiNumber = MAX(Nsquare(i,j,k),0. _d 0)
472         &         /(verticalShear(i,j)+GGL90eps)
473    CML         IDEMIX_RiNumber = 1./GGL90eps
474              IDEMIX_RiNumber = MAX( KappaM(i,j)*Nsquare(i,j,k), 0. _d 0)/
475         &     (GGL90eps+IDEMIX_tau_d(i,j,k,bi,bj)*IDEMIX_E(i,j,k,bi,bj)**2)
476              prTemp         = MIN(5.*RiNumber, 6.6 _d 0*IDEMIX_RiNumber)
477              TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
478              TKEPrandtlNumber(i,j,k) = MAX( 1. _d 0,TKEPrandtlNumber(i,j,k))
479             ENDDO
480            ENDDO
481           ELSE
482    #else /* ndef ALLOW_GGL90_IDEMIX */
483           IF (.TRUE.) THEN
484    #endif /* ALLOW_GGL90_IDEMIX */
485            DO j=jMin,jMax
486             DO i=iMin,iMax
487              RiNumber = MAX(Nsquare(i,j,k),0. _d 0)
488         &         /(verticalShear(i,j)+GGL90eps)
489              prTemp = 1. _d 0
490              IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
491              TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
492             ENDDO
493            ENDDO
494           ENDIF
495    
496  C     viscosity and diffusivity         DO j=jMin,jMax
497           KappaM = GGL90ck*GGL90mixingLength(i,j,k)*SQRTTKE(i,j,k)          DO i=iMin,iMax
498           GGL90visctmp(i,j,k) = MAX(KappaM,diffKrNrT(k))  C        diffusivity
499       &                            * maskC(i,j,k,bi,bj)           KappaH = KappaM(i,j)/TKEPrandtlNumber(i,j,k)
500  c        note: storing GGL90visctmp like this, and using it later to compute           KappaE(i,j,k) = GGL90alpha * KappaM(i,j) * maskC(i,j,k,bi,bj)
 c              GGL9rdiffKr etc. is robust in case of smoothing (e.g. see OPA)  
          KappaM = MAX(KappaM,viscArNr(k)) * maskC(i,j,k,bi,bj)  
          KappaH = KappaM/TKEPrandtlNumber(i,j,k)  
          KappaE(i,j,k) = GGL90alpha * KappaM * maskC(i,j,k,bi,bj)  
501    
502  C     dissipation term  C     dissipation term
503           TKEdissipation = explDissFac*GGL90ceps           TKEdissipation = explDissFac*GGL90ceps
# Line 424  C     dissipation term Line 506  C     dissipation term
506  C     partial update with sum of explicit contributions  C     partial update with sum of explicit contributions
507           GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)           GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
508       &        + deltaTggl90*(       &        + deltaTggl90*(
509       &        + KappaM*verticalShear       &        + KappaM(i,j)*verticalShear(i,j)
510       &        - KappaH*Nsquare(i,j,k)       &        - KappaH*Nsquare(i,j,k)
511       &        - TKEdissipation       &        - TKEdissipation
512       &        )       &        )
513          ENDDO          ENDDO
514         ENDDO         ENDDO
515    
516    #ifdef ALLOW_GGL90_IDEMIX
517           IF ( useIDEMIX ) THEN
518    C     add IDEMIX contribution to the turbulent kinetic energy
519            DO j=jMin,jMax
520             DO i=iMin,iMax
521              GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
522         &         + deltaTggl90*(
523         &         + IDEMIX_tau_d(i,j,k,bi,bj)*IDEMIX_E(i,j,k,bi,bj)**2
524         &         )
525             ENDDO
526            ENDDO
527           ENDIF
528    #endif /* ALLOW_GGL90_IDEMIX */
529    
530  #ifdef ALLOW_GGL90_HORIZDIFF  #ifdef ALLOW_GGL90_HORIZDIFF
531         IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN         IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
532  C--    Add horiz. diffusion tendency  C--    Add horiz. diffusion tendency
# Line 454  C     set up matrix Line 550  C     set up matrix
550  C--   Lower diagonal  C--   Lower diagonal
551        DO j=jMin,jMax        DO j=jMin,jMax
552         DO i=iMin,iMax         DO i=iMin,iMax
553           a(i,j,1) = 0. _d 0           a3d(i,j,1) = 0. _d 0
554         ENDDO         ENDDO
555        ENDDO        ENDDO
556        DO k=2,Nr        DO k=2,Nr
# Line 464  C--   Lower diagonal Line 560  C--   Lower diagonal
560  C-    We keep recip_hFacC in the diffusive flux calculation,  C-    We keep recip_hFacC in the diffusive flux calculation,
561  C-    but no hFacC in TKE volume control  C-    but no hFacC in TKE volume control
562  C-    No need for maskC(k-1) with recip_hFacC(k-1)  C-    No need for maskC(k-1) with recip_hFacC(k-1)
563           a(i,j,k) = -deltaTggl90           a3d(i,j,k) = -deltaTggl90
564       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)
565       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))
566       &        *recip_drC(k)*maskC(i,j,k,bi,bj)       &        *recip_drC(k)*maskC(i,j,k,bi,bj)
# Line 474  C-    No need for maskC(k-1) with recip_ Line 570  C-    No need for maskC(k-1) with recip_
570  C--   Upper diagonal  C--   Upper diagonal
571        DO j=jMin,jMax        DO j=jMin,jMax
572         DO i=iMin,iMax         DO i=iMin,iMax
573           c(i,j,1)  = 0. _d 0           c3d(i,j,1)  = 0. _d 0
574         ENDDO         ENDDO
575        ENDDO        ENDDO
576        DO k=2,Nr        DO k=2,Nr
577         DO j=jMin,jMax         DO j=jMin,jMax
578          DO i=iMin,iMax          DO i=iMin,iMax
579            kp1=MAX(1,MIN(klowC(i,j,bi,bj),k+1))           kp1=MAX(1,MIN(klowC(i,j,bi,bj),k+1))
580  C-    We keep recip_hFacC in the diffusive flux calculation,  C-    We keep recip_hFacC in the diffusive flux calculation,
581  C-    but no hFacC in TKE volume control  C-    but no hFacC in TKE volume control
582  C-    No need for maskC(k) with recip_hFacC(k)  C-    No need for maskC(k) with recip_hFacC(k)
583            c(i,j,k) = -deltaTggl90           c3d(i,j,k) = -deltaTggl90
584       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)
585       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))
586       &        *recip_drC(k)*maskC(i,j,k-1,bi,bj)       &        *recip_drC(k)*maskC(i,j,k-1,bi,bj)
587          ENDDO          ENDDO
588         ENDDO         ENDDO
589        ENDDO        ENDDO
590    
591    #ifdef ALLOW_GGL90_IDEMIX
592          IF ( useIDEMIX ) THEN
593           DO k=2,Nr
594            DO j=jMin,jMax
595             DO i=iMin,iMax
596              a3d(i,j,k) = a3d(i,j,k)*recip_hFacI(i,j,k)
597              c3d(i,j,k) = c3d(i,j,k)*recip_hFacI(i,j,k)
598             ENDDO
599            ENDDO
600           ENDDO
601          ENDIF
602    #endif /* ALLOW_GGL90_IDEMIX */
603    
604          IF (.NOT.GGL90_dirichlet) THEN
605    C      Neumann bottom boundary condition for TKE: no flux from bottom
606           DO j=jMin,jMax
607            DO i=iMin,iMax
608             kBottom   = MAX(kLowC(i,j,bi,bj),1)
609             c3d(i,j,kBottom) = 0. _d 0
610            ENDDO
611           ENDDO
612          ENDIF
613    
614  C--   Center diagonal  C--   Center diagonal
615        DO k=1,Nr        DO k=1,Nr
616         km1 = MAX(k-1,1)         km1 = MAX(k-1,1)
617         DO j=jMin,jMax         DO j=jMin,jMax
618          DO i=iMin,iMax          DO i=iMin,iMax
619            b(i,j,k) = 1. _d 0 - c(i,j,k) - a(i,j,k)           b3d(i,j,k) = 1. _d 0 - c3d(i,j,k) - a3d(i,j,k)
620       &        + implDissFac*deltaTggl90*GGL90ceps*SQRTTKE(i,j,k)       &        + implDissFac*deltaTggl90*GGL90ceps*SQRTTKE(i,j,k)
621       &        * rMixingLength(i,j,k)       &        * rMixingLength(i,j,k)
622       &        * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)       &        * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)
623           ENDDO          ENDDO
624         ENDDO         ENDDO
625        ENDDO        ENDDO
626  C     end set up matrix  C     end set up matrix
# Line 520  C     Dirichlet surface boundary conditi Line 640  C     Dirichlet surface boundary conditi
640          GGL90TKE(i,j,kSurf,bi,bj) = maskC(i,j,kSurf,bi,bj)          GGL90TKE(i,j,kSurf,bi,bj) = maskC(i,j,kSurf,bi,bj)
641       &           *MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)       &           *MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)
642          GGL90TKE(i,j,kp1,bi,bj) = GGL90TKE(i,j,kp1,bi,bj)          GGL90TKE(i,j,kp1,bi,bj) = GGL90TKE(i,j,kp1,bi,bj)
643       &               - a(i,j,kp1)*GGL90TKE(i,j,kSurf,bi,bj)       &               - a3d(i,j,kp1)*GGL90TKE(i,j,kSurf,bi,bj)
644          a(i,j,kp1) = 0. _d 0          a3d(i,j,kp1) = 0. _d 0
 C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom  
         kBottom   = MAX(kLowC(i,j,bi,bj),1)  
         GGL90TKE(i,j,kBottom,bi,bj) = GGL90TKE(i,j,kBottom,bi,bj)  
      &                              - GGL90TKEbottom*c(i,j,kBottom)  
         c(i,j,kBottom) = 0. _d 0  
645         ENDDO         ENDDO
646        ENDDO        ENDDO
647    
648          IF (GGL90_dirichlet) THEN
649    C      Dirichlet bottom boundary condition for TKE = GGL90TKEbottom
650           DO j=jMin,jMax
651            DO i=iMin,iMax
652             kBottom   = MAX(kLowC(i,j,bi,bj),1)
653             GGL90TKE(i,j,kBottom,bi,bj) = GGL90TKE(i,j,kBottom,bi,bj)
654         &                              - GGL90TKEbottom*c3d(i,j,kBottom)
655             c3d(i,j,kBottom) = 0. _d 0
656            ENDDO
657           ENDDO
658          ENDIF
659    
660  C     solve tri-diagonal system  C     solve tri-diagonal system
661        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
662       I                        a, b, c,       I                        a3d, b3d, c3d,
663       U                        GGL90TKE,       U                        GGL90TKE(1-OLx,1-OLy,1,bi,bj),
664       O                        errCode,       O                        errCode,
665       I                        bi, bj, myThid )       I                        bi, bj, myThid )
666    
# Line 554  C     =============================== Line 681  C     ===============================
681         DO j=1,sNy         DO j=1,sNy
682          DO i=1,sNx          DO i=1,sNx
683  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
684           tmpVisc=           tmpVisc = (
685       &  (       &     p4 *    GGL90visctmp(i  ,j  ,k)*mskCor(i  ,j  ,bi,bj)
686       &   p4 *  GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)       &    +p8 *( ( GGL90visctmp(i-1,j  ,k)*mskCor(i-1,j  ,bi,bj)
687       &  +p8 *( GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj)       &           + GGL90visctmp(i+1,j  ,k)*mskCor(i+1,j  ,bi,bj) )
688       &       + GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj)       &         + ( GGL90visctmp(i  ,j-1,k)*mskCor(i  ,j-1,bi,bj)
689       &       + GGL90visctmp(i+1,j  ,k) * mskCor(i+1,j  ,bi,bj)       &           + GGL90visctmp(i  ,j+1,k)*mskCor(i  ,j+1,bi,bj) ) )
690       &       + GGL90visctmp(i  ,j+1,k) * mskCor(i  ,j+1,bi,bj))       &    +p16*( ( GGL90visctmp(i+1,j+1,k)*mskCor(i+1,j+1,bi,bj)
691       &  +p16*( GGL90visctmp(i+1,j+1,k) * mskCor(i+1,j+1,bi,bj)       &           + GGL90visctmp(i-1,j-1,k)*mskCor(i-1,j-1,bi,bj) )
692       &       + GGL90visctmp(i+1,j-1,k) * mskCor(i+1,j-1,bi,bj)       &         + ( GGL90visctmp(i+1,j-1,k)*mskCor(i+1,j-1,bi,bj)
693       &       + GGL90visctmp(i-1,j+1,k) * mskCor(i-1,j+1,bi,bj)       &           + GGL90visctmp(i-1,j+1,k)*mskCor(i-1,j+1,bi,bj) ) )
694       &       + GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj))       &             )/(
695       &  )       &     p4
696       & /(p4       &    +p8 *( (  maskC(i-1,j  ,k,bi,bj)*mskCor(i-1,j  ,bi,bj)
697       &  +p8 *(       maskC(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)       &           +  maskC(i+1,j  ,k,bi,bj)*mskCor(i+1,j  ,bi,bj) )
698       &       +       maskC(i  ,j-1,k,bi,bj) * mskCor(i  ,j-1,bi,bj)       &         + (  maskC(i  ,j-1,k,bi,bj)*mskCor(i  ,j-1,bi,bj)
699       &       +       maskC(i+1,j  ,k,bi,bj) * mskCor(i+1,j  ,bi,bj)       &           +  maskC(i  ,j+1,k,bi,bj)*mskCor(i  ,j+1,bi,bj) ) )
700       &       +       maskC(i  ,j+1,k,bi,bj) * mskCor(i  ,j+1,bi,bj))       &    +p16*( (  maskC(i+1,j+1,k,bi,bj)* mskCor(i+1,j+1,bi,bj)
701       &  +p16*(       maskC(i+1,j+1,k,bi,bj) * mskCor(i+1,j+1,bi,bj)       &           +  maskC(i-1,j-1,k,bi,bj)*mskCor(i-1,j-1,bi,bj) )
702       &       +       maskC(i+1,j-1,k,bi,bj) * mskCor(i+1,j-1,bi,bj)       &         + (  maskC(i+1,j-1,k,bi,bj)*mskCor(i+1,j-1,bi,bj)
703       &       +       maskC(i-1,j+1,k,bi,bj) * mskCor(i-1,j+1,bi,bj)       &           +  maskC(i-1,j+1,k,bi,bj)*mskCor(i-1,j+1,bi,bj) ) )
704       &       +       maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj))       &               )*maskC(i,j,k,bi,bj)*mskCor(i,j,bi,bj)
      &  )*maskC(i,j,k,bi,bj)*mskCor(i,j,bi,bj)  
705  #else  #else
706           tmpVisc = GGL90visctmp(i,j,k)           tmpVisc = GGL90visctmp(i,j,k)
707  #endif  #endif
708           tmpVisc = MIN(tmpVisc/TKEPrandtlNumber(i,j,k),GGL90diffMax)           tmpVisc = MIN(tmpVisc/TKEPrandtlNumber(i,j,k),GGL90diffMax)
709           GGL90diffKr(i,j,k,bi,bj)= MAX( tmpVisc , diffKrNrT(k) )           GGL90diffKr(i,j,k,bi,bj)= MAX( tmpVisc , diffKrNrS(k) )
710          ENDDO          ENDDO
711         ENDDO         ENDDO
712        ENDDO        ENDDO
# Line 589  C     =============================== Line 715  C     ===============================
715         DO j=1,sNy         DO j=1,sNy
716          DO i=1,sNx+1          DO i=1,sNx+1
717  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
718          tmpVisc =           tmpVisc = (
719       & (       &     p4 *(   GGL90visctmp(i-1,j  ,k)*mskCor(i-1,j  ,bi,bj)
720       &   p4 *(GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)       &           + GGL90visctmp(i  ,j  ,k)*mskCor(i  ,j  ,bi,bj) )
721       &       +GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj))       &    +p8 *( ( GGL90visctmp(i-1,j-1,k)*mskCor(i-1,j-1,bi,bj)
722       &  +p8 *(GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj)       &           + GGL90visctmp(i  ,j-1,k)*mskCor(i  ,j-1,bi,bj) )
723       &       +GGL90visctmp(i-1,j+1,k) * mskCor(i-1,j+1,bi,bj)       &         + ( GGL90visctmp(i-1,j+1,k)*mskCor(i-1,j+1,bi,bj)
724       &       +GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj)       &           + GGL90visctmp(i  ,j+1,k)*mskCor(i  ,j+1,bi,bj) ) )
725       &       +GGL90visctmp(i  ,j+1,k) * mskCor(i  ,j+1,bi,bj))       &             )/(
726       &  )       &     p4 * 2. _d 0
727       & /(p4 * 2. _d 0       &    +p8 *( (  maskC(i-1,j-1,k,bi,bj)*mskCor(i-1,j-1,bi,bj)
728       &  +p8 *(      maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj)       &           +  maskC(i  ,j-1,k,bi,bj)*mskCor(i  ,j-1,bi,bj) )
729       &       +      maskC(i-1,j+1,k,bi,bj) * mskCor(i-1,j+1,bi,bj)       &         + (  maskC(i-1,j+1,k,bi,bj)*mskCor(i-1,j+1,bi,bj)
730       &       +      maskC(i  ,j-1,k,bi,bj) * mskCor(i  ,j-1,bi,bj)       &           +  maskC(i  ,j+1,k,bi,bj)*mskCor(i  ,j+1,bi,bj) ) )
731       &       +      maskC(i  ,j+1,k,bi,bj) * mskCor(i  ,j+1,bi,bj))       &               )*maskC(i-1,j,k,bi,bj)*mskCor(i-1,j,bi,bj)
732       &  )       &                *maskC(i  ,j,k,bi,bj)*mskCor(i  ,j,bi,bj)
      &  *maskC(i  ,j,k,bi,bj)*mskCor(i  ,j,bi,bj)  
      &  *maskC(i-1,j,k,bi,bj)*mskCor(i-1,j,bi,bj)  
733  #else  #else
734          tmpVisc = _maskW(i,j,k,bi,bj) *           tmpVisc = _maskW(i,j,k,bi,bj) * halfRL
735       &                   (.5 _d 0*(GGL90visctmp(i,j,k)       &          *( GGL90visctmp(i-1,j,k)
736       &                            +GGL90visctmp(i-1,j,k))       &           + GGL90visctmp(i,j,k) )
      &                   )  
737  #endif  #endif
738          tmpVisc = MIN( tmpVisc , GGL90viscMax )           tmpVisc = MIN( tmpVisc , GGL90viscMax )
739          GGL90viscArU(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )           GGL90viscArU(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
740          ENDDO          ENDDO
741         ENDDO         ENDDO
742        ENDDO        ENDDO
# Line 622  C     =============================== Line 745  C     ===============================
745         DO j=1,sNy+1         DO j=1,sNy+1
746          DO i=1,sNx          DO i=1,sNx
747  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
748          tmpVisc =           tmpVisc = (
749       & (       &     p4 *(   GGL90visctmp(i  ,j-1,k)*mskCor(i  ,j-1,bi,bj)
750       &   p4 *(GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)       &           + GGL90visctmp(i  ,j  ,k)*mskCor(i  ,j  ,bi,bj) )
751       &       +GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj))       &    +p8 *( ( GGL90visctmp(i-1,j-1,k)*mskCor(i-1,j-1,bi,bj)
752       &  +p8 *(GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj)       &           + GGL90visctmp(i-1,j  ,k)*mskCor(i-1,j  ,bi,bj) )
753       &       +GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj)       &         + ( GGL90visctmp(i+1,j-1,k)*mskCor(i+1,j-1,bi,bj)
754       &       +GGL90visctmp(i+1,j  ,k) * mskCor(i+1,j  ,bi,bj)       &           + GGL90visctmp(i+1,j  ,k)*mskCor(i+1,j  ,bi,bj) ) )
755       &       +GGL90visctmp(i+1,j-1,k) * mskCor(i+1,j-1,bi,bj))       &             )/(
756       &  )       &     p4 * 2. _d 0
757       & /(p4 * 2. _d 0       &    +p8 *( (  maskC(i-1,j-1,k,bi,bj)*mskCor(i-1,j-1,bi,bj)
758       &  +p8 *(      maskC(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)       &           +  maskC(i-1,j  ,k,bi,bj)*mskCor(i-1,j  ,bi,bj) )
759       &       +      maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj)       &         + (  maskC(i+1,j-1,k,bi,bj)*mskCor(i+1,j-1,bi,bj)
760       &       +      maskC(i+1,j  ,k,bi,bj) * mskCor(i+1,j  ,bi,bj)       &           +  maskC(i+1,j  ,k,bi,bj)*mskCor(i+1,j  ,bi,bj) ) )
761       &       +      maskC(i+1,j-1,k,bi,bj) * mskCor(i+1,j-1,bi,bj))       &               )*maskC(i,j-1,k,bi,bj)*mskCor(i,j-1,bi,bj)
762       &  )       &                *maskC(i,j  ,k,bi,bj)*mskCor(i,j  ,bi,bj)
      &   *maskC(i,j  ,k,bi,bj)*mskCor(i,j  ,bi,bj)  
      &   *maskC(i,j-1,k,bi,bj)*mskCor(i,j-1,bi,bj)  
763  #else  #else
764          tmpVisc = _maskS(i,j,k,bi,bj) *           tmpVisc = _maskS(i,j,k,bi,bj) * halfRL
765       &                   (.5 _d 0*(GGL90visctmp(i,j,k)       &          *( GGL90visctmp(i,j-1,k)
766       &                            +GGL90visctmp(i,j-1,k))       &           + GGL90visctmp(i,j,k) )
      &                   )  
   
767  #endif  #endif
768          tmpVisc = MIN( tmpVisc , GGL90viscMax )           tmpVisc = MIN( tmpVisc , GGL90viscMax )
769          GGL90viscArV(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )           GGL90viscArV(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
770          ENDDO          ENDDO
771         ENDDO         ENDDO
772        ENDDO        ENDDO
773    
774  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
775        IF ( useDiagnostics ) THEN        IF ( useDiagnostics ) THEN
776           CALL DIAGNOSTICS_FILL( GGL90TKE   ,'GGL90TKE',          CALL DIAGNOSTICS_FILL( GGL90TKE   ,'GGL90TKE',
777       &                          0,Nr, 1, bi, bj, myThid )       &                         0,Nr, 1, bi, bj, myThid )
778           CALL DIAGNOSTICS_FILL( GGL90viscArU,'GGL90ArU',          CALL DIAGNOSTICS_FILL( GGL90viscArU,'GGL90ArU',
779       &                          0,Nr, 1, bi, bj, myThid )       &                         0,Nr, 1, bi, bj, myThid )
780           CALL DIAGNOSTICS_FILL( GGL90viscArV,'GGL90ArV',          CALL DIAGNOSTICS_FILL( GGL90viscArV,'GGL90ArV',
781       &                          0,Nr, 1, bi, bj, myThid )       &                         0,Nr, 1, bi, bj, myThid )
782           CALL DIAGNOSTICS_FILL( GGL90diffKr,'GGL90Kr ',          CALL DIAGNOSTICS_FILL( GGL90diffKr,'GGL90Kr ',
783       &                          0,Nr, 1, bi, bj, myThid )       &                         0,Nr, 1, bi, bj, myThid )
784           CALL DIAGNOSTICS_FILL( TKEPrandtlNumber ,'GGL90Prl',          CALL DIAGNOSTICS_FILL( TKEPrandtlNumber ,'GGL90Prl',
785       &                          0,Nr, 2, bi, bj, myThid )       &                         0,Nr, 2, bi, bj, myThid )
786           CALL DIAGNOSTICS_FILL( GGL90mixingLength,'GGL90Lmx',          CALL DIAGNOSTICS_FILL( GGL90mixingLength,'GGL90Lmx',
787       &                          0,Nr, 2, bi, bj, myThid )       &                         0,Nr, 2, bi, bj, myThid )
788    
789            kp1 = MIN(Nr,kSurf+1)
790            DO j=jMin,jMax
791             DO i=iMin,iMax
792    C     diagnose surface flux of TKE
793              surf_flx_tke(i,j) =(GGL90TKE(i,j,kSurf,bi,bj)-
794         &                        GGL90TKE(i,j,kp1,bi,bj))
795         &        *recip_drF(kSurf)*recip_hFacC(i,j,kSurf,bi,bj)
796         &        *KappaE(i,j,kp1)
797             ENDDO
798            ENDDO
799            CALL DIAGNOSTICS_FILL( surf_flx_tke,'GGL90flx',
800         &                         0, 1, 2, bi, bj, myThid )
801    
802            k=kSurf
803            DO j=jMin,jMax
804             DO i=iMin,iMax
805    C     diagnose work done by the wind
806              surf_flx_tke(i,j) =
807         &      halfRL*( surfaceForcingU(i,  j,bi,bj)*uVel(i  ,j,k,bi,bj)
808         &              +surfaceForcingU(i+1,j,bi,bj)*uVel(i+1,j,k,bi,bj))
809         &    + halfRL*( surfaceForcingV(i,j,  bi,bj)*vVel(i,j  ,k,bi,bj)
810         &              +surfaceForcingV(i,j+1,bi,bj)*vVel(i,j+1,k,bi,bj))
811             ENDDO
812            ENDDO
813            CALL DIAGNOSTICS_FILL( surf_flx_tke,'GGL90tau',
814         &                         0, 1, 2, bi, bj, myThid )
815    
816        ENDIF        ENDIF
817  #endif  #endif /* ALLOW_DIAGNOSTICS */
818    
819  #endif /* ALLOW_GGL90 */  #endif /* ALLOW_GGL90 */
820    
Legend:
Removed from v.1.19  
changed lines
  Added in v.1.34
  ViewVC Help
Powered by ViewVC 1.1.22