/[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.22 by jmc, Thu Jun 28 15:35:52 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
# Line 196  C     mixing length Line 189  C     mixing length
189    
190  C- Impose upper and lower bound for mixing length  C- Impose upper and lower bound for mixing length
191        IF ( mxlMaxFlag .EQ. 0 ) THEN        IF ( mxlMaxFlag .EQ. 0 ) THEN
192  C-  
193         DO k=2,Nr         DO k=2,Nr
194          DO J=jMin,jMax          DO j=jMin,jMax
195           DO I=iMin,iMax           DO i=iMin,iMax
196            MaxLength=totalDepth(I,J)            MaxLength=totalDepth(i,j)
197            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
198       &                                   MaxLength)       &                                   MaxLength)
199           ENDDO           ENDDO
200          ENDDO          ENDDO
201         ENDDO         ENDDO
202    
203         DO k=2,Nr         DO k=2,Nr
204          DO J=jMin,jMax          DO j=jMin,jMax
205           DO I=iMin,iMax           DO i=iMin,iMax
206            GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
207       &                                   GGL90mixingLengthMin)       &                                   GGL90mixingLengthMin)
208            rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K)            rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
209           ENDDO           ENDDO
210          ENDDO          ENDDO
211         ENDDO         ENDDO
212    
213        ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN
214  C-  
215         DO k=2,Nr         DO k=2,Nr
216          DO J=jMin,jMax          DO j=jMin,jMax
217           DO I=iMin,iMax           DO i=iMin,iMax
218            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))
219  c         MaxLength=MAX(MaxLength,20. _d 0)  c         MaxLength=MAX(MaxLength,20. _d 0)
220            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
221       &                                   MaxLength)       &                                   MaxLength)
222           ENDDO           ENDDO
223          ENDDO          ENDDO
224         ENDDO         ENDDO
225    
226         DO k=2,Nr         DO k=2,Nr
227          DO J=jMin,jMax          DO j=jMin,jMax
228           DO I=iMin,iMax           DO i=iMin,iMax
229            GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
230       &                                   GGL90mixingLengthMin)       &                                   GGL90mixingLengthMin)
231            rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K)            rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
232           ENDDO           ENDDO
233          ENDDO          ENDDO
234         ENDDO         ENDDO
235    
236        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN
237  C-  
238  cgf ensure mixing between first and second level  cgf ensure mixing between first and second level
239  c      DO J=jMin,jMax  c      DO j=jMin,jMax
240  c        DO I=iMin,iMax  c        DO i=iMin,iMax
241  c         GGL90mixingLength(I,J,2)=drF(1)  c         GGL90mixingLength(i,j,2)=drF(1)
242  c        ENDDO  c        ENDDO
243  c      ENDDO  c      ENDDO
244  cgf  cgf
245         DO k=2,Nr         DO k=2,Nr
246          DO J=jMin,jMax          DO j=jMin,jMax
247           DO I=iMin,iMax           DO i=iMin,iMax
248            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
249       &        GGL90mixingLength(I,J,K-1)+drF(k-1))       &        GGL90mixingLength(i,j,k-1)+drF(k-1))
250           ENDDO           ENDDO
251          ENDDO          ENDDO
252         ENDDO         ENDDO
253         DO J=jMin,jMax         DO j=jMin,jMax
254          DO I=iMin,iMax          DO i=iMin,iMax
255            GGL90mixingLength(I,J,Nr) = MIN(GGL90mixingLength(I,J,Nr),            GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
256       &       GGL90mixingLengthMin+drF(Nr))       &       GGL90mixingLengthMin+drF(Nr))
257          ENDDO          ENDDO
258         ENDDO         ENDDO
259         DO k=Nr-1,2,-1         DO k=Nr-1,2,-1
260          DO J=jMin,jMax          DO j=jMin,jMax
261           DO I=iMin,iMax           DO i=iMin,iMax
262            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
263       &        GGL90mixingLength(I,J,K+1)+drF(k))       &        GGL90mixingLength(i,j,k+1)+drF(k))
264           ENDDO           ENDDO
265          ENDDO          ENDDO
266         ENDDO         ENDDO
267    
268         DO k=2,Nr         DO k=2,Nr
269          DO J=jMin,jMax          DO j=jMin,jMax
270           DO I=iMin,iMax           DO i=iMin,iMax
271            GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
272       &                                   GGL90mixingLengthMin)       &                                   GGL90mixingLengthMin)
273            rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K)            rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
274           ENDDO           ENDDO
275          ENDDO          ENDDO
276         ENDDO         ENDDO
277    
278        ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN
279  C-  
280         DO k=2,Nr         DO k=2,Nr
281          DO J=jMin,jMax          DO j=jMin,jMax
282           DO I=iMin,iMax           DO i=iMin,iMax
283            mxLength_Dn(I,J,K) = MIN(GGL90mixingLength(I,J,K),            mxLength_Dn(i,j,k) = MIN(GGL90mixingLength(i,j,k),
284       &        mxLength_Dn(I,J,K-1)+drF(k-1))       &        mxLength_Dn(i,j,k-1)+drF(k-1))
285           ENDDO           ENDDO
286          ENDDO          ENDDO
287         ENDDO         ENDDO
288         DO J=jMin,jMax         DO j=jMin,jMax
289          DO I=iMin,iMax          DO i=iMin,iMax
290            GGL90mixingLength(I,J,Nr) = MIN(GGL90mixingLength(I,J,Nr),            GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
291       &       GGL90mixingLengthMin+drF(Nr))       &       GGL90mixingLengthMin+drF(Nr))
292          ENDDO          ENDDO
293         ENDDO         ENDDO
294         DO k=Nr-1,2,-1         DO k=Nr-1,2,-1
295          DO J=jMin,jMax          DO j=jMin,jMax
296           DO I=iMin,iMax           DO i=iMin,iMax
297            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
298       &        GGL90mixingLength(I,J,K+1)+drF(k))       &        GGL90mixingLength(i,j,k+1)+drF(k))
299           ENDDO           ENDDO
300          ENDDO          ENDDO
301         ENDDO         ENDDO
302    
303         DO k=2,Nr         DO k=2,Nr
304          DO J=jMin,jMax          DO j=jMin,jMax
305           DO I=iMin,iMax           DO i=iMin,iMax
306            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
307       &                                  mxLength_Dn(I,J,K))       &                                  mxLength_Dn(i,j,k))
308            tmpmlx = SQRT( GGL90mixingLength(I,J,K)*mxLength_Dn(I,J,K) )            tmpmlx = SQRT( GGL90mixingLength(i,j,k)*mxLength_Dn(i,j,k) )
309            tmpmlx = MAX( tmpmlx, GGL90mixingLengthMin)            tmpmlx = MAX( tmpmlx, GGL90mixingLengthMin)
310            rMixingLength(I,J,K) = 1. _d 0 / tmpmlx            rMixingLength(i,j,k) = 1. _d 0 / tmpmlx
311           ENDDO           ENDDO
312          ENDDO          ENDDO
313         ENDDO         ENDDO
# Line 325  C- Line 318  C-
318    
319  C- Impose minimum mixing length (to avoid division by zero)  C- Impose minimum mixing length (to avoid division by zero)
320  c      DO k=2,Nr  c      DO k=2,Nr
321  c      DO J=jMin,jMax  c      DO j=jMin,jMax
322  c       DO I=iMin,iMax  c       DO i=iMin,iMax
323  c        GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),  c        GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
324  c    &        GGL90mixingLengthMin)  c    &        GGL90mixingLengthMin)
325  c        rMixingLength(I,J,K) = 1. _d 0 /GGL90mixingLength(I,J,K)  c        rMixingLength(i,j,k) = 1. _d 0 /GGL90mixingLength(i,j,k)
326  c       ENDDO  c       ENDDO
327  c      ENDDO  c      ENDDO
328  c     ENDDO  c     ENDDO
329    
   
