/[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.25 by mlosch, Mon Nov 5 13:11:11 2012 UTC revision 1.34 by jmc, Thu Jan 14 17:50:25 2016 UTC
# Line 10  C !INTERFACE: ========================== Line 10  C !INTERFACE: ==========================
10        SUBROUTINE GGL90_CALC(        SUBROUTINE GGL90_CALC(
11       I                 bi, bj, sigmaR, myTime, myIter, 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                                    |
# Line 36  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
# Line 52  C     myThid :: My Thread Id number Line 51  C     myThid :: My Thread Id number
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 61  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 79  C     TKEPrandtlNumber :: here, an empir Line 78  C     TKEPrandtlNumber :: here, an empir
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)
# Line 97  c     _RL     SQRTTKE Line 100  c     _RL     SQRTTKE
100        _RL     KappaE           (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     KappaE           (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
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     a3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     a3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
114        _RL     b3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     b3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
115        _RL     c3d(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
 C     hFac     :: fractional thickness of W-cell  
118  C     xA, yA   :: area of lateral faces  C     xA, yA   :: area of lateral faces
119  C     dfx, dfy :: diffusive flux across lateral faces  C     dfx, dfy :: diffusive flux across lateral faces
120  C     gTKE     :: right hand side of diffusion equation  C     gTKE     :: right hand side of diffusion equation
       _RL     hFac  
121        _RL     xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
122        _RL     yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
123        _RL     dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 116  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 133  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
# Line 146  C     Initialize local fields Line 181  C     Initialize local fields
181           b3d(i,j,k) = 1. _d 0           b3d(i,j,k) = 1. _d 0
182           c3d(i,j,k) = 0. _d 0           c3d(i,j,k) = 0. _d 0
183  #endif  #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            KappaM(i,j)        = 0. _d 0
192            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
 c      km1 = k-1  
 c      kp1 = MIN(Nr,k+1)  
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) )
# Line 169  c      kp1 = MIN(Nr,k+1) Line 217  c      kp1 = MIN(Nr,k+1)
217  C     buoyancy frequency  C     buoyancy frequency
218           Nsquare(i,j,k) = gravity*gravitySign*recip_rhoConst           Nsquare(i,j,k) = gravity*gravitySign*recip_rhoConst
219       &                  * sigmaR(i,j,k)       &                  * 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 198  C- ensure mixing between first and secon Line 235  C- ensure mixing between first and secon
235         ENDDO         ENDDO
236        ENDIF        ENDIF
237    
238  C- Impose upper and lower bound for mixing length  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    
242         DO k=2,Nr         DO k=2,Nr
# Line 320  c         MaxLength=MAX(MaxLength,20. _d Line 358  c         MaxLength=MAX(MaxLength,20. _d
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    
# Line 379  C     ... across y-faces Line 408  C     ... across y-faces
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            hFac = min(.5 _d 0,_hFacC(i,j,k-1,bi,bj) ) +            gTKE(i,j) = -recip_drC(k)*recip_rA(i,j,bi,bj)
412       &          min(.5 _d 0,_hFacC(i,j,k  ,bi,bj) )       &         *recip_hFacI(i,j,k)
           gTKE(i,j) = 0.0  
           if ( hFac .ne. 0.0 )  
      &      gTKE(i,j) = -recip_drC(k)*recip_rA(i,j,bi,bj)/hFac  
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 426  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 482  C--   Upper diagonal Line 576  C--   Upper diagonal
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            c3d(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            b3d(i,j,k) = 1. _d 0 - c3d(i,j,k) - a3d(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 524  C     Dirichlet surface boundary conditi Line 642  C     Dirichlet surface boundary conditi
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       &               - a3d(i,j,kp1)*GGL90TKE(i,j,kSurf,bi,bj)       &               - a3d(i,j,kp1)*GGL90TKE(i,j,kSurf,bi,bj)
644          a3d(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*c3d(i,j,kBottom)  
         c3d(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                        a3d, b3d, c3d,       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 556  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 591  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 624  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.25  
changed lines
  Added in v.1.34
  ViewVC Help
Powered by ViewVC 1.1.22