/[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.17 by gforget, Mon Aug 9 20:34:02 2010 UTC revision 1.21 by jmc, Wed Jun 27 22:39:09 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    C     rhoK, rhoKm1     :: density at layer k and km1 (relative to k)
79        INTEGER iMin ,iMax ,jMin ,jMax        INTEGER iMin ,iMax ,jMin ,jMax
80        INTEGER I, J, K, Kp1, Km1, kSurf, kBottom        INTEGER i, j, k, kp1, km1, kSurf, kBottom
81        _RL     ab15, ab05        _RL     explDissFac, implDissFac
82        _RL     uStarSquare        _RL     uStarSquare
83        _RL     verticalShear        _RL     verticalShear
84        _RL     KappaM, KappaH        _RL     KappaM, KappaH
# Line 93  c     _RL     SQRTTKE Line 99  c     _RL     SQRTTKE
99        _RL     rhoK             (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     rhoK             (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
100        _RL     rhoKm1           (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)
       _RL     gTKE             (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
102        _RL     GGL90visctmp     (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  C-    tri-diagonal matrix
104        _RL     a(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     a3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
105        _RL     b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     b3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
106        _RL     c(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     c3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
107        INTEGER errCode        INTEGER errCode
108  #ifdef ALLOW_GGL90_HORIZDIFF  #ifdef ALLOW_GGL90_HORIZDIFF
109  C-    xA, yA   - area of lateral faces  C     xA, yA   :: area of lateral faces
110  C-    dfx, dfy - diffusive flux across lateral faces  C     dfx, dfy :: diffusive flux across lateral faces
111    C     gTKE     :: right hand side of diffusion equation
112        _RL     xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
113        _RL     yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
114        _RL     dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
115        _RL     dfy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     dfy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
116          _RL    gTKE(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
117  #endif /* ALLOW_GGL90_HORIZDIFF */  #endif /* ALLOW_GGL90_HORIZDIFF */
118  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
119        _RL p4, p8, p16        _RL p4, p8, p16
# Line 123  C     set separate time step (should be Line 130  C     set separate time step (should be
130        deltaTggl90 = dTtracerLev(1)        deltaTggl90 = dTtracerLev(1)
131    
132        kSurf = 1        kSurf = 1
133  C     implicit timestepping weights for dissipation  C     explicit/implicit timestepping weights for dissipation
134        ab15 =  1.5 _d 0        explDissFac = 0. _d 0
135        ab05 = -0.5 _d 0        implDissFac = 1. _d 0 - explDissFac
       ab15 =  1. _d 0  
       ab05 =  0. _d 0  
136    
137  C     Initialize local fields  C     Initialize local fields
138        DO K = 1, Nr        DO k = 1, Nr
139         DO J=1-Oly,sNy+Oly         DO j=1-OLy,sNy+OLy
140          DO I=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
141           gTKE(I,J,K)              = 0. _d 0           KappaE(i,j,k)            = 0. _d 0
142           KappaE(I,J,K)            = 0. _d 0           TKEPrandtlNumber(i,j,k)  = 1. _d 0
143           TKEPrandtlNumber(I,J,K)  = 1. _d 0           GGL90mixingLength(i,j,k) = GGL90mixingLengthMin
144           GGL90mixingLength(I,J,K) = GGL90mixingLengthMin           GGL90visctmp(i,j,k)      = 0. _d 0
145           GGL90visctmp(I,J,K)      = 0. _d 0  #ifndef SOLVE_DIAGONAL_LOWMEMORY
146             a3d(i,j,k) = 0. _d 0
147             b3d(i,j,k) = 1. _d 0
148             c3d(i,j,k) = 0. _d 0
149    #endif
150          ENDDO          ENDDO
151         ENDDO         ENDDO
152        ENDDO        ENDDO
153        DO J=1-Oly,sNy+Oly        DO j=1-OLy,sNy+OLy
154         DO I=1-Olx,sNx+Olx         DO i=1-OLx,sNx+OLx
155          rhoK(I,J)          = 0. _d 0          rhoK(i,j)          = 0. _d 0
156          rhoKm1(I,J)        = 0. _d 0          rhoKm1(i,j)        = 0. _d 0
157          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)
158          rMixingLength(i,j,1) = 0. _d 0          rMixingLength(i,j,1) = 0. _d 0
159          mxLength_Dn(I,J,1) = GGL90mixingLengthMin          mxLength_Dn(i,j,1) = GGL90mixingLengthMin
160          SQRTTKE(i,j,1) = SQRT( GGL90TKE(i,j,1,bi,bj) )          SQRTTKE(i,j,1) = SQRT( GGL90TKE(i,j,1,bi,bj) )
161         ENDDO         ENDDO
162        ENDDO        ENDDO
163    
164  C     start k-loop  C     start k-loop
165        DO K = 2, Nr        DO k = 2, Nr
166         Km1 = K-1         km1 = k-1
167  c      Kp1 = MIN(Nr,K+1)  c      kp1 = MIN(Nr,k+1)
168         CALL FIND_RHO_2D(         CALL FIND_RHO_2D(
169       I      iMin, iMax, jMin, jMax, K,       I      iMin, iMax, jMin, jMax, k,
170       I      theta(1-OLx,1-OLy,Km1,bi,bj), salt(1-OLx,1-OLy,Km1,bi,bj),       I      theta(1-OLx,1-OLy,km1,bi,bj), salt(1-OLx,1-OLy,km1,bi,bj),
171       O      rhoKm1,       O      rhoKm1,
172       I      Km1, bi, bj, myThid )       I      km1, bi, bj, myThid )
173    
174         CALL FIND_RHO_2D(         CALL FIND_RHO_2D(
175       I      iMin, iMax, jMin, jMax, K,       I      iMin, iMax, jMin, jMax, k,
176       I      theta(1-OLx,1-OLy,K,bi,bj), salt(1-OLx,1-OLy,K,bi,bj),       I      theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
177       O      rhoK,       O      rhoK,
178       I      K, bi, bj, myThid )       I      k, bi, bj, myThid )
179         DO J=jMin,jMax         DO j=jMin,jMax
180          DO I=iMin,iMax          DO i=iMin,iMax
181           SQRTTKE(i,j,k)=SQRT( GGL90TKE(I,J,K,bi,bj) )           SQRTTKE(i,j,k)=SQRT( GGL90TKE(i,j,k,bi,bj) )
182    
183  C     buoyancy frequency  C     buoyancy frequency
184           Nsquare(i,j,k) = - gravity*recip_rhoConst*recip_drC(K)           Nsquare(i,j,k) = gravity*gravitySign*recip_rhoConst
185       &        * ( rhoKm1(I,J) - rhoK(I,J) )*maskC(I,J,K,bi,bj)       &                  * sigmaR(i,j,k)
186             Nsquare(i,j,k) = - gravity*recip_rhoConst*recip_drC(k)
187         &        * ( rhoKm1(i,j) - rhoK(i,j) )*maskC(i,j,k,bi,bj)
188  cC     vertical shear term (dU/dz)^2+(dV/dz)^2  cC     vertical shear term (dU/dz)^2+(dV/dz)^2
189  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)
190  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)) )
191  c     &        *recip_drC(K)  c     &        *recip_drC(k)
192  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)
193  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)) )
194  c     &        *recip_drC(K)  c     &        *recip_drC(k)
195  c         verticalShear = tempU*tempU + tempV*tempV  c         verticalShear = tempU*tempU + tempV*tempV
196  c         RiNumber   = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)  c         RiNumber   = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)
197  cC     compute Prandtl number (always greater than 0)  cC     compute Prandtl number (always greater than 0)
198  c         prTemp = 1. _d 0  c         prTemp = 1. _d 0
199  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
200  c         TKEPrandtlNumber(I,J,K) = MIN(10. _d 0,prTemp)  c         TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
201  C     mixing length  C     mixing length
202           GGL90mixingLength(I,J,K) = SQRTTWO *           GGL90mixingLength(i,j,k) = SQRTTWO *
203       &        SQRTTKE(i,j,k)/SQRT( MAX(Nsquare(i,j,k),GGL90eps) )       &        SQRTTKE(i,j,k)/SQRT( MAX(Nsquare(i,j,k),GGL90eps) )
204          ENDDO          ENDDO
205         ENDDO         ENDDO
# Line 196  C     mixing length Line 207  C     mixing length
207    
208  C- Impose upper and lower bound for mixing length  C- Impose upper and lower bound for mixing length
209        IF ( mxlMaxFlag .EQ. 0 ) THEN        IF ( mxlMaxFlag .EQ. 0 ) THEN
210  C-  
211         DO k=2,Nr         DO k=2,Nr
212          DO J=jMin,jMax          DO j=jMin,jMax
213           DO I=iMin,iMax           DO i=iMin,iMax
214            MaxLength=totalDepth(I,J)            MaxLength=totalDepth(i,j)
215            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
216       &                                   MaxLength)       &                                   MaxLength)
217           ENDDO           ENDDO
218          ENDDO          ENDDO
219         ENDDO         ENDDO
220    
221         DO k=2,Nr         DO k=2,Nr
222          DO J=jMin,jMax          DO j=jMin,jMax
223           DO I=iMin,iMax           DO i=iMin,iMax
224            GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
225       &                                   GGL90mixingLengthMin)       &                                   GGL90mixingLengthMin)
226            rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K)            rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
227           ENDDO           ENDDO
228          ENDDO          ENDDO
229         ENDDO         ENDDO
230    
231        ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN
232  C-  
233         DO k=2,Nr         DO k=2,Nr
234          DO J=jMin,jMax          DO j=jMin,jMax
235           DO I=iMin,iMax           DO i=iMin,iMax
236            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))
237  c         MaxLength=MAX(MaxLength,20. _d 0)  c         MaxLength=MAX(MaxLength,20. _d 0)
238            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
239       &                                   MaxLength)       &                                   MaxLength)
240           ENDDO           ENDDO
241          ENDDO          ENDDO
242         ENDDO         ENDDO
243    
244         DO k=2,Nr         DO k=2,Nr
245          DO J=jMin,jMax          DO j=jMin,jMax
246           DO I=iMin,iMax           DO i=iMin,iMax
247            GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
248       &                                   GGL90mixingLengthMin)       &                                   GGL90mixingLengthMin)
249            rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K)            rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
250           ENDDO           ENDDO
251          ENDDO          ENDDO
252         ENDDO         ENDDO
253    
254        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN
255  C-  
256  cgf ensure mixing between first and second level  cgf ensure mixing between first and second level
257  c      DO J=jMin,jMax  c      DO j=jMin,jMax
258  c        DO I=iMin,iMax  c        DO i=iMin,iMax
259  c         GGL90mixingLength(I,J,2)=drF(1)  c         GGL90mixingLength(i,j,2)=drF(1)
260  c        ENDDO  c        ENDDO
261  c      ENDDO  c      ENDDO
262  cgf  cgf
263         DO k=2,Nr         DO k=2,Nr
264          DO J=jMin,jMax          DO j=jMin,jMax
265           DO I=iMin,iMax           DO i=iMin,iMax
266            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
267       &        GGL90mixingLength(I,J,K-1)+drF(k-1))       &        GGL90mixingLength(i,j,k-1)+drF(k-1))
268           ENDDO           ENDDO
269          ENDDO          ENDDO
270         ENDDO         ENDDO
271         DO J=jMin,jMax         DO j=jMin,jMax
272          DO I=iMin,iMax          DO i=iMin,iMax
273            GGL90mixingLength(I,J,Nr) = MIN(GGL90mixingLength(I,J,Nr),            GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
274       &       GGL90mixingLengthMin+drF(Nr))       &       GGL90mixingLengthMin+drF(Nr))
275          ENDDO          ENDDO
276         ENDDO         ENDDO
277         DO k=Nr-1,2,-1         DO k=Nr-1,2,-1
278          DO J=jMin,jMax          DO j=jMin,jMax
279           DO I=iMin,iMax           DO i=iMin,iMax
280            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
281       &        GGL90mixingLength(I,J,K+1)+drF(k))       &        GGL90mixingLength(i,j,k+1)+drF(k))
282           ENDDO           ENDDO
283          ENDDO          ENDDO
284         ENDDO         ENDDO
285    
286         DO k=2,Nr         DO k=2,Nr
287          DO J=jMin,jMax          DO j=jMin,jMax
288           DO I=iMin,iMax           DO i=iMin,iMax
289            GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
290       &                                   GGL90mixingLengthMin)       &                                   GGL90mixingLengthMin)
291            rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K)            rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
292           ENDDO           ENDDO
293          ENDDO          ENDDO
294         ENDDO         ENDDO
295    
296        ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN
297  C-  
298         DO k=2,Nr         DO k=2,Nr
299          DO J=jMin,jMax          DO j=jMin,jMax
300           DO I=iMin,iMax           DO i=iMin,iMax
301            mxLength_Dn(I,J,K) = MIN(GGL90mixingLength(I,J,K),            mxLength_Dn(i,j,k) = MIN(GGL90mixingLength(i,j,k),
302       &        mxLength_Dn(I,J,K-1)+drF(k-1))       &        mxLength_Dn(i,j,k-1)+drF(k-1))
303           ENDDO           ENDDO
304          ENDDO          ENDDO
305         ENDDO         ENDDO
306         DO J=jMin,jMax         DO j=jMin,jMax
307          DO I=iMin,iMax          DO i=iMin,iMax
308            GGL90mixingLength(I,J,Nr) = MIN(GGL90mixingLength(I,J,Nr),            GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
309       &       GGL90mixingLengthMin+drF(Nr))       &       GGL90mixingLengthMin+drF(Nr))
310          ENDDO          ENDDO
311         ENDDO         ENDDO
312         DO k=Nr-1,2,-1         DO k=Nr-1,2,-1
313          DO J=jMin,jMax          DO j=jMin,jMax
314           DO I=iMin,iMax           DO i=iMin,iMax
315            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
316       &        GGL90mixingLength(I,J,K+1)+drF(k))       &        GGL90mixingLength(i,j,k+1)+drF(k))
317           ENDDO           ENDDO
318          ENDDO          ENDDO
319         ENDDO         ENDDO
320    
321         DO k=2,Nr         DO k=2,Nr
322          DO J=jMin,jMax          DO j=jMin,jMax
323           DO I=iMin,iMax           DO i=iMin,iMax
324            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
325       &                                  mxLength_Dn(I,J,K))       &                                  mxLength_Dn(i,j,k))
326            tmpmlx = SQRT( GGL90mixingLength(I,J,K)*mxLength_Dn(I,J,K) )            tmpmlx = SQRT( GGL90mixingLength(i,j,k)*mxLength_Dn(i,j,k) )
327            tmpmlx = MAX( tmpmlx, GGL90mixingLengthMin)            tmpmlx = MAX( tmpmlx, GGL90mixingLengthMin)
328            rMixingLength(I,J,K) = 1. _d 0 / tmpmlx            rMixingLength(i,j,k) = 1. _d 0 / tmpmlx
329           ENDDO           ENDDO
330          ENDDO          ENDDO
331         ENDDO         ENDDO
# Line 325  C- Line 336  C-
336    
337  C- Impose minimum mixing length (to avoid division by zero)  C- Impose minimum mixing length (to avoid division by zero)
338  c      DO k=2,Nr  c      DO k=2,Nr
339  c      DO J=jMin,jMax  c      DO j=jMin,jMax
340  c       DO I=iMin,iMax  c       DO i=iMin,iMax
341  c        GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),  c        GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
342  c    &        GGL90mixingLengthMin)  c    &        GGL90mixingLengthMin)
343  c        rMixingLength(I,J,K) = 1. _d 0 /GGL90mixingLength(I,J,K)  c        rMixingLength(i,j,k) = 1. _d 0 /GGL90mixingLength(i,j,k)
344  c       ENDDO  c       ENDDO
345  c      ENDDO  c      ENDDO
346  c     ENDDO  c     ENDDO
347    
   
