/[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.18 by gforget, Wed Aug 11 03:32:29 2010 UTC revision 1.23 by gforget, Wed Aug 8 22:22:42 2012 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    
14  C !DESCRIPTION: \bv  C !DESCRIPTION: \bv
15  C     *==========================================================*  C     *==========================================================*
# Line 41  C !USES: =============================== Line 42  C !USES: ===============================
42    
43  C !INPUT PARAMETERS: ===================================================  C !INPUT PARAMETERS: ===================================================
44  C Routine arguments  C Routine arguments
45  C     bi, bj :: array indices on which to apply calculations  C     bi, bj :: Current tile indices
46    C     sigmaR :: Vertical gradient of iso-neutral density
47  C     myTime :: Current time in simulation  C     myTime :: Current time in simulation
48    C     myIter :: Current time-step number
49  C     myThid :: My Thread Id number  C     myThid :: My Thread Id number
50        INTEGER bi, bj        INTEGER bi, bj
51          _RL     sigmaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
52        _RL     myTime        _RL     myTime
53          INTEGER myIter
54        INTEGER myThid        INTEGER myThid
55  CEOP  CEOP
56    
# Line 53  CEOP Line 58  CEOP
58    
59  C !LOCAL VARIABLES: ====================================================  C !LOCAL VARIABLES: ====================================================
60  C Local constants  C Local constants
61  C     iMin, iMax, jMin, jMax, I, J - array computation indices  C     iMin,iMax,jMin,jMax :: index boundaries of computation domain
62  C     K, Kp1, km1, kSurf, kBottom  - vertical loop indices  C     i, j, k, kp1,km1 :: array computation indices
63  C     ab15, ab05       - weights for implicit timestepping  C     kSurf, kBottom   :: vertical indices of domain boundaries
64  C     uStarSquare      - square of friction velocity  C     explDissFac      :: explicit Dissipation Factor (in [0-1])
65  C     verticalShear    - (squared) vertical shear of horizontal velocity  C     implDissFac      :: implicit Dissipation Factor (in [0-1])
66  C     Nsquare          - squared buoyancy freqency  C     uStarSquare      :: square of friction velocity
67  C     RiNumber         - local Richardson number  C     verticalShear    :: (squared) vertical shear of horizontal velocity
68  C     KappaM           - (local) viscosity parameter (eq.10)  C     Nsquare          :: squared buoyancy freqency
69  C     KappaH           - (local) diffusivity parameter for temperature (eq.11)  C     RiNumber         :: local Richardson number
70  C     KappaE           - (local) diffusivity parameter for TKE (eq.15)  C     KappaM           :: (local) viscosity parameter (eq.10)
71  C     TKEdissipation   - dissipation of TKE  C     KappaH           :: (local) diffusivity parameter for temperature (eq.11)
72  C     GGL90mixingLength- mixing length of scheme following Banke+Delecuse  C     KappaE           :: (local) diffusivity parameter for TKE (eq.15)
73  C         rMixingLength- inverse of mixing length  C     TKEdissipation   :: dissipation of TKE
74  C     totalDepth       - thickness of water column (inverse of recip_Rcol)  C     GGL90mixingLength:: mixing length of scheme following Banke+Delecuse
75  C     TKEPrandtlNumber - here, an empirical function of the Richardson number  C         rMixingLength:: inverse of mixing length
76  C     rhoK, rhoKm1     - density at layer K and Km1 (relative to K)  C     totalDepth       :: thickness of water column (inverse of recip_Rcol)
77  C     gTKE             - right hand side of implicit equation  C     TKEPrandtlNumber :: here, an empirical function of the Richardson number
78        INTEGER iMin ,iMax ,jMin ,jMax        INTEGER iMin ,iMax ,jMin ,jMax
79        INTEGER I, J, K, Kp1, Km1, kSurf, kBottom        INTEGER i, j, k, kp1, km1, kSurf, kBottom
80        _RL     ab15, ab05        _RL     explDissFac, implDissFac
81        _RL     uStarSquare        _RL     uStarSquare
82        _RL     verticalShear        _RL     verticalShear
83        _RL     KappaM, KappaH        _RL     KappaM, KappaH
# Line 90  c     _RL     SQRTTKE Line 95  c     _RL     SQRTTKE
95        _RL     rMixingLength    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     rMixingLength    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
96        _RL     mxLength_Dn      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     mxLength_Dn      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
97        _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)  
98        _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL     gTKE             (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
99        _RL     GGL90visctmp     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     GGL90visctmp     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
100  C-    tri-diagonal matrix  C-    tri-diagonal matrix
101        _RL     a(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     a3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
102        _RL     b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     b3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
103        _RL     c(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     c3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
104        INTEGER errCode        INTEGER errCode
105  #ifdef ALLOW_GGL90_HORIZDIFF  #ifdef ALLOW_GGL90_HORIZDIFF
106  C-    xA, yA   - area of lateral faces  C     xA, yA   :: area of lateral faces
107  C-    dfx, dfy - diffusive flux across lateral faces  C     dfx, dfy :: diffusive flux across lateral faces
108    C     gTKE     :: right hand side of diffusion equation
109        _RL     xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
110        _RL     yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
111        _RL     dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
112        _RL     dfy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     dfy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
113          _RL    gTKE(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
114  #endif /* ALLOW_GGL90_HORIZDIFF */  #endif /* ALLOW_GGL90_HORIZDIFF */
115  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
116        _RL p4, p8, p16        _RL p4, p8, p16
# Line 123  C     set separate time step (should be Line 127  C     set separate time step (should be
127        deltaTggl90 = dTtracerLev(1)        deltaTggl90 = dTtracerLev(1)
128    
129        kSurf = 1        kSurf = 1
130  C     implicit timestepping weights for dissipation  C     explicit/implicit timestepping weights for dissipation
131        ab15 =  1.5 _d 0        explDissFac = 0. _d 0
132        ab05 = -0.5 _d 0        implDissFac = 1. _d 0 - explDissFac
       ab15 =  1. _d 0  
       ab05 =  0. _d 0  
133    
134  C     Initialize local fields  C     Initialize local fields
135        DO K = 1, Nr        DO k = 1, Nr
136         DO J=1-Oly,sNy+Oly         DO j=1-OLy,sNy+OLy
137          DO I=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
138           gTKE(I,J,K)              = 0. _d 0           KappaE(i,j,k)            = 0. _d 0
139           KappaE(I,J,K)            = 0. _d 0           TKEPrandtlNumber(i,j,k)  = 1. _d 0
140           TKEPrandtlNumber(I,J,K)  = 1. _d 0           GGL90mixingLength(i,j,k) = GGL90mixingLengthMin
141           GGL90mixingLength(I,J,K) = GGL90mixingLengthMin           GGL90visctmp(i,j,k)      = 0. _d 0
142           GGL90visctmp(I,J,K)      = 0. _d 0  #ifndef SOLVE_DIAGONAL_LOWMEMORY
143             a3d(i,j,k) = 0. _d 0
144             b3d(i,j,k) = 1. _d 0
145             c3d(i,j,k) = 0. _d 0
146    #endif
147          ENDDO          ENDDO
148         ENDDO         ENDDO
149        ENDDO        ENDDO
150        DO J=1-Oly,sNy+Oly        DO j=1-OLy,sNy+OLy
151         DO I=1-Olx,sNx+Olx         DO i=1-OLx,sNx+OLx
152          rhoK(I,J)          = 0. _d 0          totalDepth(i,j)    = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
         rhoKm1(I,J)        = 0. _d 0  
         totalDepth(I,J)    = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)  
153          rMixingLength(i,j,1) = 0. _d 0          rMixingLength(i,j,1) = 0. _d 0
154          mxLength_Dn(I,J,1) = GGL90mixingLengthMin          mxLength_Dn(i,j,1) = GGL90mixingLengthMin
155          SQRTTKE(i,j,1) = SQRT( GGL90TKE(i,j,1,bi,bj) )          SQRTTKE(i,j,1) = SQRT( GGL90TKE(i,j,1,bi,bj) )
156         ENDDO         ENDDO
157        ENDDO        ENDDO
158    
159  C     start k-loop  C     start k-loop
160        DO K = 2, Nr        DO k = 2, Nr
161         Km1 = K-1  c      km1 = k-1
162  c      Kp1 = MIN(Nr,K+1)  c      kp1 = MIN(Nr,k+1)
163         CALL FIND_RHO_2D(         DO j=jMin,jMax
164       I      iMin, iMax, jMin, jMax, K,          DO i=iMin,iMax
165       I      theta(1-OLx,1-OLy,Km1,bi,bj), salt(1-OLx,1-OLy,Km1,bi,bj),           SQRTTKE(i,j,k)=SQRT( GGL90TKE(i,j,k,bi,bj) )
      O      rhoKm1,  
      I      Km1, bi, bj, myThid )  
   
        CALL FIND_RHO_2D(  
      I      iMin, iMax, jMin, jMax, K,  
      I      theta(1-OLx,1-OLy,K,bi,bj), salt(1-OLx,1-OLy,K,bi,bj),  
      O      rhoK,  
      I      K, bi, bj, myThid )  
        DO J=jMin,jMax  
         DO I=iMin,iMax  
          SQRTTKE(i,j,k)=SQRT( GGL90TKE(I,J,K,bi,bj) )  
166    
167  C     buoyancy frequency  C     buoyancy frequency
168           Nsquare(i,j,k) = - gravity*recip_rhoConst*recip_drC(K)           Nsquare(i,j,k) = gravity*gravitySign*recip_rhoConst
169       &        * ( rhoKm1(I,J) - rhoK(I,J) )*maskC(I,J,K,bi,bj)       &                  * sigmaR(i,j,k)
170  cC     vertical shear term (dU/dz)^2+(dV/dz)^2  cC     vertical shear term (dU/dz)^2+(dV/dz)^2
171  c         tempU= .5 _d 0*( uVel(I,J,Km1,bi,bj)+uVel(I+1,J,Km1,bi,bj)  c         tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)
172  c     &                 -( uVel(I,J,K  ,bi,bj)+uVel(I+1,J,K  ,bi,bj)) )  c     &                 -( uVel(i,j,k  ,bi,bj)+uVel(i+1,j,k  ,bi,bj)) )
173  c     &        *recip_drC(K)  c     &        *recip_drC(k)
174  c         tempV= .5 _d 0*( vVel(I,J,Km1,bi,bj)+vVel(I,J+1,Km1,bi,bj)  c         tempV= .5 _d 0*( vVel(i,j,km1,bi,bj)+vVel(i,j+1,km1,bi,bj)
175  c     &                 -( vVel(I,J,K  ,bi,bj)+vVel(I,J+1,K  ,bi,bj)) )  c     &                 -( vVel(i,j,k  ,bi,bj)+vVel(i,j+1,k  ,bi,bj)) )
176  c     &        *recip_drC(K)  c     &        *recip_drC(k)
177  c         verticalShear = tempU*tempU + tempV*tempV  c         verticalShear = tempU*tempU + tempV*tempV
178  c         RiNumber   = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)  c         RiNumber   = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)
179  cC     compute Prandtl number (always greater than 0)  cC     compute Prandtl number (always greater than 0)
180  c         prTemp = 1. _d 0  c         prTemp = 1. _d 0
181  c         IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber  c         IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
182  c         TKEPrandtlNumber(I,J,K) = MIN(10. _d 0,prTemp)  c         TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
183  C     mixing length  C     mixing length
184           GGL90mixingLength(I,J,K) = SQRTTWO *           GGL90mixingLength(i,j,k) = SQRTTWO *
185       &        SQRTTKE(i,j,k)/SQRT( MAX(Nsquare(i,j,k),GGL90eps) )       &        SQRTTKE(i,j,k)/SQRT( MAX(Nsquare(i,j,k),GGL90eps) )
186          ENDDO          ENDDO
187         ENDDO         ENDDO
188        ENDDO        ENDDO
189    
190    C- ensure mixing between first and second level
191          IF (mxlSurfFlag) THEN
192           DO j=jMin,jMax
193            DO i=iMin,iMax
194             GGL90mixingLength(i,j,2)=drF(1)
195            ENDDO
196           ENDDO
197          ENDIF
198    
199  C- Impose upper and lower bound for mixing length  C- Impose upper and lower bound for mixing length
200        IF ( mxlMaxFlag .EQ. 0 ) THEN        IF ( mxlMaxFlag .EQ. 0 ) THEN
201  C-  
202         DO k=2,Nr         DO k=2,Nr
203          DO J=jMin,jMax          DO j=jMin,jMax
204           DO I=iMin,iMax           DO i=iMin,iMax
205            MaxLength=totalDepth(I,J)            MaxLength=totalDepth(i,j)
206            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
207       &                                   MaxLength)       &                                   MaxLength)
208           ENDDO           ENDDO
209          ENDDO          ENDDO
210         ENDDO         ENDDO
211    
212         DO k=2,Nr         DO k=2,Nr
213          DO J=jMin,jMax          DO j=jMin,jMax
214           DO I=iMin,iMax           DO i=iMin,iMax
215            GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
216       &                                   GGL90mixingLengthMin)       &                                   GGL90mixingLengthMin)
217            rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K)            rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
218           ENDDO           ENDDO
219          ENDDO          ENDDO
220         ENDDO         ENDDO
221    
222        ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN
223  C-  
224         DO k=2,Nr         DO k=2,Nr
225          DO J=jMin,jMax          DO j=jMin,jMax
226           DO I=iMin,iMax           DO i=iMin,iMax
227            MaxLength=MIN(Ro_surf(I,J,bi,bj)-rF(k),rF(k)-R_low(I,J,bi,bj))            MaxLength=MIN(Ro_surf(i,j,bi,bj)-rF(k),rF(k)-R_low(i,j,bi,bj))
228  c         MaxLength=MAX(MaxLength,20. _d 0)  c         MaxLength=MAX(MaxLength,20. _d 0)
229            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
230       &                                   MaxLength)       &                                   MaxLength)
231           ENDDO           ENDDO
232          ENDDO          ENDDO
233         ENDDO         ENDDO
234    
235         DO k=2,Nr         DO k=2,Nr
236          DO J=jMin,jMax          DO j=jMin,jMax
237           DO I=iMin,iMax           DO i=iMin,iMax
238            GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
239       &                                   GGL90mixingLengthMin)       &                                   GGL90mixingLengthMin)
240            rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K)            rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
241           ENDDO           ENDDO
242          ENDDO          ENDDO
243         ENDDO         ENDDO
244    
245        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN
246  C-  
 cgf ensure mixing between first and second level  
 c      DO J=jMin,jMax  
 c        DO I=iMin,iMax  
 c         GGL90mixingLength(I,J,2)=drF(1)  
 c        ENDDO  
 c      ENDDO  
 cgf  
