/[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.34 by jmc, Thu Jan 14 17:50:25 2016 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"
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     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     hFac/hFacI       :: fractional thickness of W-cell
63  C     verticalShear    - (squared) vertical shear of horizontal velocity  C     explDissFac      :: explicit Dissipation Factor (in [0-1])
64  C     Nsquare          - squared buoyancy freqency  C     implDissFac      :: implicit Dissipation Factor (in [0-1])
65  C     RiNumber         - local Richardson number  C     uStarSquare      :: square of friction velocity
66  C     KappaM           - (local) viscosity parameter (eq.10)  C     verticalShear    :: (squared) vertical shear of horizontal velocity
67  C     KappaH           - (local) diffusivity parameter for temperature (eq.11)  C     Nsquare          :: squared buoyancy freqency
68  C     KappaE           - (local) diffusivity parameter for TKE (eq.15)  C     RiNumber         :: local Richardson number
69  C     buoyFreq         - buoyancy freqency  C     KappaM           :: (local) viscosity parameter (eq.10)
70  C     TKEdissipation   - dissipation of TKE  C     KappaH           :: (local) diffusivity parameter for temperature (eq.11)
71  C     GGL90mixingLength- mixing length of scheme following Banke+Delecuse  C     KappaE           :: (local) diffusivity parameter for TKE (eq.15)
72  C     totalDepth       - thickness of water column (inverse of recip_Rcol)  C     TKEdissipation   :: dissipation of TKE
73  C     TKEPrandtlNumber - here, an empirical function of the Richardson number  C     GGL90mixingLength:: mixing length of scheme following Banke+Delecuse
74  C     rhoK, rhoKm1     - density at layer K and Km1 (relative to K)  C         rMixingLength:: inverse of mixing length
75  C     gTKE             - right hand side of implicit equation  C     totalDepth       :: thickness of water column (inverse of recip_Rcol)
76    C     TKEPrandtlNumber :: here, an empirical function of the Richardson number
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     ab15, ab05        _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     Nsquare        _RL     KappaH
84    c     _RL     Nsquare
85          _RL     Nsquare(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
86        _RL     deltaTggl90        _RL     deltaTggl90
87        _RL     SQRTTKE  c     _RL     SQRTTKE
88          _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
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)
99          _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     gTKE             (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     GGL90visctmp     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
103  C     tri-diagonal matrix  #ifdef ALLOW_DIAGNOSTICS
104        _RL     a(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     surf_flx_tke     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105        _RL     b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  #endif /* ALLOW_DIAGNOSTICS */
106        _RL     c(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  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
113          _RL     a3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
114          _RL     b3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
115          _RL     c3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
116          INTEGER errCode
117    #ifdef ALLOW_GGL90_HORIZDIFF
118    C     xA, yA   :: area of lateral faces
119    C     dfx, dfy :: diffusive flux across lateral faces
120    C     gTKE     :: right hand side of diffusion equation
121          _RL     xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
122          _RL     yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
123          _RL     dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
124          _RL     dfy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
125          _RL    gTKE(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
126    #endif /* ALLOW_GGL90_HORIZDIFF */
127    #ifdef ALLOW_GGL90_SMOOTH
128          _RL p4, p8, p16
129    #endif
130  CEOP  CEOP
131        iMin = 2-OLx  
132        iMax = sNx+OLx-1        PARAMETER( iMin = 2-OLx, iMax = sNx+OLx-1 )
133        jMin = 2-OLy        PARAMETER( jMin = 2-OLy, jMax = sNy+OLy-1 )
134        jMax = sNy+OLy-1  #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 = deltaTtracer        deltaTggl90 = dTtracerLev(1)
142  C      
143        kSurf = 1        kSurf = 1
144  C     implicit timestepping weights for dissipation  C     explicit/implicit timestepping weights for dissipation
145        ab15 =  1.5 _d 0        explDissFac = 0. _d 0
146        ab05 = -0.5 _d 0        implDissFac = 1. _d 0 - explDissFac
147        ab15 =  1. _d 0  
148        ab05 =  0. _d 0  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           gTKE(I,J,K)              = 0. _d 0           rMixingLength(i,j,k)     = 0. _d 0
173           KappaE(I,J,K)            = 0. _d 0           mxLength_Dn(i,j,k)       = 0. _d 0
174           TKEPrandtlNumber(I,J,K)  = 0. _d 0           GGL90visctmp(i,j,k)      = 0. _d 0
175           GGL90mixingLength(I,J,K) = 0. _d 0           KappaE(i,j,k)            = 0. _d 0
176          ENDDO           TKEPrandtlNumber(i,j,k)  = 1. _d 0
177         ENDDO               GGL90mixingLength(i,j,k) = GGL90mixingLengthMin
178        ENDDO           GGL90visctmp(i,j,k)      = 0. _d 0
179        DO J=1-Oly,sNy+Oly  #ifndef SOLVE_DIAGONAL_LOWMEMORY
180         DO I=1-Olx,sNx+Olx           a3d(i,j,k) = 0. _d 0
181          rhoK    (I,J) = 0. _d 0           b3d(i,j,k) = 1. _d 0
182          rhoKm1  (I,J) = 0. _d 0           c3d(i,j,k) = 0. _d 0
183          totalDepth(I,J)   = 0. _d 0  #endif
184          IF ( recip_Rcol(I,J,bi,bj) .NE. 0. )           Nsquare(i,j,k) = 0. _d 0
185       &       totalDepth(I,J) = 1./recip_Rcol(I,J,bi,bj)           SQRTTKE(i,j,k) = 0. _d 0
186            ENDDO
187           ENDDO
188          ENDDO
189          DO j=1-OLy,sNy+OLy
190           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)
194            rMixingLength(i,j,1) = 0. _d 0
195            mxLength_Dn(i,j,1) = GGL90mixingLengthMin
196            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        DO K = 2, Nr        IF ( useIDEMIX) CALL GGL90_IDEMIX(
209         Km1 = K-1       &   bi, bj, hFacI, recip_hFacI, sigmaR, myTime, myIter, myThid )
210         Kp1 = MIN(Nr,K+1)  #endif /* ALLOW_GGL90_IDEMIX */
211         CALL FIND_RHO(  
212       I      bi, bj, iMin, iMax, jMin, jMax, Km1, K,        DO k = 2, Nr
213       I      theta, salt,         DO j=jMin,jMax
214       O      rhoKm1,          DO i=iMin,iMax
215       I      myThid )           SQRTTKE(i,j,k)=SQRT( GGL90TKE(i,j,k,bi,bj) )
216         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  
217  C     buoyancy frequency  C     buoyancy frequency
218  C           Nsquare(i,j,k) = gravity*gravitySign*recip_rhoConst
219           Nsquare = - gravity*recip_rhoConst*recip_drC(K)       &                  * sigmaR(i,j,k)
220       &        * ( rhoKm1(I,J) - rhoK(I,J) )*maskC(I,J,K,bi,bj)  C     vertical shear term (dU/dz)^2+(dV/dz)^2 is computed later
221  C     vertical shear term (dU/dz)^2+(dV/dz)^2  C     to save some memory
          tempu= .5*(  uVel(I,J,Km1,bi,bj)+uVel(I+1,J,Km1,bi,bj)  
      &             - (uVel(I,J,K  ,bi,bj)+uVel(I+1,J,K  ,bi,bj)) )  
      &        *recip_drC(K)  
          tempv= .5*(  vVel(I,J,Km1,bi,bj)+vVel(I,J+1,Km1,bi,bj)  
      &             - (vVel(I,J,K  ,bi,bj)+vVel(I,J+1,K  ,bi,bj)) )  
      &        *recip_drC(K)  
          verticalShear = tempU*tempU + tempV*tempV  
          RiNumber   = MAX(Nsquare,0. _d 0)/(verticalShear+GGL90eps)  
 C     compute Prandtl number (always greater than 0)  
          prTemp = 1. _d 0  
          IF ( RiNumber .GE. 0.2 ) prTemp = 5.0 * RiNumber  
          TKEPrandtlNumber(I,J,K) = MIN(10.,prTemp)  
222  C     mixing length  C     mixing length
223           GGL90mixingLength(I,J,K) =           GGL90mixingLength(i,j,k) = SQRTTWO *
224       &        SQRTTKE/SQRT( MAX(Nsquare,GGL90eps) )       &        SQRTTKE(i,j,k)/SQRT( MAX(Nsquare(i,j,k),GGL90eps) )
225  C     impose upper bound for mixing length (total depth)          ENDDO
226           GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),         ENDDO
227       &        totalDepth(I,J))        ENDDO
228  C     impose minimum mixing length (to avoid division by zero)  
229           GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),  C- ensure mixing between first and second level
230       &        GGL90mixingLengthMin)        IF (mxlSurfFlag) THEN
231  C     viscosity of last timestep         DO j=jMin,jMax
232           KappaM = GGL90ck*GGL90mixingLength(I,J,K)*SQRTTKE          DO i=iMin,iMax
233           KappaE(I,J,K) = KappaM*GGL90alpha           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
241    
242           DO k=2,Nr
243            DO j=jMin,jMax
244             DO i=iMin,iMax
245              MaxLength=totalDepth(i,j)
246              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
247         &                                   MaxLength)
248             ENDDO
249            ENDDO
250           ENDDO
251    
252           DO k=2,Nr
253            DO j=jMin,jMax
254             DO i=iMin,iMax
255              GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
256         &                                   GGL90mixingLengthMin)
257              rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
258             ENDDO
259            ENDDO
260           ENDDO
261    
262          ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN
263    
264           DO k=2,Nr
265            DO j=jMin,jMax
266             DO i=iMin,iMax
267              MaxLength=MIN(Ro_surf(i,j,bi,bj)-rF(k),rF(k)-R_low(i,j,bi,bj))
268    c         MaxLength=MAX(MaxLength,20. _d 0)
269              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
270         &                                   MaxLength)
271             ENDDO
272            ENDDO
273           ENDDO
274    
275           DO k=2,Nr
276            DO j=jMin,jMax
277             DO i=iMin,iMax
278              GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
279         &                                   GGL90mixingLengthMin)
280              rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
281             ENDDO
282            ENDDO
283           ENDDO
284    
285          ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN
286    
287           DO k=2,Nr
288            DO j=jMin,jMax
289             DO i=iMin,iMax
290              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
291         &        GGL90mixingLength(i,j,k-1)+drF(k-1))
292             ENDDO
293            ENDDO
294           ENDDO
295           DO j=jMin,jMax
296            DO i=iMin,iMax
297              GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
298         &       GGL90mixingLengthMin+drF(Nr))
299            ENDDO
300           ENDDO
301           DO k=Nr-1,2,-1
302            DO j=jMin,jMax
303             DO i=iMin,iMax
304              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
305         &        GGL90mixingLength(i,j,k+1)+drF(k))
306             ENDDO
307            ENDDO
308           ENDDO
309    
310           DO k=2,Nr
311            DO j=jMin,jMax
312             DO i=iMin,iMax
313              GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
314         &                                   GGL90mixingLengthMin)
315              rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
316             ENDDO
317            ENDDO
318           ENDDO
319    
320          ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN
321    
322           DO k=2,Nr
323            DO j=jMin,jMax
324             DO i=iMin,iMax
325              mxLength_Dn(i,j,k) = MIN(GGL90mixingLength(i,j,k),
326         &        mxLength_Dn(i,j,k-1)+drF(k-1))
327             ENDDO
328            ENDDO
329           ENDDO
330           DO j=jMin,jMax
331            DO i=iMin,iMax
332              GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
333         &       GGL90mixingLengthMin+drF(Nr))
334            ENDDO
335           ENDDO
336           DO k=Nr-1,2,-1
337            DO j=jMin,jMax
338             DO i=iMin,iMax
339              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
340         &        GGL90mixingLength(i,j,k+1)+drF(k))
341             ENDDO
342            ENDDO
343           ENDDO
344    
345           DO k=2,Nr
346            DO j=jMin,jMax
347             DO i=iMin,iMax
348              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
349         &                                  mxLength_Dn(i,j,k))
350              tmpmlx = SQRT( GGL90mixingLength(i,j,k)*mxLength_Dn(i,j,k) )
351              tmpmlx = MAX( tmpmlx, GGL90mixingLengthMin)
352              rMixingLength(i,j,k) = 1. _d 0 / tmpmlx
353             ENDDO
354            ENDDO
355           ENDDO
356    
357          ELSE
358           STOP 'GGL90_CALC: Wrong mxlMaxFlag (mixing length limit)'
359          ENDIF
360    
361    C     start "proper" k-loop (the code above was moved out and up to
362    C     implemement various mixing parameters efficiently)
363          DO k=2,Nr
364           km1 = k-1
365    
366    #ifdef ALLOW_GGL90_HORIZDIFF
367          IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
368    C     horizontal diffusion of TKE (requires an exchange in
369    C      do_fields_blocking_exchanges)
370    C     common factors
371            DO j=1-OLy,sNy+OLy
372             DO i=1-OLx,sNx+OLx
373              xA(i,j) = _dyG(i,j,bi,bj)*drC(k)*
374         &                 (min(.5 _d 0,_hFacW(i,j,k-1,bi,bj) ) +
375         &                  min(.5 _d 0,_hFacW(i,j,k  ,bi,bj) ) )
376              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
380            ENDDO
381    C     Compute diffusive fluxes
382    C     ... across x-faces
383            DO j=1-OLy,sNy+OLy
384             dfx(1-OLx,j)=0. _d 0
385             DO i=1-OLx+1,sNx+OLx
386              dfx(i,j) = -GGL90diffTKEh*xA(i,j)
387         &      *_recip_dxC(i,j,bi,bj)
388         &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i-1,j,k,bi,bj))
389    #ifdef ISOTROPIC_COS_SCALING
390         &      *CosFacU(j,bi,bj)
391    #endif /* ISOTROPIC_COS_SCALING */
392             ENDDO
393            ENDDO
394    C     ... across y-faces
395            DO i=1-OLx,sNx+OLx
396             dfy(i,1-OLy)=0. _d 0
397            ENDDO
398            DO j=1-OLy+1,sNy+OLy
399             DO i=1-OLx,sNx+OLx
400              dfy(i,j) = -GGL90diffTKEh*yA(i,j)
401         &      *_recip_dyC(i,j,bi,bj)
402         &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i,j-1,k,bi,bj))
403    #ifdef ISOTROPIC_COS_SCALING
404         &      *CosFacV(j,bi,bj)
405    #endif /* ISOTROPIC_COS_SCALING */
406             ENDDO
407            ENDDO
408    C     Compute divergence of fluxes
409            DO j=1-OLy,sNy+OLy-1
410             DO i=1-OLx,sNx+OLx-1
411              gTKE(i,j) = -recip_drC(k)*recip_rA(i,j,bi,bj)
412         &         *recip_hFacI(i,j,k)
413         &         *((dfx(i+1,j)-dfx(i,j))
414         &         + (dfy(i,j+1)-dfy(i,j)) )
415             ENDDO
416            ENDDO
417    C     end if GGL90diffTKEh .eq. 0.
418           ENDIF
419    #endif /* ALLOW_GGL90_HORIZDIFF */
420    
421    C     viscosity and diffusivity
422           DO j=jMin,jMax
423            DO i=iMin,iMax
424             KappaM(i,j) = GGL90ck*GGL90mixingLength(i,j,k)*SQRTTKE(i,j,k)
425             GGL90visctmp(i,j,k) = MAX(KappaM(i,j),diffKrNrS(k))
426         &                            * maskC(i,j,k,bi,bj)
427    C        note: storing GGL90visctmp like this, and using it later to compute
428    C              GGL9rdiffKr etc. is robust in case of smoothing (e.g. see OPA)
429             KappaM(i,j) = MAX(KappaM(i,j),viscArNr(k)) * maskC(i,j,k,bi,bj)
430            ENDDO
431           ENDDO
432    
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)
464    #ifdef ALLOW_GGL90_IDEMIX
465           IF ( useIDEMIX ) THEN
466            DO j=jMin,jMax
467             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           DO j=jMin,jMax
497            DO i=iMin,iMax
498    C        diffusivity
499             KappaH = KappaM(i,j)/TKEPrandtlNumber(i,j,k)
500             KappaE(i,j,k) = GGL90alpha * KappaM(i,j) * maskC(i,j,k,bi,bj)
501    
502  C     dissipation term  C     dissipation term
503           TKEdissipation = ab05*GGL90ceps           TKEdissipation = explDissFac*GGL90ceps
504       &        *SQRTTKE/GGL90mixingLength(I,J,K)       &        *SQRTTKE(i,j,k)*rMixingLength(i,j,k)
505       &        *GGL90TKE(I,J,K,bi,bj)             &        *GGL90TKE(i,j,k,bi,bj)
506  C     sum up contributions to form the right hand side  C     partial update with sum of explicit contributions
507           gTKE(I,J,K) = 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       &        - KappaM*Nsquare/TKEPrandtlNumber(I,J,K)       &        - KappaH*Nsquare(i,j,k)
511       &        - TKEdissipation       &        - TKEdissipation
512       &        )       &        )
513          ENDDO            ENDDO
514         ENDDO         ENDDO
515        ENDDO      
516  C      #ifdef ALLOW_GGL90_IDEMIX
517  C     Implicit time step to update TKE for k=1,Nr; TKE(Nr+1)=0 by default         IF ( useIDEMIX ) THEN
518  C      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
531           IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
532    C--    Add horiz. diffusion tendency
533            DO j=jMin,jMax
534             DO i=iMin,iMax
535              GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
536         &                          + gTKE(i,j)*deltaTggl90
537             ENDDO
538            ENDDO
539           ENDIF
540    #endif /* ALLOW_GGL90_HORIZDIFF */
541    
542    C--   end of k loop
543          ENDDO
544    
545    C     ============================================
546    C     Implicit time step to update TKE for k=1,Nr;
547    C     TKE(Nr+1)=0 by default
548    C     ============================================
549  C     set up matrix  C     set up matrix
550  C--   Old aLower  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
557         km1=MAX(1,k-1)         km1=MAX(2,k-1)
558         DO j=jMin,jMax         DO j=jMin,jMax
559          DO i=iMin,iMax          DO i=iMin,iMax
560           a(i,j,k) = -deltaTggl90  C-    We keep recip_hFacC in the diffusive flux calculation,
561       &        *recip_drF(km1)*recip_hFacI(i,j,k,bi,bj)  C-    but no hFacC in TKE volume control
562       &        *.5*(KappaE(i,j, k )+KappaE(i,j,km1))  C-    No need for maskC(k-1) with recip_hFacC(k-1)
563       &        *recip_drC(k)           a3d(i,j,k) = -deltaTggl90
564            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)
565         &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))
566         &        *recip_drC(k)*maskC(i,j,k,bi,bj)
567          ENDDO          ENDDO
568         ENDDO         ENDDO
569        ENDDO        ENDDO
570    C--   Upper diagonal
 C--   Old aUpper  
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
          c(i,j,Nr) = 0. _d 0  