330        DO k=2,Nr        DO k=2,Nr
331         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  
332    
333    #ifdef ALLOW_GGL90_HORIZDIFF
334           IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
335  C     horizontal diffusion of TKE (requires an exchange in  C     horizontal diffusion of TKE (requires an exchange in
336  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  
337  C     common factors  C     common factors
338          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
339           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
340            xA(i,j) = _dyG(i,j,bi,bj)            xA(i,j) = _dyG(i,j,bi,bj)
341       &         *drF(k)*_hFacW(i,j,k,bi,bj)       &         *drF(k)*_hFacW(i,j,k,bi,bj)
342            yA(i,j) = _dxG(i,j,bi,bj)            yA(i,j) = _dxG(i,j,bi,bj)
# Line 395  C     common factors Line 345  C     common factors
345          ENDDO          ENDDO
346  C     Compute diffusive fluxes  C     Compute diffusive fluxes
347  C     ... across x-faces  C     ... across x-faces
348          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
349           dfx(1-Olx,j)=0. _d 0           dfx(1-OLx,j)=0. _d 0
350           DO i=1-Olx+1,sNx+Olx           DO i=1-OLx+1,sNx+OLx
351            dfx(i,j) = -GGL90diffTKEh*xA(i,j)            dfx(i,j) = -GGL90diffTKEh*xA(i,j)
352       &      *_recip_dxC(i,j,bi,bj)       &      *_recip_dxC(i,j,bi,bj)
353       &      *(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 355  C     ... across x-faces
355           ENDDO           ENDDO
356          ENDDO          ENDDO
357  C     ... across y-faces  C     ... across y-faces
358          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
359           dfy(i,1-Oly)=0. _d 0           dfy(i,1-OLy)=0. _d 0
360          ENDDO          ENDDO
361          DO j=1-Oly+1,sNy+Oly          DO j=1-OLy+1,sNy+OLy
362           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
363            dfy(i,j) = -GGL90diffTKEh*yA(i,j)            dfy(i,j) = -GGL90diffTKEh*yA(i,j)
364       &      *_recip_dyC(i,j,bi,bj)       &      *_recip_dyC(i,j,bi,bj)
365       &      *(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 369  C     ... across y-faces
369           ENDDO           ENDDO
370          ENDDO          ENDDO
371  C     Compute divergence of fluxes  C     Compute divergence of fluxes
372          DO j=1-Oly,sNy+Oly-1          DO j=1-OLy,sNy+OLy-1
373           DO i=1-Olx,sNx+Olx-1           DO i=1-OLx,sNx+OLx-1
374            gTKE(i,j,k)=gTKE(i,j,k)            gTKE(i,j) =
375       &   -_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)
376       &         *( (dfx(i+1,j)-dfx(i,j))       &         *( (dfx(i+1,j)-dfx(i,j))
377       &           +(dfy(i,j+1)-dfy(i,j))       &           +(dfy(i,j+1)-dfy(i,j))
378       &           )*deltaTggl90       &          )
379           ENDDO           ENDDO
380          ENDDO          ENDDO
381  C       end of k-loop  C      end if GGL90diffTKEh .eq. 0.
382           ENDIF
383    #endif /* ALLOW_GGL90_HORIZDIFF */
384    
385           DO j=jMin,jMax
386            DO i=iMin,iMax
387    C     vertical shear term (dU/dz)^2+(dV/dz)^2
388             tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)
389         &                 -( uVel(i,j,k  ,bi,bj)+uVel(i+1,j,k  ,bi,bj)) )
390         &        *recip_drC(k)
391             tempV= .5 _d 0*( vVel(i,j,km1,bi,bj)+vVel(i,j+1,km1,bi,bj)
392         &                 -( vVel(i,j,k  ,bi,bj)+vVel(i,j+1,k  ,bi,bj)) )
393         &        *recip_drC(k)
394             verticalShear = tempU*tempU + tempV*tempV
395             RiNumber = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)
396    C     compute Prandtl number (always greater than 0)
397             prTemp = 1. _d 0
398             IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
399             TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
400    c         TKEPrandtlNumber(i,j,k) = 1. _d 0
401    
402    C     viscosity and diffusivity
403             KappaM = GGL90ck*GGL90mixingLength(i,j,k)*SQRTTKE(i,j,k)
404             GGL90visctmp(i,j,k) = MAX(KappaM,diffKrNrT(k))
405         &                            * maskC(i,j,k,bi,bj)
406    c        note: storing GGL90visctmp like this, and using it later to compute
407    c              GGL9rdiffKr etc. is robust in case of smoothing (e.g. see OPA)
408             KappaM = MAX(KappaM,viscArNr(k)) * maskC(i,j,k,bi,bj)
409             KappaH = KappaM/TKEPrandtlNumber(i,j,k)
410             KappaE(i,j,k) = GGL90alpha * KappaM * maskC(i,j,k,bi,bj)
411    
412    C     dissipation term
413             TKEdissipation = explDissFac*GGL90ceps
414         &        *SQRTTKE(i,j,k)*rMixingLength(i,j,k)
415         &        *GGL90TKE(i,j,k,bi,bj)
416    C     partial update with sum of explicit contributions
417             GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
418         &        + deltaTggl90*(
419         &        + KappaM*verticalShear
420         &        - KappaH*Nsquare(i,j,k)
421         &        - TKEdissipation
422         &        )
423            ENDDO
424         ENDDO         ENDDO
425  C     end if GGL90diffTKEh .eq. 0.  
426        ENDIF  #ifdef ALLOW_GGL90_HORIZDIFF
427           IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
428    C--    Add horiz. diffusion tendency
429            DO j=jMin,jMax
430             DO i=iMin,iMax
431              GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
432         &                          + gTKE(i,j)*deltaTggl90
433             ENDDO
434            ENDDO
435           ENDIF
436  #endif /* ALLOW_GGL90_HORIZDIFF */  #endif /* ALLOW_GGL90_HORIZDIFF */
437    
438    C--   end of k loop
439          ENDDO
440    
441  C     ============================================  C     ============================================
442  C     Implicit time step to update TKE for k=1,Nr;  C     Implicit time step to update TKE for k=1,Nr;
443  C     TKE(Nr+1)=0 by default  C     TKE(Nr+1)=0 by default
# Line 442  C     set up matrix Line 446  C     set up matrix
446  C--   Lower diagonal  C--   Lower diagonal
447        DO j=jMin,jMax        DO j=jMin,jMax
448         DO i=iMin,iMax         DO i=iMin,iMax
449           a(i,j,1) = 0. _d 0           a3d(i,j,1) = 0. _d 0
450         ENDDO         ENDDO
451        ENDDO        ENDDO
452        DO k=2,Nr        DO k=2,Nr
# Line 452  C--   Lower diagonal Line 456  C--   Lower diagonal
456  C-    We keep recip_hFacC in the diffusive flux calculation,  C-    We keep recip_hFacC in the diffusive flux calculation,
457  C-    but no hFacC in TKE volume control  C-    but no hFacC in TKE volume control
458  C-    No need for maskC(k-1) with recip_hFacC(k-1)  C-    No need for maskC(k-1) with recip_hFacC(k-1)
459           a(i,j,k) = -deltaTggl90           a3d(i,j,k) = -deltaTggl90
460       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)
461       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))
462       &        *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 466  C-    No need for maskC(k-1) with recip_
466  C--   Upper diagonal  C--   Upper diagonal
467        DO j=jMin,jMax        DO j=jMin,jMax
468         DO i=iMin,iMax         DO i=iMin,iMax
469           c(i,j,1)  = 0. _d 0           c3d(i,j,1)  = 0. _d 0
470         ENDDO         ENDDO
471        ENDDO        ENDDO
472        DO k=2,Nr        DO k=2,Nr
# Line 472  C--   Upper diagonal Line 476  C--   Upper diagonal
476  C-    We keep recip_hFacC in the diffusive flux calculation,  C-    We keep recip_hFacC in the diffusive flux calculation,
477  C-    but no hFacC in TKE volume control  C-    but no hFacC in TKE volume control
478  C-    No need for maskC(k) with recip_hFacC(k)  C-    No need for maskC(k) with recip_hFacC(k)
479            c(i,j,k) = -deltaTggl90            c3d(i,j,k) = -deltaTggl90
480       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)
481       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))
482       &        *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 488  C--   Center diagonal
488         km1 = MAX(k-1,1)         km1 = MAX(k-1,1)
489         DO j=jMin,jMax         DO j=jMin,jMax
490          DO i=iMin,iMax          DO i=iMin,iMax
491            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)
492       &        + ab15*deltaTggl90*GGL90ceps*SQRTTKE(I,J,K)       &        + implDissFac*deltaTggl90*GGL90ceps*SQRTTKE(i,j,k)
493       &        * rMixingLength(I,J,K)       &        * rMixingLength(i,j,k)
494       &        * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)       &        * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)
495           ENDDO           ENDDO
496         ENDDO         ENDDO
# Line 495  C     end set up matrix Line 499  C     end set up matrix
499    
500  C     Apply boundary condition  C     Apply boundary condition
501        kp1 = MIN(Nr,kSurf+1)        kp1 = MIN(Nr,kSurf+1)
502        DO J=jMin,jMax        DO j=jMin,jMax
503         DO I=iMin,iMax         DO i=iMin,iMax
504  C     estimate friction velocity uStar from surface forcing  C     estimate friction velocity uStar from surface forcing
505          uStarSquare = SQRT(          uStarSquare = SQRT(
506       &    ( .5 _d 0*( surfaceForcingU(I,  J,  bi,bj)       &    ( .5 _d 0*( surfaceForcingU(i,  j,  bi,bj)
507       &              + surfaceForcingU(I+1,J,  bi,bj) ) )**2       &              + surfaceForcingU(i+1,j,  bi,bj) ) )**2
508       &  + ( .5 _d 0*( surfaceForcingV(I,  J,  bi,bj)       &  + ( .5 _d 0*( surfaceForcingV(i,  j,  bi,bj)
509       &              + surfaceForcingV(I,  J+1,bi,bj) ) )**2       &              + surfaceForcingV(i,  j+1,bi,bj) ) )**2
510       &                     )       &                     )
511  C     Dirichlet surface boundary condition for TKE  C     Dirichlet surface boundary condition for TKE
512          gTKE(I,J,kSurf) = MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)          GGL90TKE(i,j,kSurf,bi,bj) = maskC(i,j,kSurf,bi,bj)
513       &                     *maskC(I,J,kSurf,bi,bj)       &           *MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)
514          gTKE(i,j,kp1) = gTKE(i,j,kp1)          GGL90TKE(i,j,kp1,bi,bj) = GGL90TKE(i,j,kp1,bi,bj)
515       &                - a(i,j,kp1)*gTKE(i,j,kSurf)       &               - a3d(i,j,kp1)*GGL90TKE(i,j,kSurf,bi,bj)
516          a(i,j,kp1) = 0. _d 0          a3d(i,j,kp1) = 0. _d 0
517  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom
518          kBottom   = MAX(kLowC(I,J,bi,bj),1)          kBottom   = MAX(kLowC(i,j,bi,bj),1)
519          gTKE(I,J,kBottom) = gTKE(I,J,kBottom)          GGL90TKE(i,j,kBottom,bi,bj) = GGL90TKE(i,j,kBottom,bi,bj)
520       &                    - GGL90TKEbottom*c(I,J,kBottom)       &                              - GGL90TKEbottom*c3d(i,j,kBottom)
521          c(I,J,kBottom) = 0. _d 0          c3d(i,j,kBottom) = 0. _d 0
522         ENDDO         ENDDO
523        ENDDO        ENDDO
524    
525  C     solve tri-diagonal system, and store solution on gTKE (previously rhs)  C     solve tri-diagonal system
526        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
527       I                        a, b, c,       I                        a3d, b3d, c3d,
528       U                        gTKE,       U                        GGL90TKE,
529       O                        errCode,       O                        errCode,
530       I                        1, 1, myThid )       I                        bi, bj, myThid )
531    
532  C     now update TKE        DO k=1,Nr
533        DO K=1,Nr         DO j=jMin,jMax
534         DO J=jMin,jMax          DO i=iMin,iMax
         DO I=iMin,iMax  