247         DO k=2,Nr         DO k=2,Nr
248          DO J=jMin,jMax          DO j=jMin,jMax
249           DO I=iMin,iMax           DO i=iMin,iMax
250            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
251       &        GGL90mixingLength(I,J,K-1)+drF(k-1))       &        GGL90mixingLength(i,j,k-1)+drF(k-1))
252           ENDDO           ENDDO
253          ENDDO          ENDDO
254         ENDDO         ENDDO
255         DO J=jMin,jMax         DO j=jMin,jMax
256          DO I=iMin,iMax          DO i=iMin,iMax
257            GGL90mixingLength(I,J,Nr) = MIN(GGL90mixingLength(I,J,Nr),            GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
258       &       GGL90mixingLengthMin+drF(Nr))       &       GGL90mixingLengthMin+drF(Nr))
259          ENDDO          ENDDO
260         ENDDO         ENDDO
261         DO k=Nr-1,2,-1         DO k=Nr-1,2,-1
262          DO J=jMin,jMax          DO j=jMin,jMax
263           DO I=iMin,iMax           DO i=iMin,iMax
264            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
265       &        GGL90mixingLength(I,J,K+1)+drF(k))       &        GGL90mixingLength(i,j,k+1)+drF(k))
266           ENDDO           ENDDO
267          ENDDO          ENDDO
268         ENDDO         ENDDO
269    
270         DO k=2,Nr         DO k=2,Nr
271          DO J=jMin,jMax          DO j=jMin,jMax
272           DO I=iMin,iMax           DO i=iMin,iMax
273            GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
274       &                                   GGL90mixingLengthMin)       &                                   GGL90mixingLengthMin)
275            rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K)            rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
276           ENDDO           ENDDO
277          ENDDO          ENDDO
278         ENDDO         ENDDO
279    
280        ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN
281  C-  
282         DO k=2,Nr         DO k=2,Nr
283          DO J=jMin,jMax          DO j=jMin,jMax
284           DO I=iMin,iMax           DO i=iMin,iMax
285            mxLength_Dn(I,J,K) = MIN(GGL90mixingLength(I,J,K),            mxLength_Dn(i,j,k) = MIN(GGL90mixingLength(i,j,k),
286       &        mxLength_Dn(I,J,K-1)+drF(k-1))       &        mxLength_Dn(i,j,k-1)+drF(k-1))
287           ENDDO           ENDDO
288          ENDDO          ENDDO
289         ENDDO         ENDDO
290         DO J=jMin,jMax         DO j=jMin,jMax
291          DO I=iMin,iMax          DO i=iMin,iMax
292            GGL90mixingLength(I,J,Nr) = MIN(GGL90mixingLength(I,J,Nr),            GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
293       &       GGL90mixingLengthMin+drF(Nr))       &       GGL90mixingLengthMin+drF(Nr))
294          ENDDO          ENDDO
295         ENDDO         ENDDO
296         DO k=Nr-1,2,-1         DO k=Nr-1,2,-1
297          DO J=jMin,jMax          DO j=jMin,jMax
298           DO I=iMin,iMax           DO i=iMin,iMax
299            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
300       &        GGL90mixingLength(I,J,K+1)+drF(k))       &        GGL90mixingLength(i,j,k+1)+drF(k))
301           ENDDO           ENDDO
302          ENDDO          ENDDO
303         ENDDO         ENDDO
304    
305         DO k=2,Nr         DO k=2,Nr
306          DO J=jMin,jMax          DO j=jMin,jMax
307           DO I=iMin,iMax           DO i=iMin,iMax
308            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
309       &                                  mxLength_Dn(I,J,K))       &                                  mxLength_Dn(i,j,k))
310            tmpmlx = SQRT( GGL90mixingLength(I,J,K)*mxLength_Dn(I,J,K) )            tmpmlx = SQRT( GGL90mixingLength(i,j,k)*mxLength_Dn(i,j,k) )
311            tmpmlx = MAX( tmpmlx, GGL90mixingLengthMin)            tmpmlx = MAX( tmpmlx, GGL90mixingLengthMin)
312            rMixingLength(I,J,K) = 1. _d 0 / tmpmlx            rMixingLength(i,j,k) = 1. _d 0 / tmpmlx
313           ENDDO           ENDDO
314          ENDDO          ENDDO
315         ENDDO         ENDDO
# Line 325  C- Line 320  C-
320    
321  C- Impose minimum mixing length (to avoid division by zero)  C- Impose minimum mixing length (to avoid division by zero)
322  c      DO k=2,Nr  c      DO k=2,Nr
323  c      DO J=jMin,jMax  c      DO j=jMin,jMax
324  c       DO I=iMin,iMax  c       DO i=iMin,iMax
325  c        GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),  c        GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
326  c    &        GGL90mixingLengthMin)  c    &        GGL90mixingLengthMin)
327  c        rMixingLength(I,J,K) = 1. _d 0 /GGL90mixingLength(I,J,K)  c        rMixingLength(i,j,k) = 1. _d 0 /GGL90mixingLength(i,j,k)
328  c       ENDDO  c       ENDDO
329  c      ENDDO  c      ENDDO
330  c     ENDDO  c     ENDDO
331    
   