348        DO k=2,Nr        DO k=2,Nr
349         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  
350    
351    #ifdef ALLOW_GGL90_HORIZDIFF
352           IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
353  C     horizontal diffusion of TKE (requires an exchange in  C     horizontal diffusion of TKE (requires an exchange in
354  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  
355  C     common factors  C     common factors
356          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
357           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
358            xA(i,j) = _dyG(i,j,bi,bj)            xA(i,j) = _dyG(i,j,bi,bj)
359       &         *drF(k)*_hFacW(i,j,k,bi,bj)       &         *drF(k)*_hFacW(i,j,k,bi,bj)
360            yA(i,j) = _dxG(i,j,bi,bj)            yA(i,j) = _dxG(i,j,bi,bj)
# Line 395  C     common factors Line 363  C     common factors
363          ENDDO          ENDDO
364  C     Compute diffusive fluxes  C     Compute diffusive fluxes
365  C     ... across x-faces  C     ... across x-faces
366          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
367           dfx(1-Olx,j)=0. _d 0           dfx(1-OLx,j)=0. _d 0
368           DO i=1-Olx+1,sNx+Olx           DO i=1-OLx+1,sNx+OLx
369            dfx(i,j) = -GGL90diffTKEh*xA(i,j)            dfx(i,j) = -GGL90diffTKEh*xA(i,j)
370       &      *_recip_dxC(i,j,bi,bj)       &      *_recip_dxC(i,j,bi,bj)
371       &      *(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 373  C     ... across x-faces
373           ENDDO           ENDDO
374          ENDDO          ENDDO
375  C     ... across y-faces  C     ... across y-faces
376          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
377           dfy(i,1-Oly)=0. _d 0           dfy(i,1-OLy)=0. _d 0
378          ENDDO          ENDDO
379          DO j=1-Oly+1,sNy+Oly          DO j=1-OLy+1,sNy+OLy
380           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
381            dfy(i,j) = -GGL90diffTKEh*yA(i,j)            dfy(i,j) = -GGL90diffTKEh*yA(i,j)
382       &      *_recip_dyC(i,j,bi,bj)       &      *_recip_dyC(i,j,bi,bj)
383       &      *(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 387  C     ... across y-faces
387           ENDDO           ENDDO
388          ENDDO          ENDDO
389  C     Compute divergence of fluxes  C     Compute divergence of fluxes
390          DO j=1-Oly,sNy+Oly-1          DO j=1-OLy,sNy+OLy-1
391           DO i=1-Olx,sNx+Olx-1           DO i=1-OLx,sNx+OLx-1
392            gTKE(i,j,k)=gTKE(i,j,k)            gTKE(i,j) =
393       &   -_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)
394       &         *( (dfx(i+1,j)-dfx(i,j))       &         *( (dfx(i+1,j)-dfx(i,j))
395       &           +(dfy(i,j+1)-dfy(i,j))       &           +(dfy(i,j+1)-dfy(i,j))
396       &           )*deltaTggl90       &          )
397           ENDDO           ENDDO
398          ENDDO          ENDDO
399  C       end of k-loop  C      end if GGL90diffTKEh .eq. 0.
400           ENDIF
401    #endif /* ALLOW_GGL90_HORIZDIFF */
402    
403           DO j=jMin,jMax
404            DO i=iMin,iMax
405    C     vertical shear term (dU/dz)^2+(dV/dz)^2
406             tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)
407         &                 -( uVel(i,j,k  ,bi,bj)+uVel(i+1,j,k  ,bi,bj)) )
408         &        *recip_drC(k)
409             tempV= .5 _d 0*( vVel(i,j,km1,bi,bj)+vVel(i,j+1,km1,bi,bj)
410         &                 -( vVel(i,j,k  ,bi,bj)+vVel(i,j+1,k  ,bi,bj)) )
411         &        *recip_drC(k)
412             verticalShear = tempU*tempU + tempV*tempV
413             RiNumber = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)
414    C     compute Prandtl number (always greater than 0)
415             prTemp = 1. _d 0
416             IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
417             TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
418    c         TKEPrandtlNumber(i,j,k) = 1. _d 0
419    
420    C     viscosity and diffusivity
421             KappaM = GGL90ck*GGL90mixingLength(i,j,k)*SQRTTKE(i,j,k)
422             GGL90visctmp(i,j,k) = MAX(KappaM,diffKrNrT(k))
423         &                            * maskC(i,j,k,bi,bj)
424    c        note: storing GGL90visctmp like this, and using it later to compute
425    c              GGL9rdiffKr etc. is robust in case of smoothing (e.g. see OPA)
426             KappaM = MAX(KappaM,viscArNr(k)) * maskC(i,j,k,bi,bj)
427             KappaH = KappaM/TKEPrandtlNumber(i,j,k)
428             KappaE(i,j,k) = GGL90alpha * KappaM * maskC(i,j,k,bi,bj)
429    
430    C     dissipation term
431             TKEdissipation = explDissFac*GGL90ceps
432         &        *SQRTTKE(i,j,k)*rMixingLength(i,j,k)
433         &        *GGL90TKE(i,j,k,bi,bj)
434    C     partial update with sum of explicit contributions
435             GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
436         &        + deltaTggl90*(
437         &        + KappaM*verticalShear
438         &        - KappaH*Nsquare(i,j,k)
439         &        - TKEdissipation
440         &        )
441            ENDDO
442         ENDDO         ENDDO
443  C     end if GGL90diffTKEh .eq. 0.  
444        ENDIF  #ifdef ALLOW_GGL90_HORIZDIFF
445           IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
446    C--    Add horiz. diffusion tendency
447            DO j=jMin,jMax
448             DO i=iMin,iMax
449              GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
450         &                          + gTKE(i,j)*deltaTggl90
451             ENDDO
452            ENDDO
453           ENDIF
454  #endif /* ALLOW_GGL90_HORIZDIFF */  #endif /* ALLOW_GGL90_HORIZDIFF */
455    
456    C--   end of k loop
457          ENDDO
458    
459  C     ============================================  C     ============================================
460  C     Implicit time step to update TKE for k=1,Nr;  C     Implicit time step to update TKE for k=1,Nr;
461  C     TKE(Nr+1)=0 by default  C     TKE(Nr+1)=0 by default
# Line 442  C     set up matrix Line 464  C     set up matrix
464  C--   Lower diagonal  C--   Lower diagonal
465        DO j=jMin,jMax        DO j=jMin,jMax
466         DO i=iMin,iMax         DO i=iMin,iMax
467           a(i,j,1) = 0. _d 0           a3d(i,j,1) = 0. _d 0
468         ENDDO         ENDDO
469        ENDDO        ENDDO
470        DO k=2,Nr        DO k=2,Nr
# Line 452  C--   Lower diagonal Line 474  C--   Lower diagonal
474  C-    We keep recip_hFacC in the diffusive flux calculation,  C-    We keep recip_hFacC in the diffusive flux calculation,
475  C-    but no hFacC in TKE volume control  C-    but no hFacC in TKE volume control
476  C-    No need for maskC(k-1) with recip_hFacC(k-1)  C-    No need for maskC(k-1) with recip_hFacC(k-1)
477           a(i,j,k) = -deltaTggl90           a3d(i,j,k) = -deltaTggl90
478       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)
479       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))
480       &        *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 484  C-    No need for maskC(k-1) with recip_
484  C--   Upper diagonal  C--   Upper diagonal
485        DO j=jMin,jMax        DO j=jMin,jMax
486         DO i=iMin,iMax         DO i=iMin,iMax
487           c(i,j,1)  = 0. _d 0           c3d(i,j,1)  = 0. _d 0
488         ENDDO         ENDDO
489        ENDDO        ENDDO
490        DO k=2,Nr        DO k=2,Nr
# Line 472  C--   Upper diagonal Line 494  C--   Upper diagonal
494  C-    We keep recip_hFacC in the diffusive flux calculation,  C-    We keep recip_hFacC in the diffusive flux calculation,
495  C-    but no hFacC in TKE volume control  C-    but no hFacC in TKE volume control
496  C-    No need for maskC(k) with recip_hFacC(k)  C-    No need for maskC(k) with recip_hFacC(k)
497            c(i,j,k) = -deltaTggl90            c3d(i,j,k) = -deltaTggl90
498       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)
499       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))
500       &        *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 506  C--   Center diagonal
506         km1 = MAX(k-1,1)         km1 = MAX(k-1,1)
507         DO j=jMin,jMax         DO j=jMin,jMax
508          DO i=iMin,iMax          DO i=iMin,iMax
509            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)
510       &        + ab15*deltaTggl90*GGL90ceps*SQRTTKE(I,J,K)       &        + implDissFac*deltaTggl90*GGL90ceps*SQRTTKE(i,j,k)
511       &        * rMixingLength(I,J,K)       &        * rMixingLength(i,j,k)
512       &        * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)       &        * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)
513           ENDDO           ENDDO
514         ENDDO         ENDDO
# Line 495  C     end set up matrix Line 517  C     end set up matrix
517    
518  C     Apply boundary condition  C     Apply boundary condition
519        kp1 = MIN(Nr,kSurf+1)        kp1 = MIN(Nr,kSurf+1)
520        DO J=jMin,jMax        DO j=jMin,jMax
521         DO I=iMin,iMax         DO i=iMin,iMax
522  C     estimate friction velocity uStar from surface forcing  C     estimate friction velocity uStar from surface forcing
523          uStarSquare = SQRT(          uStarSquare = SQRT(
524       &    ( .5 _d 0*( surfaceForcingU(I,  J,  bi,bj)       &    ( .5 _d 0*( surfaceForcingU(i,  j,  bi,bj)
525       &              + surfaceForcingU(I+1,J,  bi,bj) ) )**2       &              + surfaceForcingU(i+1,j,  bi,bj) ) )**2
526       &  + ( .5 _d 0*( surfaceForcingV(I,  J,  bi,bj)       &  + ( .5 _d 0*( surfaceForcingV(i,  j,  bi,bj)
527       &              + surfaceForcingV(I,  J+1,bi,bj) ) )**2       &              + surfaceForcingV(i,  j+1,bi,bj) ) )**2
528       &                     )       &                     )
529  C     Dirichlet surface boundary condition for TKE  C     Dirichlet surface boundary condition for TKE
530          gTKE(I,J,kSurf) = MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)          GGL90TKE(i,j,kSurf,bi,bj) = maskC(i,j,kSurf,bi,bj)
531       &                     *maskC(I,J,kSurf,bi,bj)       &           *MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)
532          gTKE(i,j,kp1) = gTKE(i,j,kp1)          GGL90TKE(i,j,kp1,bi,bj) = GGL90TKE(i,j,kp1,bi,bj)
533       &                - a(i,j,kp1)*gTKE(i,j,kSurf)       &               - a3d(i,j,kp1)*GGL90TKE(i,j,kSurf,bi,bj)
534          a(i,j,kp1) = 0. _d 0          a3d(i,j,kp1) = 0. _d 0
535  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom
536          kBottom   = MAX(kLowC(I,J,bi,bj),1)          kBottom   = MAX(kLowC(i,j,bi,bj),1)
537          gTKE(I,J,kBottom) = gTKE(I,J,kBottom)          GGL90TKE(i,j,kBottom,bi,bj) = GGL90TKE(i,j,kBottom,bi,bj)
538       &                    - GGL90TKEbottom*c(I,J,kBottom)       &                              - GGL90TKEbottom*c3d(i,j,kBottom)
539          c(I,J,kBottom) = 0. _d 0          c3d(i,j,kBottom) = 0. _d 0
540         ENDDO         ENDDO
541        ENDDO        ENDDO
542    
543  C     solve tri-diagonal system, and store solution on gTKE (previously rhs)  C     solve tri-diagonal system
544        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
545       I                        a, b, c,       I                        a3d, b3d, c3d,
546       U                        gTKE,       U                        GGL90TKE,
547       O                        errCode,       O                        errCode,
548       I                        bi, bj, myThid )       I                        bi, bj, myThid )
549    
550  C     now update TKE        DO k=1,Nr
551        DO K=1,Nr         DO j=jMin,jMax
552         DO J=jMin,jMax          DO i=iMin,iMax
         DO I=iMin,iMax  