535  C     impose minimum TKE to avoid numerical undershoots below zero  C     impose minimum TKE to avoid numerical undershoots below zero
536           GGL90TKE(I,J,K,bi,bj) = MAX( gTKE(I,J,K), GGL90TKEmin )           GGL90TKE(i,j,k,bi,bj) = maskC(i,j,k,bi,bj)
537       &        * maskC(I,J,K,bi,bj)       &                  *MAX( GGL90TKE(i,j,k,bi,bj), GGL90TKEmin )
538          ENDDO          ENDDO
539         ENDDO         ENDDO
540        ENDDO        ENDDO
# Line 539  C     impose minimum TKE to avoid numeri Line 542  C     impose minimum TKE to avoid numeri
542  C     end of time step  C     end of time step
543  C     ===============================  C     ===============================
544    
545        DO K=2,Nr        DO k=2,Nr
546         DO J=1,sNy         DO j=1,sNy
547          DO I=1,sNx          DO i=1,sNx
548  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
549           tmpVisc=           tmpVisc=
550       &  (       &  (
# Line 566  C     =============================== Line 569  C     ===============================
569       &       +       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))
570       &  )*maskC(i,j,k,bi,bj)*mskCor(i,j,bi,bj)       &  )*maskC(i,j,k,bi,bj)*mskCor(i,j,bi,bj)
571  #else  #else
572           tmpVisc = GGL90visctmp(I,J,K)           tmpVisc = GGL90visctmp(i,j,k)
573  #endif  #endif
574           tmpVisc = MIN(tmpVisc/TKEPrandtlNumber(i,j,k),GGL90diffMax)           tmpVisc = MIN(tmpVisc/TKEPrandtlNumber(i,j,k),GGL90diffMax)
575           GGL90diffKr(I,J,K,bi,bj)= MAX( tmpVisc , diffKrNrT(k) )           GGL90diffKr(i,j,k,bi,bj)= MAX( tmpVisc , diffKrNrT(k) )
576          ENDDO          ENDDO
577         ENDDO         ENDDO
578        ENDDO        ENDDO
579    
580          DO k=2,Nr
581           DO j=1,sNy
582        DO K=2,Nr          DO i=1,sNx+1
        DO J=1,sNy  
         DO I=1,sNx+1  