332        DO k=2,Nr        DO k=2,Nr
333         Km1 = K-1         km1 = k-1
        DO J=jMin,jMax  
         DO I=iMin,iMax  
 C     vertical shear term (dU/dz)^2+(dV/dz)^2  
          tempU= .5 _d 0*( 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 _d 0*( 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(i,j,k),0. _d 0)/(verticalShear+GGL90eps)  
 C     compute Prandtl number (always greater than 0)  
          prTemp = 1. _d 0  
          IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber  
          TKEPrandtlNumber(I,J,K) = MIN(10. _d 0,prTemp)  
 c         TKEPrandtlNumber(I,J,K) = 1. _d 0  
   
 C     viscosity and diffusivity  
          KappaM = GGL90ck*GGL90mixingLength(I,J,K)*SQRTTKE(i,j,k)  
          GGL90visctmp(I,J,K) = MAX(KappaM,diffKrNrT(k))  
      &                            * maskC(I,J,K,bi,bj)  
 c        note: storing GGL90visctmp like this, and using it later to compute  
 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)  
   
 C     dissipation term  
          TKEdissipation = ab05*GGL90ceps  
      &        *SQRTTKE(i,j,k)*rMixingLength(I,J,K)  
      &        *GGL90TKE(I,J,K,bi,bj)  
 C     sum up contributions to form the right hand side  
          gTKE(I,J,K) = GGL90TKE(I,J,K,bi,bj)  
      &        + deltaTggl90*(  
      &        + KappaM*verticalShear  
      &        - KappaH*Nsquare(i,j,k)  
      &        - TKEdissipation  
      &        )  
         ENDDO  
        ENDDO  
       ENDDO  
