/[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.1 by mlosch, Thu Sep 16 11:27:18 2004 UTC revision 1.31 by jmc, Mon Feb 23 21:20:15 2015 UTC
# Line 7  CBOP Line 7  CBOP
7  C !ROUTINE: GGL90_CALC  C !ROUTINE: GGL90_CALC
8    
9  C !INTERFACE: ======================================================  C !INTERFACE: ======================================================
10        subroutine GGL90_CALC(        SUBROUTINE GGL90_CALC(
11       I     bi, bj, myTime, myThid )       I                 bi, bj, sigmaR, myTime, myIter, myThid )
12    
13  C !DESCRIPTION: \bv  C !DESCRIPTION: \bv
14  C     /==========================================================\  C     *==========================================================*
15  C     | SUBROUTINE GGL90_CALC                                    |  C     | SUBROUTINE GGL90_CALC                                    |
16  C     | o Compute all GGL90 fields defined in GGL90.h            |  C     | o Compute all GGL90 fields defined in GGL90.h            |
17  C     |==========================================================|  C     *==========================================================*
18  C     | Equation numbers refer to                                |  C     | Equation numbers refer to                                |
19  C     | Gaspar et al. (1990), JGR 95 (C9), pp 16,179             |  C     | Gaspar et al. (1990), JGR 95 (C9), pp 16,179             |
20  C     | Some parts of the implementation follow Blanke and       |  C     | Some parts of the implementation follow Blanke and       |
21  C     | Delecuse (1993), JPO, and OPA code, in particular the    |  C     | Delecuse (1993), JPO, and OPA code, in particular the    |
22  C     | computation of the                                       |  C     | computation of the                                       |
23  C     | mixing length = max(min(lk,depth),lkmin)                 |  C     | mixing length = max(min(lk,depth),lkmin)                 |
24  C     \==========================================================/  C     *==========================================================*
       IMPLICIT NONE  
 C  
 C--------------------------------------------------------------------  
25    
26  C global parameters updated by ggl90_calc  C global parameters updated by ggl90_calc
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: ============================================================
33          IMPLICIT NONE
34  #include "SIZE.h"  #include "SIZE.h"
35  #include "EEPARAMS.h"  #include "EEPARAMS.h"
36  #include "PARAMS.h"  #include "PARAMS.h"
# Line 43  C !USES: =============================== Line 40  C !USES: ===============================
40  #include "GRID.h"  #include "GRID.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     myTime - Current time in simulation  C     sigmaR :: Vertical gradient of iso-neutral density
46    C     myTime :: Current time in simulation
47    C     myIter :: Current time-step number
48    C     myThid :: My Thread Id number
49        INTEGER bi, bj        INTEGER bi, bj
50        INTEGER myThid        _RL     sigmaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
51        _RL     myTime        _RL     myTime
52          INTEGER myIter
53          INTEGER myThid
54    
55  #ifdef ALLOW_GGL90  #ifdef ALLOW_GGL90
56    
57  C !LOCAL VARIABLES: ====================================================  C !LOCAL VARIABLES: ====================================================
58  c Local constants  C Local constants
59  C     iMin, iMax, jMin, jMax, I, J - array computation indices  C     iMin,iMax,jMin,jMax :: index boundaries of computation domain
60  C     K, Kp1, km1, kSurf, kBottom  - vertical loop indices  C     i, j, k, kp1,km1 :: array computation indices
61  C     ab15, ab05       - weights for implicit timestepping  C     kSurf, kBottom   :: vertical indices of domain boundaries
62  C     uStarSquare      - square of friction velocity  C     explDissFac      :: explicit Dissipation Factor (in [0-1])
63  C     verticalShear    - (squared) vertical shear of horizontal velocity  C     implDissFac      :: implicit Dissipation Factor (in [0-1])
64  C     Nsquare          - squared buoyancy freqency  C     uStarSquare      :: square of friction velocity
65  C     RiNumber         - local Richardson number  C     verticalShear    :: (squared) vertical shear of horizontal velocity
66  C     KappaM           - (local) viscosity parameter (eq.10)  C     Nsquare          :: squared buoyancy freqency
67  C     KappaH           - (local) diffusivity parameter for temperature (eq.11)  C     RiNumber         :: local Richardson number
68  C     KappaE           - (local) diffusivity parameter for TKE (eq.15)  C     KappaM           :: (local) viscosity parameter (eq.10)
69  C     buoyFreq         - buoyancy freqency  C     KappaH           :: (local) diffusivity parameter for temperature (eq.11)
70  C     TKEdissipation   - dissipation of TKE  C     KappaE           :: (local) diffusivity parameter for TKE (eq.15)
71  C     GGL90mixingLength- mixing length of scheme following Banke+Delecuse  C     TKEdissipation   :: dissipation of TKE
72  C     totalDepth       - thickness of water column (inverse of recip_Rcol)  C     GGL90mixingLength:: mixing length of scheme following Banke+Delecuse
73  C     TKEPrandtlNumber - here, an empirical function of the Richardson number  C         rMixingLength:: inverse of mixing length
74  C     rhoK, rhoKm1     - density at layer K and Km1 (relative to K)  C     totalDepth       :: thickness of water column (inverse of recip_Rcol)
75  C     gTKE             - right hand side of implicit equation  C     TKEPrandtlNumber :: here, an empirical function of the Richardson number
76        INTEGER iMin ,iMax ,jMin ,jMax        INTEGER iMin ,iMax ,jMin ,jMax
77        INTEGER I, J, K, Kp1, Km1, kSurf, kBottom        INTEGER i, j, k, kp1, km1, kSurf, kBottom
78        _RL     ab15, ab05        _RL     explDissFac, implDissFac
79        _RL     uStarSquare        _RL     uStarSquare
80        _RL     verticalShear        _RL     verticalShear
81        _RL     KappaM, KappaH        _RL     KappaM, KappaH
82        _RL     Nsquare  c     _RL     Nsquare
83          _RL     Nsquare(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
84        _RL     deltaTggl90        _RL     deltaTggl90
85        _RL     SQRTTKE  c     _RL     SQRTTKE
86          _RL     SQRTTKE(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
87        _RL     RiNumber        _RL     RiNumber
88    #ifdef ALLOW_GGL90_IDEMIX
89          _RL     IDEMIX_RiNumber
90    #endif
91        _RL     TKEdissipation        _RL     TKEdissipation
92        _RL     tempU, tempV, prTemp        _RL     tempU, tempV, prTemp
93          _RL     MaxLength, tmpmlx, tmpVisc
94        _RL     TKEPrandtlNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     TKEPrandtlNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
95        _RL     GGL90mixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     GGL90mixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
96          _RL     rMixingLength    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
97          _RL     mxLength_Dn      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
98        _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)  
99        _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
100        _RL     gTKE             (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     GGL90visctmp     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
101  C     tri-diagonal matrix  #ifdef ALLOW_DIAGNOSTICS
102        _RL     a(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     surf_flx_tke     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
103        _RL     b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  #endif /* ALLOW_DIAGNOSTICS */
104        _RL     c(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  C-    tri-diagonal matrix
105          _RL     a3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
106          _RL     b3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
107          _RL     c3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
108          INTEGER errCode
109    #ifdef ALLOW_GGL90_HORIZDIFF
110    C     hFac     :: fractional thickness of W-cell
111    C     xA, yA   :: area of lateral faces
112    C     dfx, dfy :: diffusive flux across lateral faces
113    C     gTKE     :: right hand side of diffusion equation
114          _RL     hFac
115          _RL     xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
116          _RL     yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
117          _RL     dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
118          _RL     dfy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
119          _RL    gTKE(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
120    #endif /* ALLOW_GGL90_HORIZDIFF */
121    #ifdef ALLOW_GGL90_SMOOTH
122          _RL p4, p8, p16
123  CEOP  CEOP
124          p4=0.25 _d 0
125          p8=0.125 _d 0
126          p16=0.0625 _d 0
127    #endif
128        iMin = 2-OLx        iMin = 2-OLx
129        iMax = sNx+OLx-1        iMax = sNx+OLx-1
130        jMin = 2-OLy        jMin = 2-OLy
131        jMax = sNy+OLy-1        jMax = sNy+OLy-1
132    
133  C     set separate time step (should be deltaTtracer)  C     set separate time step (should be deltaTtracer)
134        deltaTggl90 = deltaTtracer        deltaTggl90 = dTtracerLev(1)
135  C      
136        kSurf = 1        kSurf = 1
137  C     implicit timestepping weights for dissipation  C     explicit/implicit timestepping weights for dissipation
138        ab15 =  1.5 _d 0        explDissFac = 0. _d 0
139        ab05 = -0.5 _d 0        implDissFac = 1. _d 0 - explDissFac
       ab15 =  1. _d 0  
       ab05 =  0. _d 0  
140    
141  C     Initialize local fields  C     Initialize local fields
142        DO K = 1, Nr        DO k = 1, Nr
143         DO J=1-Oly,sNy+Oly         DO j=1-OLy,sNy+OLy
144          DO I=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
145           gTKE(I,J,K)              = 0. _d 0           rMixingLength(i,j,k)     = 0. _d 0
146           KappaE(I,J,K)            = 0. _d 0           mxLength_Dn(i,j,k)       = 0. _d 0
147           TKEPrandtlNumber(I,J,K)  = 0. _d 0           GGL90visctmp(i,j,k)      = 0. _d 0
148           GGL90mixingLength(I,J,K) = 0. _d 0           KappaE(i,j,k)            = 0. _d 0
149          ENDDO           TKEPrandtlNumber(i,j,k)  = 1. _d 0
150         ENDDO               GGL90mixingLength(i,j,k) = GGL90mixingLengthMin
151        ENDDO           GGL90visctmp(i,j,k)      = 0. _d 0
152        DO J=1-Oly,sNy+Oly  #ifndef SOLVE_DIAGONAL_LOWMEMORY
153         DO I=1-Olx,sNx+Olx           a3d(i,j,k) = 0. _d 0
154          rhoK    (I,J) = 0. _d 0           b3d(i,j,k) = 1. _d 0
155          rhoKm1  (I,J) = 0. _d 0           c3d(i,j,k) = 0. _d 0
156          totalDepth(I,J)   = 0. _d 0  #endif
157          IF ( recip_Rcol(I,J,bi,bj) .NE. 0. )           Nsquare(i,j,k) = 0. _d 0
158       &       totalDepth(I,J) = 1./recip_Rcol(I,J,bi,bj)           SQRTTKE(i,j,k) = 0. _d 0
159            ENDDO
160         ENDDO         ENDDO
161        ENDDO        ENDDO
162          DO j=1-OLy,sNy+OLy
163           DO i=1-OLx,sNx+OLx
164            totalDepth(i,j)    = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
165            rMixingLength(i,j,1) = 0. _d 0
166            mxLength_Dn(i,j,1) = GGL90mixingLengthMin
167            SQRTTKE(i,j,1) = SQRT( GGL90TKE(i,j,1,bi,bj) )
168    #ifdef ALLOW_GGL90_HORIZDIFF
169            xA(i,j)  = 0. _d 0
170            yA(i,j)  = 0. _d 0
171            dfx(i,j) = 0. _d 0
172            dfy(i,j) = 0. _d 0
173            gTKE(i,j) = 0. _d 0
174    #endif /* ALLOW_GGL90_HORIZDIFF */
175           ENDDO
176          ENDDO
177    
178    #ifdef ALLOW_GGL90_IDEMIX
179          IF ( useIDEMIX) CALL GGL90_IDEMIX(
180         &   bi, bj, sigmaR, myTime, myIter, myThid )
181    #endif /* ALLOW_GGL90_IDEMIX */
182    
183  C     start k-loop  C     start k-loop
184        DO K = 2, Nr        DO k = 2, Nr
185         Km1 = K-1  c      km1 = k-1
186         Kp1 = MIN(Nr,K+1)  c      kp1 = MIN(Nr,k+1)
187         CALL FIND_RHO(         DO j=jMin,jMax
188       I      bi, bj, iMin, iMax, jMin, jMax, Km1, K,          DO i=iMin,iMax
189       I      theta, salt,           SQRTTKE(i,j,k)=SQRT( GGL90TKE(i,j,k,bi,bj) )
190       O      rhoKm1,  
      I      myThid )  
        CALL FIND_RHO(  
      I      bi, bj, iMin, iMax, jMin, jMax, K, K,  
      I      theta, salt,  
      O      rhoK,  
      I      myThid )  
        DO J=jMin,jMax  
         DO I=iMin,iMax  
          SQRTTKE=SQRT( GGL90TKE(I,J,K,bi,bj) )  
 C  
191  C     buoyancy frequency  C     buoyancy frequency
192  C           Nsquare(i,j,k) = gravity*gravitySign*recip_rhoConst
193           Nsquare = - gravity*recip_rhoConst*recip_drC(K)       &                  * sigmaR(i,j,k)
194       &        * ( rhoKm1(I,J) - rhoK(I,J) )*maskC(I,J,K,bi,bj)  cC     vertical shear term (dU/dz)^2+(dV/dz)^2
195    c         tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)
196    c     &                 -( uVel(i,j,k  ,bi,bj)+uVel(i+1,j,k  ,bi,bj)) )
197    c     &        *recip_drC(k)
198    c         tempV= .5 _d 0*( vVel(i,j,km1,bi,bj)+vVel(i,j+1,km1,bi,bj)
199    c     &                 -( vVel(i,j,k  ,bi,bj)+vVel(i,j+1,k  ,bi,bj)) )
200    c     &        *recip_drC(k)
201    c         verticalShear = tempU*tempU + tempV*tempV
202    c         RiNumber   = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)
203    cC     compute Prandtl number (always greater than 0)
204    c         prTemp = 1. _d 0
205    c         IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
206    c         TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
207    C     mixing length
208             GGL90mixingLength(i,j,k) = SQRTTWO *
209         &        SQRTTKE(i,j,k)/SQRT( MAX(Nsquare(i,j,k),GGL90eps) )
210            ENDDO
211           ENDDO
212          ENDDO
213    
214    C- ensure mixing between first and second level
215          IF (mxlSurfFlag) THEN
216           DO j=jMin,jMax
217            DO i=iMin,iMax
218             GGL90mixingLength(i,j,2)=drF(1)
219            ENDDO
220           ENDDO
221          ENDIF
222    
223    C- Impose upper and lower bound for mixing length
224          IF ( mxlMaxFlag .EQ. 0 ) THEN
225    
226           DO k=2,Nr
227            DO j=jMin,jMax
228             DO i=iMin,iMax
229              MaxLength=totalDepth(i,j)
230              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
231         &                                   MaxLength)
232             ENDDO
233            ENDDO
234           ENDDO
235    
236           DO k=2,Nr
237            DO j=jMin,jMax
238             DO i=iMin,iMax
239              GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
240         &                                   GGL90mixingLengthMin)
241              rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
242             ENDDO
243            ENDDO
244           ENDDO
245    
246          ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN
247    
248           DO k=2,Nr
249            DO j=jMin,jMax
250             DO i=iMin,iMax
251              MaxLength=MIN(Ro_surf(i,j,bi,bj)-rF(k),rF(k)-R_low(i,j,bi,bj))
252    c         MaxLength=MAX(MaxLength,20. _d 0)
253              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
254         &                                   MaxLength)
255             ENDDO
256            ENDDO
257           ENDDO
258    
259           DO k=2,Nr
260            DO j=jMin,jMax
261             DO i=iMin,iMax
262              GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
263         &                                   GGL90mixingLengthMin)
264              rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
265             ENDDO
266            ENDDO
267           ENDDO
268    
269          ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN
270    
271           DO k=2,Nr
272            DO j=jMin,jMax
273             DO i=iMin,iMax
274              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
275         &        GGL90mixingLength(i,j,k-1)+drF(k-1))
276             ENDDO
277            ENDDO
278           ENDDO
279           DO j=jMin,jMax
280            DO i=iMin,iMax
281              GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
282         &       GGL90mixingLengthMin+drF(Nr))
283            ENDDO
284           ENDDO
285           DO k=Nr-1,2,-1
286            DO j=jMin,jMax
287             DO i=iMin,iMax
288              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
289         &        GGL90mixingLength(i,j,k+1)+drF(k))
290             ENDDO
291            ENDDO
292           ENDDO
293    
294           DO k=2,Nr
295            DO j=jMin,jMax
296             DO i=iMin,iMax
297              GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
298         &                                   GGL90mixingLengthMin)
299              rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
300             ENDDO
301            ENDDO
302           ENDDO
303    
304          ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN
305    
306           DO k=2,Nr
307            DO j=jMin,jMax
308             DO i=iMin,iMax
309              mxLength_Dn(i,j,k) = MIN(GGL90mixingLength(i,j,k),
310         &        mxLength_Dn(i,j,k-1)+drF(k-1))
311             ENDDO
312            ENDDO
313           ENDDO
314           DO j=jMin,jMax
315            DO i=iMin,iMax
316              GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
317         &       GGL90mixingLengthMin+drF(Nr))
318            ENDDO
319           ENDDO
320           DO k=Nr-1,2,-1
321            DO j=jMin,jMax
322             DO i=iMin,iMax
323              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
324         &        GGL90mixingLength(i,j,k+1)+drF(k))
325             ENDDO
326            ENDDO
327           ENDDO
328    
329           DO k=2,Nr
330            DO j=jMin,jMax
331             DO i=iMin,iMax
332              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
333         &                                  mxLength_Dn(i,j,k))
334              tmpmlx = SQRT( GGL90mixingLength(i,j,k)*mxLength_Dn(i,j,k) )
335              tmpmlx = MAX( tmpmlx, GGL90mixingLengthMin)
336              rMixingLength(i,j,k) = 1. _d 0 / tmpmlx
337             ENDDO
338            ENDDO
339           ENDDO
340    
341          ELSE
342           STOP 'GGL90_CALC: Wrong mxlMaxFlag (mixing length limit)'
343          ENDIF
344    
345    C- Impose minimum mixing length (to avoid division by zero)
346    c      DO k=2,Nr
347    c      DO j=jMin,jMax
348    c       DO i=iMin,iMax
349    c        GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
350    c    &        GGL90mixingLengthMin)
351    c        rMixingLength(i,j,k) = 1. _d 0 /GGL90mixingLength(i,j,k)
352    c       ENDDO
353    c      ENDDO
354    c     ENDDO
355    
356          DO k=2,Nr
357           km1 = k-1
358    
359    #ifdef ALLOW_GGL90_HORIZDIFF
360          IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
361    C     horizontal diffusion of TKE (requires an exchange in
362    C      do_fields_blocking_exchanges)
363    C     common factors
364            DO j=1-OLy,sNy+OLy
365             DO i=1-OLx,sNx+OLx
366              xA(i,j) = _dyG(i,j,bi,bj)*drC(k)*
367         &                 (min(.5 _d 0,_hFacW(i,j,k-1,bi,bj) ) +
368         &                  min(.5 _d 0,_hFacW(i,j,k  ,bi,bj) ) )
369              yA(i,j) = _dxG(i,j,bi,bj)*drC(k)*
370         &                 (min(.5 _d 0,_hFacS(i,j,k-1,bi,bj) ) +
371         &                  min(.5 _d 0,_hFacS(i,j,k  ,bi,bj) ) )
372             ENDDO
373            ENDDO
374    C     Compute diffusive fluxes
375    C     ... across x-faces
376            DO j=1-OLy,sNy+OLy
377             dfx(1-OLx,j)=0. _d 0
378             DO i=1-OLx+1,sNx+OLx
379              dfx(i,j) = -GGL90diffTKEh*xA(i,j)
380         &      *_recip_dxC(i,j,bi,bj)
381         &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i-1,j,k,bi,bj))
382    #ifdef ISOTROPIC_COS_SCALING
383         &      *CosFacU(j,bi,bj)
384    #endif /* ISOTROPIC_COS_SCALING */
385             ENDDO
386            ENDDO
387    C     ... across y-faces
388            DO i=1-OLx,sNx+OLx
389             dfy(i,1-OLy)=0. _d 0
390            ENDDO
391            DO j=1-OLy+1,sNy+OLy
392             DO i=1-OLx,sNx+OLx
393              dfy(i,j) = -GGL90diffTKEh*yA(i,j)
394         &      *_recip_dyC(i,j,bi,bj)
395         &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i,j-1,k,bi,bj))
396    #ifdef ISOTROPIC_COS_SCALING
397         &      *CosFacV(j,bi,bj)
398    #endif /* ISOTROPIC_COS_SCALING */
399             ENDDO
400            ENDDO
401    C     Compute divergence of fluxes
402            DO j=1-OLy,sNy+OLy-1
403             DO i=1-OLx,sNx+OLx-1
404    #ifdef ALLOW_GGL90_IDEMIX
405              gTKE(i,j) = -recip_drC(k)*recip_rA(i,j,bi,bj)
406         &         *recip_hFacI(i,j,k,bi,bj)
407    #else
408              hFac = MIN(.5 _d 0,_hFacC(i,j,k-1,bi,bj) ) +
409         &           MIN(.5 _d 0,_hFacC(i,j,k  ,bi,bj) )
410              gTKE(i,j) = 0.0
411              IF ( hFac .ne. 0.0 )
412         &      gTKE(i,j) = -recip_drC(k)*recip_rA(i,j,bi,bj)/hFac
413    #endif
414         &         *((dfx(i+1,j)-dfx(i,j))
415         &          +(dfy(i,j+1)-dfy(i,j)) )
416             ENDDO
417            ENDDO
418    C      end if GGL90diffTKEh .eq. 0.
419           ENDIF
420    #endif /* ALLOW_GGL90_HORIZDIFF */
421    
422           DO j=jMin,jMax
423            DO i=iMin,iMax
424  C     vertical shear term (dU/dz)^2+(dV/dz)^2  C     vertical shear term (dU/dz)^2+(dV/dz)^2
425           tempu= .5*(  uVel(I,J,Km1,bi,bj)+uVel(I+1,J,Km1,bi,bj)           tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)
426       &             - (uVel(I,J,K  ,bi,bj)+uVel(I+1,J,K  ,bi,bj)) )       &                 -( uVel(i,j,k  ,bi,bj)+uVel(i+1,j,k  ,bi,bj)) )
427       &        *recip_drC(K)       &        *recip_drC(k)
428           tempv= .5*(  vVel(I,J,Km1,bi,bj)+vVel(I,J+1,Km1,bi,bj)  #ifdef ALLOW_GGL90_IDEMIX
429       &             - (vVel(I,J,K  ,bi,bj)+vVel(I,J+1,K  ,bi,bj)) )       &        *recip_hFacI(i,j,k,bi,bj)
430       &        *recip_drC(K)  #endif
431             tempV= .5 _d 0*( vVel(i,j,km1,bi,bj)+vVel(i,j+1,km1,bi,bj)
432         &                 -( vVel(i,j,k  ,bi,bj)+vVel(i,j+1,k  ,bi,bj)) )
433         &        *recip_drC(k)
434    #ifdef ALLOW_GGL90_IDEMIX
435         &        *recip_hFacI(i,j,k,bi,bj)
436    #endif
437           verticalShear = tempU*tempU + tempV*tempV           verticalShear = tempU*tempU + tempV*tempV
438           RiNumber   = MAX(Nsquare,0. _d 0)/(verticalShear+GGL90eps)  
439    C     viscosity and diffusivity
440             KappaM = GGL90ck*GGL90mixingLength(i,j,k)*SQRTTKE(i,j,k)
441             GGL90visctmp(i,j,k) = MAX(KappaM,diffKrNrS(k))
442         &                            * maskC(i,j,k,bi,bj)
443    C        note: storing GGL90visctmp like this, and using it later to compute
444    C              GGL9rdiffKr etc. is robust in case of smoothing (e.g. see OPA)
445             KappaM = MAX(KappaM,viscArNr(k)) * maskC(i,j,k,bi,bj)
446    
447  C     compute Prandtl number (always greater than 0)  C     compute Prandtl number (always greater than 0)
448             RiNumber = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)
449    #ifdef ALLOW_GGL90_IDEMIX
450    CML         IDEMIX_RiNumber = 1./GGL90eps
451             IDEMIX_RiNumber = MAX( KappaM*Nsquare(i,j,k), 0. _d 0)/
452         &    (GGL90eps+IDEMIX_tau_d(i,j,k,bi,bj)*IDEMIX_E(i,j,k,bi,bj)**2)
453             prTemp         = MIN(5.*RiNumber, 6.6 _d 0*IDEMIX_RiNumber)
454    #else
455           prTemp = 1. _d 0           prTemp = 1. _d 0
456           IF ( RiNumber .GE. 0.2 ) prTemp = 5.0 * RiNumber           IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
457           TKEPrandtlNumber(I,J,K) = MIN(10.,prTemp)  #endif /* ALLOW_GGL90_IDEMIX */
458  C     mixing length           TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
459           GGL90mixingLength(I,J,K) =           TKEPrandtlNumber(i,j,k) = MAX( 1. _d 0,TKEPrandtlNumber(i,j,k))
460       &        SQRTTKE/SQRT( MAX(Nsquare,GGL90eps) )  
461  C     impose upper bound for mixing length (total depth)  C        diffusivity
462           GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),           KappaH = KappaM/TKEPrandtlNumber(i,j,k)
463       &        totalDepth(I,J))           KappaE(i,j,k) = GGL90alpha * KappaM * maskC(i,j,k,bi,bj)
464  C     impose minimum mixing length (to avoid division by zero)  
          GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),  
      &        GGL90mixingLengthMin)  
 C     viscosity of last timestep  
          KappaM = GGL90ck*GGL90mixingLength(I,J,K)*SQRTTKE  
          KappaE(I,J,K) = KappaM*GGL90alpha  