553  C     impose minimum TKE to avoid numerical undershoots below zero  C     impose minimum TKE to avoid numerical undershoots below zero
554           GGL90TKE(I,J,K,bi,bj) = MAX( gTKE(I,J,K), GGL90TKEmin )           GGL90TKE(i,j,k,bi,bj) = maskC(i,j,k,bi,bj)
555       &        * maskC(I,J,K,bi,bj)       &                  *MAX( GGL90TKE(i,j,k,bi,bj), GGL90TKEmin )
556          ENDDO          ENDDO
557         ENDDO         ENDDO
558        ENDDO        ENDDO
# Line 539  C     impose minimum TKE to avoid numeri Line 560  C     impose minimum TKE to avoid numeri
560  C     end of time step  C     end of time step
561  C     ===============================  C     ===============================
562    
563        DO K=2,Nr        DO k=2,Nr
564         DO J=1,sNy         DO j=1,sNy
565          DO I=1,sNx          DO i=1,sNx
566  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
567           tmpVisc=           tmpVisc=
568       &  (       &  (
# Line 566  C     =============================== Line 587  C     ===============================
587       &       +       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))
588       &  )*maskC(i,j,k,bi,bj)*mskCor(i,j,bi,bj)       &  )*maskC(i,j,k,bi,bj)*mskCor(i,j,bi,bj)
589  #else  #else
590           tmpVisc = GGL90visctmp(I,J,K)           tmpVisc = GGL90visctmp(i,j,k)
591  #endif  #endif
592           tmpVisc = MIN(tmpVisc/TKEPrandtlNumber(i,j,k),GGL90diffMax)           tmpVisc = MIN(tmpVisc/TKEPrandtlNumber(i,j,k),GGL90diffMax)
593           GGL90diffKr(I,J,K,bi,bj)= MAX( tmpVisc , diffKrNrT(k) )           GGL90diffKr(i,j,k,bi,bj)= MAX( tmpVisc , diffKrNrT(k) )
594          ENDDO          ENDDO
595         ENDDO         ENDDO
596        ENDDO        ENDDO
597    
598          DO k=2,Nr
599           DO j=1,sNy
600        DO K=2,Nr          DO i=1,sNx+1
        DO J=1,sNy  
         DO I=1,sNx+1  