334    
335    #ifdef ALLOW_GGL90_HORIZDIFF
336           IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
337  C     horizontal diffusion of TKE (requires an exchange in  C     horizontal diffusion of TKE (requires an exchange in
338  C      do_fields_blocking_exchanges)  C      do_fields_blocking_exchanges)
 #ifdef ALLOW_GGL90_HORIZDIFF  
       IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN  
        DO K=2,Nr  
339  C     common factors  C     common factors
340          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
341           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
342            xA(i,j) = _dyG(i,j,bi,bj)            xA(i,j) = _dyG(i,j,bi,bj)
343       &         *drF(k)*_hFacW(i,j,k,bi,bj)       &         *drF(k)*_hFacW(i,j,k,bi,bj)
344            yA(i,j) = _dxG(i,j,bi,bj)            yA(i,j) = _dxG(i,j,bi,bj)
# Line 395  C     common factors Line 347  C     common factors
347          ENDDO          ENDDO
348  C     Compute diffusive fluxes  C     Compute diffusive fluxes
349  C     ... across x-faces  C     ... across x-faces
350          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
351           dfx(1-Olx,j)=0. _d 0           dfx(1-OLx,j)=0. _d 0
352           DO i=1-Olx+1,sNx+Olx           DO i=1-OLx+1,sNx+OLx
353            dfx(i,j) = -GGL90diffTKEh*xA(i,j)            dfx(i,j) = -GGL90diffTKEh*xA(i,j)
354       &      *_recip_dxC(i,j,bi,bj)       &      *_recip_dxC(i,j,bi,bj)
355       &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i-1,j,k,bi,bj))       &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i-1,j,k,bi,bj))
# Line 405  C     ... across x-faces Line 357  C     ... across x-faces
357           ENDDO           ENDDO
358          ENDDO          ENDDO
359  C     ... across y-faces  C     ... across y-faces
360          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
361           dfy(i,1-Oly)=0. _d 0           dfy(i,1-OLy)=0. _d 0
362          ENDDO          ENDDO
363          DO j=1-Oly+1,sNy+Oly          DO j=1-OLy+1,sNy+OLy
364           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
365            dfy(i,j) = -GGL90diffTKEh*yA(i,j)            dfy(i,j) = -GGL90diffTKEh*yA(i,j)
366       &      *_recip_dyC(i,j,bi,bj)       &      *_recip_dyC(i,j,bi,bj)
367       &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i,j-1,k,bi,bj))       &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i,j-1,k,bi,bj))
# Line 419  C     ... across y-faces Line 371  C     ... across y-faces
371           ENDDO           ENDDO
372          ENDDO          ENDDO
373  C     Compute divergence of fluxes  C     Compute divergence of fluxes
374          DO j=1-Oly,sNy+Oly-1          DO j=1-OLy,sNy+OLy-1
375           DO i=1-Olx,sNx+Olx-1           DO i=1-OLx,sNx+OLx-1
376            gTKE(i,j,k)=gTKE(i,j,k)            gTKE(i,j) =
377       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj)       &    -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj)
378       &         *( (dfx(i+1,j)-dfx(i,j))       &         *( (dfx(i+1,j)-dfx(i,j))
379       &           +(dfy(i,j+1)-dfy(i,j))       &           +(dfy(i,j+1)-dfy(i,j))
380       &           )*deltaTggl90       &          )
381           ENDDO           ENDDO
382          ENDDO          ENDDO
383  C       end of k-loop  C      end if GGL90diffTKEh .eq. 0.
384           ENDIF
385    #endif /* ALLOW_GGL90_HORIZDIFF */
386    
387           DO j=jMin,jMax
388            DO i=iMin,iMax
389    C     vertical shear term (dU/dz)^2+(dV/dz)^2
390             tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)
391         &                 -( uVel(i,j,k  ,bi,bj)+uVel(i+1,j,k  ,bi,bj)) )
392         &        *recip_drC(k)
393             tempV= .5 _d 0*( vVel(i,j,km1,bi,bj)+vVel(i,j+1,km1,bi,bj)
394         &                 -( vVel(i,j,k  ,bi,bj)+vVel(i,j+1,k  ,bi,bj)) )
395         &        *recip_drC(k)
396             verticalShear = tempU*tempU + tempV*tempV
397             RiNumber = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)
398    C     compute Prandtl number (always greater than 0)
399             prTemp = 1. _d 0
400             IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
401             TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
402    c         TKEPrandtlNumber(i,j,k) = 1. _d 0
403    
404    C     viscosity and diffusivity
405             KappaM = GGL90ck*GGL90mixingLength(i,j,k)*SQRTTKE(i,j,k)
406             GGL90visctmp(i,j,k) = MAX(KappaM,diffKrNrT(k))
407         &                            * maskC(i,j,k,bi,bj)
408    c        note: storing GGL90visctmp like this, and using it later to compute
409    c              GGL9rdiffKr etc. is robust in case of smoothing (e.g. see OPA)
410             KappaM = MAX(KappaM,viscArNr(k)) * maskC(i,j,k,bi,bj)
411             KappaH = KappaM/TKEPrandtlNumber(i,j,k)
412             KappaE(i,j,k) = GGL90alpha * KappaM * maskC(i,j,k,bi,bj)
413    
414    C     dissipation term
415             TKEdissipation = explDissFac*GGL90ceps
416         &        *SQRTTKE(i,j,k)*rMixingLength(i,j,k)
417         &        *GGL90TKE(i,j,k,bi,bj)
418    C     partial update with sum of explicit contributions
419             GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
420         &        + deltaTggl90*(
421         &        + KappaM*verticalShear
422         &        - KappaH*Nsquare(i,j,k)
423         &        - TKEdissipation
424         &        )
425            ENDDO
426         ENDDO         ENDDO
427  C     end if GGL90diffTKEh .eq. 0.  
428        ENDIF  #ifdef ALLOW_GGL90_HORIZDIFF
429           IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
430    C--    Add horiz. diffusion tendency
431            DO j=jMin,jMax
432             DO i=iMin,iMax
433              GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
434         &                          + gTKE(i,j)*deltaTggl90
435             ENDDO
436            ENDDO
437           ENDIF
438  #endif /* ALLOW_GGL90_HORIZDIFF */  #endif /* ALLOW_GGL90_HORIZDIFF */
439    
440    C--   end of k loop
441          ENDDO
442    
443  C     ============================================  C     ============================================
444  C     Implicit time step to update TKE for k=1,Nr;  C     Implicit time step to update TKE for k=1,Nr;
445  C     TKE(Nr+1)=0 by default  C     TKE(Nr+1)=0 by default
# Line 442  C     set up matrix Line 448  C     set up matrix
448  C--   Lower diagonal  C--   Lower diagonal
449        DO j=jMin,jMax        DO j=jMin,jMax
450         DO i=iMin,iMax         DO i=iMin,iMax
451           a(i,j,1) = 0. _d 0           a3d(i,j,1) = 0. _d 0
452         ENDDO         ENDDO
453        ENDDO        ENDDO
454        DO k=2,Nr        DO k=2,Nr
# Line 452  C--   Lower diagonal Line 458  C--   Lower diagonal
458  C-    We keep recip_hFacC in the diffusive flux calculation,  C-    We keep recip_hFacC in the diffusive flux calculation,
459  C-    but no hFacC in TKE volume control  C-    but no hFacC in TKE volume control
460  C-    No need for maskC(k-1) with recip_hFacC(k-1)  C-    No need for maskC(k-1) with recip_hFacC(k-1)
461           a(i,j,k) = -deltaTggl90           a3d(i,j,k) = -deltaTggl90
462       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)
463       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))
464       &        *recip_drC(k)*maskC(i,j,k,bi,bj)       &        *recip_drC(k)*maskC(i,j,k,bi,bj)
# Line 462  C-    No need for maskC(k-1) with recip_ Line 468  C-    No need for maskC(k-1) with recip_
468  C--   Upper diagonal  C--   Upper diagonal
469        DO j=jMin,jMax        DO j=jMin,jMax
470         DO i=iMin,iMax         DO i=iMin,iMax
471           c(i,j,1)  = 0. _d 0           c3d(i,j,1)  = 0. _d 0
472         ENDDO         ENDDO
473        ENDDO        ENDDO
474        DO k=2,Nr        DO k=2,Nr
# Line 472  C--   Upper diagonal Line 478  C--   Upper diagonal
478  C-    We keep recip_hFacC in the diffusive flux calculation,  C-    We keep recip_hFacC in the diffusive flux calculation,
479  C-    but no hFacC in TKE volume control  C-    but no hFacC in TKE volume control
480  C-    No need for maskC(k) with recip_hFacC(k)  C-    No need for maskC(k) with recip_hFacC(k)
481            c(i,j,k) = -deltaTggl90            c3d(i,j,k) = -deltaTggl90
482       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)
483       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))
484       &        *recip_drC(k)*maskC(i,j,k-1,bi,bj)       &        *recip_drC(k)*maskC(i,j,k-1,bi,bj)
# Line 484  C--   Center diagonal Line 490  C--   Center diagonal
490         km1 = MAX(k-1,1)         km1 = MAX(k-1,1)
491         DO j=jMin,jMax         DO j=jMin,jMax
492          DO i=iMin,iMax          DO i=iMin,iMax
493            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)
494       &        + ab15*deltaTggl90*GGL90ceps*SQRTTKE(I,J,K)       &        + implDissFac*deltaTggl90*GGL90ceps*SQRTTKE(i,j,k)
495       &        * rMixingLength(I,J,K)       &        * rMixingLength(i,j,k)
496       &        * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)       &        * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)
497           ENDDO           ENDDO
498         ENDDO         ENDDO
# Line 495  C     end set up matrix Line 501  C     end set up matrix
501    
502  C     Apply boundary condition  C     Apply boundary condition
503        kp1 = MIN(Nr,kSurf+1)        kp1 = MIN(Nr,kSurf+1)
504        DO J=jMin,jMax        DO j=jMin,jMax
505         DO I=iMin,iMax         DO i=iMin,iMax
506  C     estimate friction velocity uStar from surface forcing  C     estimate friction velocity uStar from surface forcing
507          uStarSquare = SQRT(          uStarSquare = SQRT(
508       &    ( .5 _d 0*( surfaceForcingU(I,  J,  bi,bj)       &    ( .5 _d 0*( surfaceForcingU(i,  j,  bi,bj)
509       &              + surfaceForcingU(I+1,J,  bi,bj) ) )**2       &              + surfaceForcingU(i+1,j,  bi,bj) ) )**2
510       &  + ( .5 _d 0*( surfaceForcingV(I,  J,  bi,bj)       &  + ( .5 _d 0*( surfaceForcingV(i,  j,  bi,bj)
511       &              + surfaceForcingV(I,  J+1,bi,bj) ) )**2       &              + surfaceForcingV(i,  j+1,bi,bj) ) )**2
512       &                     )       &                     )
513  C     Dirichlet surface boundary condition for TKE  C     Dirichlet surface boundary condition for TKE
514          gTKE(I,J,kSurf) = MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)          GGL90TKE(i,j,kSurf,bi,bj) = maskC(i,j,kSurf,bi,bj)
515       &                     *maskC(I,J,kSurf,bi,bj)       &           *MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)
516          gTKE(i,j,kp1) = gTKE(i,j,kp1)          GGL90TKE(i,j,kp1,bi,bj) = GGL90TKE(i,j,kp1,bi,bj)
517       &                - a(i,j,kp1)*gTKE(i,j,kSurf)       &               - a3d(i,j,kp1)*GGL90TKE(i,j,kSurf,bi,bj)
518          a(i,j,kp1) = 0. _d 0          a3d(i,j,kp1) = 0. _d 0
519  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom
520          kBottom   = MAX(kLowC(I,J,bi,bj),1)          kBottom   = MAX(kLowC(i,j,bi,bj),1)
521          gTKE(I,J,kBottom) = gTKE(I,J,kBottom)          GGL90TKE(i,j,kBottom,bi,bj) = GGL90TKE(i,j,kBottom,bi,bj)
522       &                    - GGL90TKEbottom*c(I,J,kBottom)       &                              - GGL90TKEbottom*c3d(i,j,kBottom)
523          c(I,J,kBottom) = 0. _d 0          c3d(i,j,kBottom) = 0. _d 0
524         ENDDO         ENDDO
525        ENDDO        ENDDO
526    
527  C     solve tri-diagonal system, and store solution on gTKE (previously rhs)  C     solve tri-diagonal system
528        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
529       I                        a, b, c,       I                        a3d, b3d, c3d,
530       U                        gTKE,       U                        GGL90TKE,
531       O                        errCode,       O                        errCode,
532       I                        1, 1, myThid )       I                        bi, bj, myThid )
533    
534  C     now update TKE        DO k=1,Nr
535        DO K=1,Nr         DO j=jMin,jMax
536         DO J=jMin,jMax          DO i=iMin,iMax
         DO I=iMin,iMax  