465  C     dissipation term  C     dissipation term
466           TKEdissipation = ab05*GGL90ceps           TKEdissipation = explDissFac*GGL90ceps
467       &        *SQRTTKE/GGL90mixingLength(I,J,K)       &        *SQRTTKE(i,j,k)*rMixingLength(i,j,k)
468       &        *GGL90TKE(I,J,K,bi,bj)             &        *GGL90TKE(i,j,k,bi,bj)
469  C     sum up contributions to form the right hand side  C     partial update with sum of explicit contributions
470           gTKE(I,J,K) = GGL90TKE(I,J,K,bi,bj)           GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
471       &        + deltaTggl90*(       &        + deltaTggl90*(
472       &        + KappaM*verticalShear       &        + KappaM*verticalShear
473       &        - KappaM*Nsquare/TKEPrandtlNumber(I,J,K)       &        - KappaH*Nsquare(i,j,k)
474       &        - TKEdissipation       &        - TKEdissipation
475    #ifdef ALLOW_GGL90_IDEMIX
476         &        + IDEMIX_tau_d(i,j,k,bi,bj)*IDEMIX_E(i,j,k,bi,bj)**2
477    #endif
478       &        )       &        )
479          ENDDO            ENDDO
480         ENDDO         ENDDO
481        ENDDO      
482  C      #ifdef ALLOW_GGL90_HORIZDIFF
483  C     Implicit time step to update TKE for k=1,Nr; TKE(Nr+1)=0 by default         IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
484  C      C--    Add horiz. diffusion tendency
485            DO j=jMin,jMax
486             DO i=iMin,iMax
487              GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
488         &                          + gTKE(i,j)*deltaTggl90
489             ENDDO
490            ENDDO
491           ENDIF
492    #endif /* ALLOW_GGL90_HORIZDIFF */
493    
494    C--   end of k loop
495          ENDDO
496    
497    C     ============================================
498    C     Implicit time step to update TKE for k=1,Nr;
499    C     TKE(Nr+1)=0 by default
500    C     ============================================
501  C     set up matrix  C     set up matrix
502  C--   Old aLower  C--   Lower diagonal
503        DO j=jMin,jMax        DO j=jMin,jMax
504         DO i=iMin,iMax         DO i=iMin,iMax
505           a(i,j,1) = 0. _d 0           a3d(i,j,1) = 0. _d 0
506         ENDDO         ENDDO
507        ENDDO        ENDDO
508        DO k=2,Nr        DO k=2,Nr
509         km1=MAX(1,k-1)         km1=MAX(2,k-1)
510         DO j=jMin,jMax         DO j=jMin,jMax
511          DO i=iMin,iMax          DO i=iMin,iMax
512           a(i,j,k) = -deltaTggl90  C-    We keep recip_hFacC in the diffusive flux calculation,
513       &        *recip_drF(km1)*recip_hFacI(i,j,k,bi,bj)  C-    but no hFacC in TKE volume control
514       &        *.5*(KappaE(i,j, k )+KappaE(i,j,km1))  C-    No need for maskC(k-1) with recip_hFacC(k-1)
515       &        *recip_drC(k)           a3d(i,j,k) = -deltaTggl90
516            IF (recip_hFacI(i,j,km1,bi,bj).EQ.0.) a(i,j,k)=0.       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)
517         &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))
518         &        *recip_drC(k)*maskC(i,j,k,bi,bj)
519    #ifdef ALLOW_GGL90_IDEMIX
520         &        *recip_hFacI(i,j,k,bi,bj)
521    #endif
522          ENDDO          ENDDO
523         ENDDO         ENDDO
524        ENDDO        ENDDO
525    C--   Upper diagonal
 C--   Old aUpper  