601  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
602          tmpVisc =          tmpVisc =
603       & (       & (
604       &   p4 *(GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)       &   p4 *(GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
605       &       +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 622  C     ===============================
622       &                            +GGL90visctmp(i-1,j,k))       &                            +GGL90visctmp(i-1,j,k))
623       &                   )       &                   )
624  #endif  #endif
625          tmpVisc = MIN( tmpVisc , GGL90viscMax )          tmpVisc = MIN( tmpVisc , GGL90viscMax )
626          GGL90viscArU(i,j,k,bi,bj) = MAX( tmpVisc , viscArNr(k)  )          GGL90viscArU(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
627          ENDDO          ENDDO
628         ENDDO         ENDDO
629        ENDDO        ENDDO
630    
631          DO k=2,Nr
632        DO K=2,Nr         DO j=1,sNy+1
633         DO J=1,sNy+1          DO i=1,sNx
         DO I=1,sNx  
634  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
635          tmpVisc =          tmpVisc =
636       & (       & (
# Line 639  C     =============================== Line 657  C     ===============================
657    
658  #endif  #endif
659          tmpVisc = MIN( tmpVisc , GGL90viscMax )          tmpVisc = MIN( tmpVisc , GGL90viscMax )
660          GGL90viscArV(i,j,k,bi,bj) = MAX( tmpVisc , viscArNr(k)  )          GGL90viscArV(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
661          ENDDO          ENDDO
662         ENDDO         ENDDO
663        ENDDO        ENDDO
Legend:
Removed from v.1.17  
changed lines
  Added in v.1.21
  ViewVC Help
Powered by ViewVC 1.1.22