537  C     impose minimum TKE to avoid numerical undershoots below zero  C     impose minimum TKE to avoid numerical undershoots below zero
538           GGL90TKE(I,J,K,bi,bj) = MAX( gTKE(I,J,K), GGL90TKEmin )           GGL90TKE(i,j,k,bi,bj) = maskC(i,j,k,bi,bj)
539       &        * maskC(I,J,K,bi,bj)       &                  *MAX( GGL90TKE(i,j,k,bi,bj), GGL90TKEmin )
540          ENDDO          ENDDO
541         ENDDO         ENDDO
542        ENDDO        ENDDO
# Line 539  C     impose minimum TKE to avoid numeri Line 544  C     impose minimum TKE to avoid numeri
544  C     end of time step  C     end of time step
545  C     ===============================  C     ===============================
546    
547        DO K=2,Nr        DO k=2,Nr
548         DO J=1,sNy         DO j=1,sNy
549          DO I=1,sNx          DO i=1,sNx
550  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
551           tmpVisc=           tmpVisc=
552       &  (       &  (
# Line 566  C     =============================== Line 571  C     ===============================
571       &       +       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))
572       &  )*maskC(i,j,k,bi,bj)*mskCor(i,j,bi,bj)       &  )*maskC(i,j,k,bi,bj)*mskCor(i,j,bi,bj)
573  #else  #else
574           tmpVisc = GGL90visctmp(I,J,K)           tmpVisc = GGL90visctmp(i,j,k)
575  #endif  #endif
576           tmpVisc = MIN(tmpVisc/TKEPrandtlNumber(i,j,k),GGL90diffMax)           tmpVisc = MIN(tmpVisc/TKEPrandtlNumber(i,j,k),GGL90diffMax)
577           GGL90diffKr(I,J,K,bi,bj)= MAX( tmpVisc , diffKrNrT(k) )           GGL90diffKr(i,j,k,bi,bj)= MAX( tmpVisc , diffKrNrT(k) )
578          ENDDO          ENDDO
579         ENDDO         ENDDO
580        ENDDO        ENDDO
581    
582          DO k=2,Nr
583           DO j=1,sNy
584        DO K=2,Nr          DO i=1,sNx+1
        DO J=1,sNy  
         DO I=1,sNx+1  