574         ENDDO         ENDDO
575        ENDDO        ENDDO
 CML      DO k=1,Nr-1  
576        DO k=2,Nr        DO k=2,Nr
        kp1=min(Nr,k+1)  
577         DO j=jMin,jMax         DO j=jMin,jMax
578          DO i=iMin,iMax          DO i=iMin,iMax
579            c(i,j,k) = -deltaTggl90           kp1=MAX(1,MIN(klowC(i,j,bi,bj),k+1))
580       &        *recip_drF( k )*recip_hFacI(i,j,k,bi,bj)  C-    We keep recip_hFacC in the diffusive flux calculation,
581       &               *.5*(KappaE(i,j,k)+KappaE(i,j,kp1))  C-    but no hFacC in TKE volume control
582       &        *recip_drC(k)  C-    No need for maskC(k) with recip_hFacC(k)
583  C          IF (recip_hFacI(i,j,kp1,bi,bj).EQ.0.) c(i,j,k)=0.           c3d(i,j,k) = -deltaTggl90
584         &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)
585         &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))
586         &        *recip_drC(k)*maskC(i,j,k-1,bi,bj)
587          ENDDO          ENDDO
588         ENDDO         ENDDO
589        ENDDO        ENDDO
590    
591  C--   Old aCenter  #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
615        DO k=1,Nr        DO k=1,Nr
616           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       &        + ab15*deltaTggl90*GGL90ceps*SQRT(GGL90TKE(I,J,K,bi,bj))       &        + implDissFac*deltaTggl90*GGL90ceps*SQRTTKE(i,j,k)
621       &        /GGL90mixingLength(I,J,K)       &        * rMixingLength(i,j,k)
622           ENDDO       &        * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)
623            ENDDO
624         ENDDO         ENDDO
625        ENDDO        ENDDO
626  C     end set up matrix  C     end set up matrix
627    
 C  