583  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
584          tmpVisc =          tmpVisc =
585       & (       & (
586       &   p4 *(GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)       &   p4 *(GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
587       &       +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 604  C     ===============================
604       &                            +GGL90visctmp(i-1,j,k))       &                            +GGL90visctmp(i-1,j,k))
605       &                   )       &                   )
606  #endif  #endif
607          tmpVisc = MIN( tmpVisc , GGL90viscMax )          tmpVisc = MIN( tmpVisc , GGL90viscMax )
608          GGL90viscArU(i,j,k,bi,bj) = MAX( tmpVisc , viscArNr(k)  )          GGL90viscArU(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
609          ENDDO          ENDDO
610         ENDDO         ENDDO
611        ENDDO        ENDDO
612    
613          DO k=2,Nr
614        DO K=2,Nr         DO j=1,sNy+1
615         DO J=1,sNy+1          DO i=1,sNx
         DO I=1,sNx  
616  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
617          tmpVisc =          tmpVisc =
618       & (       & (
# Line 639  C     =============================== Line 639  C     ===============================
639    
640  #endif  #endif
641          tmpVisc = MIN( tmpVisc , GGL90viscMax )          tmpVisc = MIN( tmpVisc , GGL90viscMax )
642          GGL90viscArV(i,j,k,bi,bj) = MAX( tmpVisc , viscArNr(k)  )          GGL90viscArV(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
643          ENDDO          ENDDO
644         ENDDO         ENDDO
645        ENDDO        ENDDO
Legend:
Removed from v.1.18  
changed lines
  Added in v.1.22
  ViewVC Help
Powered by ViewVC 1.1.22