585  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
586          tmpVisc =          tmpVisc =
587       & (       & (
588       &   p4 *(GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)       &   p4 *(GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
589       &       +GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj))       &       +GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj))
# Line 603  C     =============================== Line 606  C     ===============================
606       &                            +GGL90visctmp(i-1,j,k))       &                            +GGL90visctmp(i-1,j,k))
607       &                   )       &                   )
608  #endif  #endif
609          tmpVisc = MIN( tmpVisc , GGL90viscMax )          tmpVisc = MIN( tmpVisc , GGL90viscMax )
610          GGL90viscArU(i,j,k,bi,bj) = MAX( tmpVisc , viscArNr(k)  )          GGL90viscArU(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
611          ENDDO          ENDDO
612         ENDDO         ENDDO
613        ENDDO        ENDDO
614    
615          DO k=2,Nr
616        DO K=2,Nr         DO j=1,sNy+1
617         DO J=1,sNy+1          DO i=1,sNx
         DO I=1,sNx  
618  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
619          tmpVisc =          tmpVisc =
620       & (       & (
# Line 639  C     =============================== Line 641  C     ===============================
641    
642  #endif  #endif
643          tmpVisc = MIN( tmpVisc , GGL90viscMax )          tmpVisc = MIN( tmpVisc , GGL90viscMax )
644          GGL90viscArV(i,j,k,bi,bj) = MAX( tmpVisc , viscArNr(k)  )          GGL90viscArV(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
645          ENDDO          ENDDO
646         ENDDO         ENDDO
647        ENDDO        ENDDO
Legend:
Removed from v.1.18  
changed lines
  Added in v.1.23
  ViewVC Help
Powered by ViewVC 1.1.22