628  C     Apply boundary condition  C     Apply boundary condition
629  C            kp1 = MIN(Nr,kSurf+1)
630        DO J=jMin,jMax        DO j=jMin,jMax
631         DO I=iMin,iMax         DO i=iMin,iMax
632  C     estimate friction velocity uStar from surface forcing  C     estimate friction velocity uStar from surface forcing
633          uStarSquare = SQRT(          uStarSquare = SQRT(
634       &         ( .5*( surfaceForcingU(I,  J,  bi,bj)       &    ( .5 _d 0*( surfaceForcingU(i,  j,  bi,bj)
635       &              + surfaceForcingU(I+1,J,  bi,bj) ) )**2       &              + surfaceForcingU(i+1,j,  bi,bj) ) )**2
636       &       + ( .5*( surfaceForcingV(I,  J,  bi,bj)       &  + ( .5 _d 0*( surfaceForcingV(i,  j,  bi,bj)
637       &              + surfaceForcingV(I,  J+1,bi,bj) ) )**2       &              + surfaceForcingV(i,  j+1,bi,bj) ) )**2
638       &                     )       &                     )
639  C     Dirichlet surface boundary condition for TKE  C     Dirichlet surface boundary condition for TKE
640          gTKE(I,J,kSurf) = MAX(GGL90TKEmin,GGL90m2*uStarSquare)          GGL90TKE(i,j,kSurf,bi,bj) = maskC(i,j,kSurf,bi,bj)
641       &                     *maskC(I,J,kSurf,bi,bj)       &           *MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)
642  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom          GGL90TKE(i,j,kp1,bi,bj) = GGL90TKE(i,j,kp1,bi,bj)
643          kBottom   = MIN(MAX(kLowC(I,J,bi,bj),1),Nr)       &               - a3d(i,j,kp1)*GGL90TKE(i,j,kSurf,bi,bj)
644          gTKE(I,J,kBottom) = gTKE(I,J,kBottom)          a3d(i,j,kp1) = 0. _d 0
645       &       - GGL90TKEbottom*c(I,J,kBottom)         ENDDO
646         ENDDO        ENDDO
647        ENDDO      
648  C        IF (GGL90_dirichlet) THEN
649  C     solve tri-diagonal system, and store solution on gTKE (previously rhs)  C      Dirichlet bottom boundary condition for TKE = GGL90TKEbottom
650  C         DO j=jMin,jMax
651        CALL GGL90_SOLVE( bi, bj, iMin, iMax, jMin, jMax,          DO i=iMin,iMax
652       I     a, b, c,           kBottom   = MAX(kLowC(i,j,bi,bj),1)
653       U     gTKE,           GGL90TKE(i,j,kBottom,bi,bj) = GGL90TKE(i,j,kBottom,bi,bj)
654       I     myThid )       &                              - GGL90TKEbottom*c3d(i,j,kBottom)
655  C           c3d(i,j,kBottom) = 0. _d 0
656  C     now update TKE          ENDDO
657  C             ENDDO
658        DO K=1,Nr        ENDIF
659         DO J=jMin,jMax  
660          DO I=iMin,iMax  C     solve tri-diagonal system
661          CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
662         I                        a3d, b3d, c3d,
663         U                        GGL90TKE(1-OLx,1-OLy,1,bi,bj),
664         O                        errCode,
665         I                        bi, bj, myThid )
666    
667          DO k=1,Nr
668           DO j=jMin,jMax
669            DO i=iMin,iMax
670  C     impose minimum TKE to avoid numerical undershoots below zero  C     impose minimum TKE to avoid numerical undershoots below zero
671           GGL90TKE(I,J,K,bi,bj) = MAX( gTKE(I,J,K), GGL90TKEmin )           GGL90TKE(i,j,k,bi,bj) = maskC(i,j,k,bi,bj)
672       &        * maskC(I,J,K,bi,bj)       &                  *MAX( GGL90TKE(i,j,k,bi,bj), GGL90TKEmin )
673  C          ENDDO
674           ENDDO
675          ENDDO
676    
677  C     end of time step  C     end of time step
678  C  C     ===============================
679    
680          DO k=2,Nr
681           DO j=1,sNy
682            DO i=1,sNx
683    #ifdef ALLOW_GGL90_SMOOTH
684             tmpVisc = (
685         &     p4 *    GGL90visctmp(i  ,j  ,k)*mskCor(i  ,j  ,bi,bj)
686         &    +p8 *( ( GGL90visctmp(i-1,j  ,k)*mskCor(i-1,j  ,bi,bj)
687         &           + GGL90visctmp(i+1,j  ,k)*mskCor(i+1,j  ,bi,bj) )
688         &         + ( GGL90visctmp(i  ,j-1,k)*mskCor(i  ,j-1,bi,bj)
689         &           + GGL90visctmp(i  ,j+1,k)*mskCor(i  ,j+1,bi,bj) ) )
690         &    +p16*( ( GGL90visctmp(i+1,j+1,k)*mskCor(i+1,j+1,bi,bj)
691         &           + 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)
693         &           + GGL90visctmp(i-1,j+1,k)*mskCor(i-1,j+1,bi,bj) ) )
694         &             )/(
695         &     p4
696         &    +p8 *( (  maskC(i-1,j  ,k,bi,bj)*mskCor(i-1,j  ,bi,bj)
697         &           +  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)
699         &           +  maskC(i  ,j+1,k,bi,bj)*mskCor(i  ,j+1,bi,bj) ) )
700         &    +p16*( (  maskC(i+1,j+1,k,bi,bj)* mskCor(i+1,j+1,bi,bj)
701         &           +  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)
703         &           +  maskC(i-1,j+1,k,bi,bj)*mskCor(i-1,j+1,bi,bj) ) )
704         &               )*maskC(i,j,k,bi,bj)*mskCor(i,j,bi,bj)
705    #else
706             tmpVisc = GGL90visctmp(i,j,k)
707    #endif
708             tmpVisc = MIN(tmpVisc/TKEPrandtlNumber(i,j,k),GGL90diffMax)
709             GGL90diffKr(i,j,k,bi,bj)= MAX( tmpVisc , diffKrNrS(k) )
710          ENDDO          ENDDO
711         ENDDO         ENDDO
712        ENDDO            ENDDO
713  C  
714  C     compute viscosity coefficients        DO k=2,Nr
715  C             DO j=1,sNy
716        DO K=2,Nr          DO i=1,sNx+1
717         DO J=jMin,jMax  #ifdef ALLOW_GGL90_SMOOTH
718          DO I=iMin,iMax           tmpVisc = (
719  C     Eq. (11), (18) and (21)       &     p4 *(   GGL90visctmp(i-1,j  ,k)*mskCor(i-1,j  ,bi,bj)
720           KappaM = GGL90ck*GGL90mixingLength(I,J,K)*       &           + GGL90visctmp(i  ,j  ,k)*mskCor(i  ,j  ,bi,bj) )
721       &                  SQRT( GGL90TKE(I,J,K,bi,bj) )       &    +p8 *( ( GGL90visctmp(i-1,j-1,k)*mskCor(i-1,j-1,bi,bj)
722           KappaH = KappaM/TKEPrandtlNumber(I,J,K)       &           + GGL90visctmp(i  ,j-1,k)*mskCor(i  ,j-1,bi,bj) )
723  C     Set a minium (= background) value       &         + ( GGL90visctmp(i-1,j+1,k)*mskCor(i-1,j+1,bi,bj)
724           KappaM = MAX(KappaM,viscAr)       &           + GGL90visctmp(i  ,j+1,k)*mskCor(i  ,j+1,bi,bj) ) )
725           KappaH = MAX(KappaH,diffKrT)       &             )/(
726  C     Set a maximum and mask land point       &     p4 * 2. _d 0
727           GGL90viscAr(I,J,K,bi,bj) = MIN(KappaM,GGL90viscMax)       &    +p8 *( (  maskC(i-1,j-1,k,bi,bj)*mskCor(i-1,j-1,bi,bj)
728       &        * maskC(I,J,K,bi,bj)       &           +  maskC(i  ,j-1,k,bi,bj)*mskCor(i  ,j-1,bi,bj) )
729           GGL90diffKr(I,J,K,bi,bj) = MIN(KappaH,GGL90diffMax)       &         + (  maskC(i-1,j+1,k,bi,bj)*mskCor(i-1,j+1,bi,bj)
730       &        * maskC(I,J,K,bi,bj)       &           +  maskC(i  ,j+1,k,bi,bj)*mskCor(i  ,j+1,bi,bj) ) )
731          ENDDO       &               )*maskC(i-1,j,k,bi,bj)*mskCor(i-1,j,bi,bj)
732         ENDDO       &                *maskC(i  ,j,k,bi,bj)*mskCor(i  ,j,bi,bj)
733  C     end third k-loop  #else
734        ENDDO               tmpVisc = _maskW(i,j,k,bi,bj) * halfRL
735         &          *( GGL90visctmp(i-1,j,k)
736         &           + GGL90visctmp(i,j,k) )
737    #endif
738             tmpVisc = MIN( tmpVisc , GGL90viscMax )
739             GGL90viscArU(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
740            ENDDO
741           ENDDO
742          ENDDO
743    
744          DO k=2,Nr
745           DO j=1,sNy+1
746            DO i=1,sNx
747    #ifdef ALLOW_GGL90_SMOOTH
748             tmpVisc = (
749         &     p4 *(   GGL90visctmp(i  ,j-1,k)*mskCor(i  ,j-1,bi,bj)
750         &           + GGL90visctmp(i  ,j  ,k)*mskCor(i  ,j  ,bi,bj) )
751         &    +p8 *( ( GGL90visctmp(i-1,j-1,k)*mskCor(i-1,j-1,bi,bj)
752         &           + 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)
754         &           + GGL90visctmp(i+1,j  ,k)*mskCor(i+1,j  ,bi,bj) ) )
755         &             )/(
756         &     p4 * 2. _d 0
757         &    +p8 *( (  maskC(i-1,j-1,k,bi,bj)*mskCor(i-1,j-1,bi,bj)
758         &           +  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)
760         &           +  maskC(i+1,j  ,k,bi,bj)*mskCor(i+1,j  ,bi,bj) ) )
761         &               )*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)
763    #else
764             tmpVisc = _maskS(i,j,k,bi,bj) * halfRL
765         &          *( GGL90visctmp(i,j-1,k)
766         &           + GGL90visctmp(i,j,k) )
767    #endif
768             tmpVisc = MIN( tmpVisc , GGL90viscMax )
769             GGL90viscArV(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
770            ENDDO
771           ENDDO
772          ENDDO
773    
774    #ifdef ALLOW_DIAGNOSTICS
775          IF ( useDiagnostics ) THEN
776            CALL DIAGNOSTICS_FILL( GGL90TKE   ,'GGL90TKE',
777         &                         0,Nr, 1, bi, bj, myThid )
778            CALL DIAGNOSTICS_FILL( GGL90viscArU,'GGL90ArU',
779         &                         0,Nr, 1, bi, bj, myThid )
780            CALL DIAGNOSTICS_FILL( GGL90viscArV,'GGL90ArV',
781         &                         0,Nr, 1, bi, bj, myThid )
782            CALL DIAGNOSTICS_FILL( GGL90diffKr,'GGL90Kr ',
783         &                         0,Nr, 1, bi, bj, myThid )
784            CALL DIAGNOSTICS_FILL( TKEPrandtlNumber ,'GGL90Prl',
785         &                         0,Nr, 2, bi, bj, myThid )
786            CALL DIAGNOSTICS_FILL( GGL90mixingLength,'GGL90Lmx',
787         &                         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
817    #endif /* ALLOW_DIAGNOSTICS */
818    
819  #endif /* ALLOW_GGL90 */  #endif /* ALLOW_GGL90 */
820    
821        RETURN        RETURN
822        END        END
   
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.34
  ViewVC Help
Powered by ViewVC 1.1.22