526        DO j=jMin,jMax        DO j=jMin,jMax
527         DO i=iMin,iMax         DO i=iMin,iMax
528           c(i,j,1)  = 0. _d 0           c3d(i,j,1)  = 0. _d 0
          c(i,j,Nr) = 0. _d 0  
529         ENDDO         ENDDO
530        ENDDO        ENDDO
 CML      DO k=1,Nr-1  
531        DO k=2,Nr        DO k=2,Nr
        kp1=min(Nr,k+1)  
532         DO j=jMin,jMax         DO j=jMin,jMax
533          DO i=iMin,iMax          DO i=iMin,iMax
534            c(i,j,k) = -deltaTggl90            kp1=MAX(1,MIN(klowC(i,j,bi,bj),k+1))
535       &        *recip_drF( k )*recip_hFacI(i,j,k,bi,bj)  C-    We keep recip_hFacC in the diffusive flux calculation,
536       &               *.5*(KappaE(i,j,k)+KappaE(i,j,kp1))  C-    but no hFacC in TKE volume control
537       &        *recip_drC(k)  C-    No need for maskC(k) with recip_hFacC(k)
538  C          IF (recip_hFacI(i,j,kp1,bi,bj).EQ.0.) c(i,j,k)=0.            c3d(i,j,k) = -deltaTggl90
539         &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)
540         &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))
541         &        *recip_drC(k)*maskC(i,j,k-1,bi,bj)
542    #ifdef ALLOW_GGL90_IDEMIX
543         &        *recip_hFacI(i,j,k,bi,bj)
544    #endif
545          ENDDO          ENDDO
546         ENDDO         ENDDO
547        ENDDO        ENDDO
548    
549  C--   Old aCenter        IF (.NOT.GGL90_dirichlet) THEN
550    C      Neumann bottom boundary condition for TKE: no flux from bottom
551           DO j=jMin,jMax
552            DO i=iMin,iMax
553             kBottom   = MAX(kLowC(i,j,bi,bj),1)
554             c3d(i,j,kBottom) = 0. _d 0
555            ENDDO
556           ENDDO
557          ENDIF
558    
559    C--   Center diagonal
560        DO k=1,Nr        DO k=1,Nr
561           km1 = MAX(k-1,1)
562         DO j=jMin,jMax         DO j=jMin,jMax
563          DO i=iMin,iMax          DO i=iMin,iMax
564            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)
565       &        + ab15*deltaTggl90*GGL90ceps*SQRT(GGL90TKE(I,J,K,bi,bj))       &        + implDissFac*deltaTggl90*GGL90ceps*SQRTTKE(i,j,k)
566       &        /GGL90mixingLength(I,J,K)       &        * rMixingLength(i,j,k)
567         &        * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)
568           ENDDO           ENDDO
569         ENDDO         ENDDO
570        ENDDO        ENDDO
571  C     end set up matrix  C     end set up matrix
572    
 C  
573  C     Apply boundary condition  C     Apply boundary condition
574  C            kp1 = MIN(Nr,kSurf+1)
575        DO J=jMin,jMax        DO j=jMin,jMax
576         DO I=iMin,iMax         DO i=iMin,iMax
577  C     estimate friction velocity uStar from surface forcing  C     estimate friction velocity uStar from surface forcing
578          uStarSquare = SQRT(          uStarSquare = SQRT(
579       &         ( .5*( surfaceForcingU(I,  J,  bi,bj)       &    ( .5 _d 0*( surfaceForcingU(i,  j,  bi,bj)
580       &              + surfaceForcingU(I+1,J,  bi,bj) ) )**2       &              + surfaceForcingU(i+1,j,  bi,bj) ) )**2
581       &       + ( .5*( surfaceForcingV(I,  J,  bi,bj)       &  + ( .5 _d 0*( surfaceForcingV(i,  j,  bi,bj)
582       &              + surfaceForcingV(I,  J+1,bi,bj) ) )**2       &              + surfaceForcingV(i,  j+1,bi,bj) ) )**2
583       &                     )       &                     )
584  C     Dirichlet surface boundary condition for TKE  C     Dirichlet surface boundary condition for TKE
585          gTKE(I,J,kSurf) = MAX(GGL90TKEmin,GGL90m2*uStarSquare)          GGL90TKE(i,j,kSurf,bi,bj) = maskC(i,j,kSurf,bi,bj)
586       &                     *maskC(I,J,kSurf,bi,bj)       &           *MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)
587  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom          GGL90TKE(i,j,kp1,bi,bj) = GGL90TKE(i,j,kp1,bi,bj)
588          kBottom   = MIN(MAX(kLowC(I,J,bi,bj),1),Nr)       &               - a3d(i,j,kp1)*GGL90TKE(i,j,kSurf,bi,bj)
589          gTKE(I,J,kBottom) = gTKE(I,J,kBottom)          a3d(i,j,kp1) = 0. _d 0
590       &       - GGL90TKEbottom*c(I,J,kBottom)         ENDDO
591         ENDDO        ENDDO
592        ENDDO      
593  C        IF (GGL90_dirichlet) THEN
594  C     solve tri-diagonal system, and store solution on gTKE (previously rhs)  C      Dirichlet bottom boundary condition for TKE = GGL90TKEbottom
595  C         DO j=jMin,jMax
596        CALL GGL90_SOLVE( bi, bj, iMin, iMax, jMin, jMax,          DO i=iMin,iMax
597       I     a, b, c,           kBottom   = MAX(kLowC(i,j,bi,bj),1)
598       U     gTKE,           GGL90TKE(i,j,kBottom,bi,bj) = GGL90TKE(i,j,kBottom,bi,bj)
599       I     myThid )       &                              - GGL90TKEbottom*c3d(i,j,kBottom)
600  C           c3d(i,j,kBottom) = 0. _d 0
601  C     now update TKE          ENDDO
602  C             ENDDO
603        DO K=1,Nr        ENDIF
604         DO J=jMin,jMax  
605          DO I=iMin,iMax  C     solve tri-diagonal system
606          CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
607         I                        a3d, b3d, c3d,
608         U                        GGL90TKE(1-OLx,1-OLy,1,bi,bj),
609         O                        errCode,
610         I                        bi, bj, myThid )
611    
612          DO k=1,Nr
613           DO j=jMin,jMax
614            DO i=iMin,iMax
615  C     impose minimum TKE to avoid numerical undershoots below zero  C     impose minimum TKE to avoid numerical undershoots below zero
616           GGL90TKE(I,J,K,bi,bj) = MAX( gTKE(I,J,K), GGL90TKEmin )           GGL90TKE(i,j,k,bi,bj) = maskC(i,j,k,bi,bj)
617       &        * maskC(I,J,K,bi,bj)       &                  *MAX( GGL90TKE(i,j,k,bi,bj), GGL90TKEmin )
618  C          ENDDO
619           ENDDO
620          ENDDO
621    
622  C     end of time step  C     end of time step
623  C  C     ===============================
624    
625          DO k=2,Nr
626           DO j=1,sNy
627            DO i=1,sNx
628    #ifdef ALLOW_GGL90_SMOOTH
629             tmpVisc=
630         &  (
631         &   p4 *  GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
632         &  +p8 *( GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj)
633         &       + GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj)
634         &       + GGL90visctmp(i+1,j  ,k) * mskCor(i+1,j  ,bi,bj)
635         &       + GGL90visctmp(i  ,j+1,k) * mskCor(i  ,j+1,bi,bj))
636         &  +p16*( GGL90visctmp(i+1,j+1,k) * mskCor(i+1,j+1,bi,bj)
637         &       + GGL90visctmp(i+1,j-1,k) * mskCor(i+1,j-1,bi,bj)
638         &       + GGL90visctmp(i-1,j+1,k) * mskCor(i-1,j+1,bi,bj)
639         &       + GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj))
640         &  )
641         & /(p4
642         &  +p8 *(       maskC(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)
643         &       +       maskC(i  ,j-1,k,bi,bj) * mskCor(i  ,j-1,bi,bj)
644         &       +       maskC(i+1,j  ,k,bi,bj) * mskCor(i+1,j  ,bi,bj)
645         &       +       maskC(i  ,j+1,k,bi,bj) * mskCor(i  ,j+1,bi,bj))
646         &  +p16*(       maskC(i+1,j+1,k,bi,bj) * mskCor(i+1,j+1,bi,bj)
647         &       +       maskC(i+1,j-1,k,bi,bj) * mskCor(i+1,j-1,bi,bj)
648         &       +       maskC(i-1,j+1,k,bi,bj) * mskCor(i-1,j+1,bi,bj)
649         &       +       maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj))
650         &  )*maskC(i,j,k,bi,bj)*mskCor(i,j,bi,bj)
651    #else
652             tmpVisc = GGL90visctmp(i,j,k)
653    #endif
654             tmpVisc = MIN(tmpVisc/TKEPrandtlNumber(i,j,k),GGL90diffMax)
655             GGL90diffKr(i,j,k,bi,bj)= MAX( tmpVisc , diffKrNrS(k) )
656          ENDDO          ENDDO
657         ENDDO         ENDDO
658        ENDDO            ENDDO
659  C  
660  C     compute viscosity coefficients        DO k=2,Nr
661  C             DO j=1,sNy
662        DO K=2,Nr          DO i=1,sNx+1
663         DO J=jMin,jMax  #ifdef ALLOW_GGL90_SMOOTH
664          DO I=iMin,iMax          tmpVisc =
665  C     Eq. (11), (18) and (21)       & (
666           KappaM = GGL90ck*GGL90mixingLength(I,J,K)*       &   p4 *(GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
667       &                  SQRT( GGL90TKE(I,J,K,bi,bj) )       &       +GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj))
668           KappaH = KappaM/TKEPrandtlNumber(I,J,K)       &  +p8 *(GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj)
669  C     Set a minium (= background) value       &       +GGL90visctmp(i-1,j+1,k) * mskCor(i-1,j+1,bi,bj)
670           KappaM = MAX(KappaM,viscAr)       &       +GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj)
671           KappaH = MAX(KappaH,diffKrT)       &       +GGL90visctmp(i  ,j+1,k) * mskCor(i  ,j+1,bi,bj))
672  C     Set a maximum and mask land point       &  )
673           GGL90viscAr(I,J,K,bi,bj) = MIN(KappaM,GGL90viscMax)       & /(p4 * 2. _d 0
674       &        * maskC(I,J,K,bi,bj)       &  +p8 *(      maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj)
675           GGL90diffKr(I,J,K,bi,bj) = MIN(KappaH,GGL90diffMax)       &       +      maskC(i-1,j+1,k,bi,bj) * mskCor(i-1,j+1,bi,bj)
676       &        * maskC(I,J,K,bi,bj)       &       +      maskC(i  ,j-1,k,bi,bj) * mskCor(i  ,j-1,bi,bj)
677          ENDDO       &       +      maskC(i  ,j+1,k,bi,bj) * mskCor(i  ,j+1,bi,bj))
678         ENDDO       &  )
679  C     end third k-loop       &  *maskC(i  ,j,k,bi,bj)*mskCor(i  ,j,bi,bj)
680        ENDDO           &  *maskC(i-1,j,k,bi,bj)*mskCor(i-1,j,bi,bj)
681    #else
682            tmpVisc = _maskW(i,j,k,bi,bj) *
683         &                   (.5 _d 0*(GGL90visctmp(i,j,k)
684         &                            +GGL90visctmp(i-1,j,k))
685         &                   )
686    #endif
687            tmpVisc = MIN( tmpVisc , GGL90viscMax )
688            GGL90viscArU(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
689            ENDDO
690           ENDDO
691          ENDDO
692    
693          DO k=2,Nr
694           DO j=1,sNy+1
695            DO i=1,sNx
696    #ifdef ALLOW_GGL90_SMOOTH
697            tmpVisc =
698         & (
699         &   p4 *(GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
700         &       +GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj))
701         &  +p8 *(GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj)
702         &       +GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj)
703         &       +GGL90visctmp(i+1,j  ,k) * mskCor(i+1,j  ,bi,bj)
704         &       +GGL90visctmp(i+1,j-1,k) * mskCor(i+1,j-1,bi,bj))
705         &  )
706         & /(p4 * 2. _d 0
707         &  +p8 *(      maskC(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)
708         &       +      maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj)
709         &       +      maskC(i+1,j  ,k,bi,bj) * mskCor(i+1,j  ,bi,bj)
710         &       +      maskC(i+1,j-1,k,bi,bj) * mskCor(i+1,j-1,bi,bj))
711         &  )
712         &   *maskC(i,j  ,k,bi,bj)*mskCor(i,j  ,bi,bj)
713         &   *maskC(i,j-1,k,bi,bj)*mskCor(i,j-1,bi,bj)
714    #else
715            tmpVisc = _maskS(i,j,k,bi,bj) *
716         &                   (.5 _d 0*(GGL90visctmp(i,j,k)
717         &                            +GGL90visctmp(i,j-1,k))
718         &                   )
719    
720    #endif
721            tmpVisc = MIN( tmpVisc , GGL90viscMax )
722            GGL90viscArV(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
723            ENDDO
724           ENDDO
725          ENDDO
726    
727    #ifdef ALLOW_DIAGNOSTICS
728          IF ( useDiagnostics ) THEN
729            CALL DIAGNOSTICS_FILL( GGL90TKE   ,'GGL90TKE',
730         &                         0,Nr, 1, bi, bj, myThid )
731            CALL DIAGNOSTICS_FILL( GGL90viscArU,'GGL90ArU',
732         &                         0,Nr, 1, bi, bj, myThid )
733            CALL DIAGNOSTICS_FILL( GGL90viscArV,'GGL90ArV',
734         &                         0,Nr, 1, bi, bj, myThid )
735            CALL DIAGNOSTICS_FILL( GGL90diffKr,'GGL90Kr ',
736         &                         0,Nr, 1, bi, bj, myThid )
737            CALL DIAGNOSTICS_FILL( TKEPrandtlNumber ,'GGL90Prl',
738         &                         0,Nr, 2, bi, bj, myThid )
739            CALL DIAGNOSTICS_FILL( GGL90mixingLength,'GGL90Lmx',
740         &                         0,Nr, 2, bi, bj, myThid )
741    
742            kp1 = MIN(Nr,kSurf+1)
743            DO j=jMin,jMax
744             DO i=iMin,iMax
745    C     diagnose surface flux of TKE
746              surf_flx_tke(i,j) =(GGL90TKE(i,j,kSurf,bi,bj)-
747         &                        GGL90TKE(i,j,kp1,bi,bj))
748         &        *recip_drF(kSurf)*recip_hFacC(i,j,kSurf,bi,bj)
749         &        *KappaE(i,j,kp1)
750             ENDDO
751            ENDDO
752            CALL DIAGNOSTICS_FILL( surf_flx_tke,'GGL90flx',
753         &                         0, 1, 2, bi, bj, myThid )
754    
755            k=kSurf
756            DO j=jMin,jMax
757             DO i=iMin,iMax
758    C     diagnose work done by the wind
759              surf_flx_tke(i,j) =
760         &      halfRL*( surfaceForcingU(i,  j,bi,bj)*uVel(i  ,j,k,bi,bj)
761         &              +surfaceForcingU(i+1,j,bi,bj)*uVel(i+1,j,k,bi,bj))
762         &    + halfRL*( surfaceForcingV(i,j,  bi,bj)*vVel(i,j  ,k,bi,bj)
763         &              +surfaceForcingV(i,j+1,bi,bj)*vVel(i,j+1,k,bi,bj))
764             ENDDO
765            ENDDO
766            CALL DIAGNOSTICS_FILL( surf_flx_tke,'GGL90tau',
767         &                         0, 1, 2, bi, bj, myThid )
768    
769          ENDIF
770    #endif /* ALLOW_DIAGNOSTICS */
771    
772  #endif /* ALLOW_GGL90 */  #endif /* ALLOW_GGL90 */
773    
774        RETURN        RETURN
775        END        END
   
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.31
  ViewVC Help
Powered by ViewVC 1.1.22