/[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.15 by dfer, Sat Nov 14 14:27:56 2009 UTC revision 1.19 by jmc, Thu Aug 19 23:52:37 2010 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, myTime, myIter, myThid )
12    
13  C !DESCRIPTION: \bv  C !DESCRIPTION: \bv
14  C     *==========================================================*  C     *==========================================================*
# Line 43  C !INPUT PARAMETERS: =================== Line 43  C !INPUT PARAMETERS: ===================
43  C Routine arguments  C Routine arguments
44  C     bi, bj :: array indices on which to apply calculations  C     bi, bj :: array indices on which to apply calculations
45  C     myTime :: Current time in simulation  C     myTime :: Current time in simulation
46    C     myIter :: Current time-step number
47  C     myThid :: My Thread Id number  C     myThid :: My Thread Id number
48        INTEGER bi, bj        INTEGER bi, bj
49        _RL     myTime        _RL     myTime
50          INTEGER myIter
51        INTEGER myThid        INTEGER myThid
52  CEOP  CEOP
53    
# Line 53  CEOP Line 55  CEOP
55    
56  C !LOCAL VARIABLES: ====================================================  C !LOCAL VARIABLES: ====================================================
57  C Local constants  C Local constants
58  C     iMin, iMax, jMin, jMax, I, J - array computation indices  C     iMin,iMax,jMin,jMax :: index boundaries of computation domain
59  C     K, Kp1, km1, kSurf, kBottom  - vertical loop indices  C     i, j, k, kp1,km1 :: array computation indices
60  C     ab15, ab05       - weights for implicit timestepping  C     kSurf, kBottom   :: vertical indices of domain boundaries
61  C     uStarSquare      - square of friction velocity  C     explDissFac      :: explicit Dissipation Factor (in [0-1])
62  C     verticalShear    - (squared) vertical shear of horizontal velocity  C     implDissFac      :: implicit Dissipation Factor (in [0-1])
63  C     Nsquare          - squared buoyancy freqency  C     uStarSquare      :: square of friction velocity
64  C     RiNumber         - local Richardson number  C     verticalShear    :: (squared) vertical shear of horizontal velocity
65  C     KappaM           - (local) viscosity parameter (eq.10)  C     Nsquare          :: squared buoyancy freqency
66  C     KappaH           - (local) diffusivity parameter for temperature (eq.11)  C     RiNumber         :: local Richardson number
67  C     KappaE           - (local) diffusivity parameter for TKE (eq.15)  C     KappaM           :: (local) viscosity parameter (eq.10)
68  C     TKEdissipation   - dissipation of TKE  C     KappaH           :: (local) diffusivity parameter for temperature (eq.11)
69  C     GGL90mixingLength- mixing length of scheme following Banke+Delecuse  C     KappaE           :: (local) diffusivity parameter for TKE (eq.15)
70  C         rMixingLength- inverse of mixing length  C     TKEdissipation   :: dissipation of TKE
71  C     totalDepth       - thickness of water column (inverse of recip_Rcol)  C     GGL90mixingLength:: mixing length of scheme following Banke+Delecuse
72  C     TKEPrandtlNumber - here, an empirical function of the Richardson number  C         rMixingLength:: inverse of mixing length
73  C     rhoK, rhoKm1     - density at layer K and Km1 (relative to K)  C     totalDepth       :: thickness of water column (inverse of recip_Rcol)
74  C     gTKE             - right hand side of implicit equation  C     TKEPrandtlNumber :: here, an empirical function of the Richardson number
75    C     rhoK, rhoKm1     :: density at layer k and km1 (relative to k)
76        INTEGER iMin ,iMax ,jMin ,jMax        INTEGER iMin ,iMax ,jMin ,jMax
77        INTEGER I, J, K, Kp1, Km1, kSurf, kBottom        INTEGER i, j, k, kp1, km1, kSurf, kBottom
78        _RL     ab15, ab05        _RL     explDissFac, implDissFac
79        _RL     uStarSquare        _RL     uStarSquare
80        _RL     verticalShear        _RL     verticalShear
81        _RL     KappaM, KappaH        _RL     KappaM, KappaH
# Line 84  c     _RL     SQRTTKE Line 87  c     _RL     SQRTTKE
87        _RL     RiNumber        _RL     RiNumber
88        _RL     TKEdissipation        _RL     TKEdissipation
89        _RL     tempU, tempV, prTemp        _RL     tempU, tempV, prTemp
90        _RL     MaxLength, tmpmlx        _RL     MaxLength, tmpmlx, tmpVisc
91        _RL     TKEPrandtlNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     TKEPrandtlNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
92        _RL     GGL90mixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     GGL90mixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
93        _RL     rMixingLength    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     rMixingLength    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
# Line 93  c     _RL     SQRTTKE Line 96  c     _RL     SQRTTKE
96        _RL     rhoK             (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     rhoK             (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
97        _RL     rhoKm1           (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)
99        _RL     gTKE             (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     a(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
102        _RL     b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
103        _RL     c(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     c(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, tmpdiffKrS        _RL p4, p8, p16
117        p4=0.25 _d 0        p4=0.25 _d 0
118        p8=0.125 _d 0        p8=0.125 _d 0
119        p16=0.0625 _d 0        p16=0.0625 _d 0
# Line 122  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          ENDDO          ENDDO
143         ENDDO         ENDDO
144        ENDDO        ENDDO
145        DO J=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
146         DO I=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
147          rhoK(I,J)          = 0. _d 0          rhoK(i,j)          = 0. _d 0
148          rhoKm1(I,J)        = 0. _d 0          rhoKm1(i,j)        = 0. _d 0
149          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)
150          rMixingLength(i,j,1) = 0. _d 0          rMixingLength(i,j,1) = 0. _d 0
151          mxLength_Dn(I,J,1) = GGL90mixingLengthMin          mxLength_Dn(i,j,1) = GGL90mixingLengthMin
152          SQRTTKE(i,j,1) = SQRT( GGL90TKE(i,j,1,bi,bj) )          SQRTTKE(i,j,1) = SQRT( GGL90TKE(i,j,1,bi,bj) )
153         ENDDO         ENDDO
154        ENDDO        ENDDO
155    
156  C     start k-loop  C     start k-loop
157        DO K = 2, Nr        DO k = 2, Nr
158         Km1 = K-1         km1 = k-1
159  c      Kp1 = MIN(Nr,K+1)  c      kp1 = MIN(Nr,k+1)
160         CALL FIND_RHO_2D(         CALL FIND_RHO_2D(
161       I      iMin, iMax, jMin, jMax, K,       I      iMin, iMax, jMin, jMax, k,
162       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),
163       O      rhoKm1,       O      rhoKm1,
164       I      Km1, bi, bj, myThid )       I      km1, bi, bj, myThid )
165    
166         CALL FIND_RHO_2D(         CALL FIND_RHO_2D(
167       I      iMin, iMax, jMin, jMax, K,       I      iMin, iMax, jMin, jMax, k,
168       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),
169       O      rhoK,       O      rhoK,
170       I      K, bi, bj, myThid )       I      k, bi, bj, myThid )
171         DO J=jMin,jMax         DO j=jMin,jMax
172          DO I=iMin,iMax          DO i=iMin,iMax
173           SQRTTKE(i,j,k)=SQRT( GGL90TKE(I,J,K,bi,bj) )           SQRTTKE(i,j,k)=SQRT( GGL90TKE(i,j,k,bi,bj) )
174    
175  C     buoyancy frequency  C     buoyancy frequency
176           Nsquare(i,j,k) = - gravity*recip_rhoConst*recip_drC(K)           Nsquare(i,j,k) = - gravity*recip_rhoConst*recip_drC(k)
177       &        * ( rhoKm1(I,J) - rhoK(I,J) )*maskC(I,J,K,bi,bj)       &        * ( rhoKm1(i,j) - rhoK(i,j) )*maskC(i,j,k,bi,bj)
178  cC     vertical shear term (dU/dz)^2+(dV/dz)^2  cC     vertical shear term (dU/dz)^2+(dV/dz)^2
179  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)
180  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)) )
181  c     &        *recip_drC(K)  c     &        *recip_drC(k)
182  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)
183  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)) )
184  c     &        *recip_drC(K)  c     &        *recip_drC(k)
185  c         verticalShear = tempU*tempU + tempV*tempV  c         verticalShear = tempU*tempU + tempV*tempV
186  c         RiNumber   = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)  c         RiNumber   = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)
187  cC     compute Prandtl number (always greater than 0)  cC     compute Prandtl number (always greater than 0)
188  c         prTemp = 1. _d 0  c         prTemp = 1. _d 0
189  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
190  c         TKEPrandtlNumber(I,J,K) = MIN(10. _d 0,prTemp)  c         TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
191  C     mixing length  C     mixing length
192           GGL90mixingLength(I,J,K) = SQRTTWO *           GGL90mixingLength(i,j,k) = SQRTTWO *
193       &        SQRTTKE(i,j,k)/SQRT( MAX(Nsquare(i,j,k),GGL90eps) )       &        SQRTTKE(i,j,k)/SQRT( MAX(Nsquare(i,j,k),GGL90eps) )
194          ENDDO          ENDDO
195         ENDDO         ENDDO
# Line 196  C- Impose upper and lower bound for mixi Line 199  C- Impose upper and lower bound for mixi
199        IF ( mxlMaxFlag .EQ. 0 ) THEN        IF ( mxlMaxFlag .EQ. 0 ) THEN
200  C-  C-
201         DO k=2,Nr         DO k=2,Nr
202          DO J=jMin,jMax          DO j=jMin,jMax
203           DO I=iMin,iMax           DO i=iMin,iMax
204            MaxLength=totalDepth(I,J)            MaxLength=totalDepth(i,j)
205            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
206       &                                   MaxLength)       &                                   MaxLength)
207           ENDDO           ENDDO
208          ENDDO          ENDDO
209         ENDDO         ENDDO
210    
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            GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
215       &                                   GGL90mixingLengthMin)       &                                   GGL90mixingLengthMin)
216            rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K)            rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
217           ENDDO           ENDDO
218          ENDDO          ENDDO
219         ENDDO         ENDDO
# Line 218  C- Line 221  C-
221        ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN
222  C-  C-
223         DO k=2,Nr         DO k=2,Nr
224          DO J=jMin,jMax          DO j=jMin,jMax
225           DO I=iMin,iMax           DO i=iMin,iMax
226            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))
227  c         MaxLength=MAX(MaxLength,20. _d 0)  c         MaxLength=MAX(MaxLength,20. _d 0)
228            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
229       &                                   MaxLength)       &                                   MaxLength)
230           ENDDO           ENDDO
231          ENDDO          ENDDO
232         ENDDO         ENDDO
233    
234         DO k=2,Nr         DO k=2,Nr
235          DO J=jMin,jMax          DO j=jMin,jMax
236           DO I=iMin,iMax           DO i=iMin,iMax
237            GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
238       &                                   GGL90mixingLengthMin)       &                                   GGL90mixingLengthMin)
239            rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K)            rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
240           ENDDO           ENDDO
241          ENDDO          ENDDO
242         ENDDO         ENDDO
243    
244        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN
245  C-  C-
246    cgf ensure mixing between first and second level
247    c      DO j=jMin,jMax
248    c        DO i=iMin,iMax
249    c         GGL90mixingLength(i,j,2)=drF(1)
250    c        ENDDO
251    c      ENDDO
252    cgf
253         DO k=2,Nr         DO k=2,Nr
254          DO J=jMin,jMax          DO j=jMin,jMax
255           DO I=iMin,iMax           DO i=iMin,iMax
256            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
257       &        GGL90mixingLength(I,J,K-1)+drF(k-1))       &        GGL90mixingLength(i,j,k-1)+drF(k-1))
258           ENDDO           ENDDO
259          ENDDO          ENDDO
260         ENDDO         ENDDO
261         DO J=jMin,jMax         DO j=jMin,jMax
262          DO I=iMin,iMax          DO i=iMin,iMax
263            GGL90mixingLength(I,J,Nr) = MIN(GGL90mixingLength(I,J,Nr),            GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
264       &       GGL90mixingLengthMin+drF(Nr))       &       GGL90mixingLengthMin+drF(Nr))
265          ENDDO          ENDDO
266         ENDDO         ENDDO
267         DO k=Nr-1,2,-1         DO k=Nr-1,2,-1
268          DO J=jMin,jMax          DO j=jMin,jMax
269           DO I=iMin,iMax           DO i=iMin,iMax
270            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
271       &        GGL90mixingLength(I,J,K+1)+drF(k))       &        GGL90mixingLength(i,j,k+1)+drF(k))
272           ENDDO           ENDDO
273          ENDDO          ENDDO
274         ENDDO         ENDDO
275    
276         DO k=2,Nr         DO k=2,Nr
277          DO J=jMin,jMax          DO j=jMin,jMax
278           DO I=iMin,iMax           DO i=iMin,iMax
279            GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
280       &                                   GGL90mixingLengthMin)       &                                   GGL90mixingLengthMin)
281            rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K)            rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
282           ENDDO           ENDDO
283          ENDDO          ENDDO
284         ENDDO         ENDDO
# Line 276  C- Line 286  C-
286        ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN
287  C-  C-
288         DO k=2,Nr         DO k=2,Nr
289          DO J=jMin,jMax          DO j=jMin,jMax
290           DO I=iMin,iMax           DO i=iMin,iMax
291            mxLength_Dn(I,J,K) = MIN(GGL90mixingLength(I,J,K),            mxLength_Dn(i,j,k) = MIN(GGL90mixingLength(i,j,k),
292       &        mxLength_Dn(I,J,K-1)+drF(k-1))       &        mxLength_Dn(i,j,k-1)+drF(k-1))
293           ENDDO           ENDDO
294          ENDDO          ENDDO
295         ENDDO         ENDDO
296         DO J=jMin,jMax         DO j=jMin,jMax
297          DO I=iMin,iMax          DO i=iMin,iMax
298            GGL90mixingLength(I,J,Nr) = MIN(GGL90mixingLength(I,J,Nr),            GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
299       &       GGL90mixingLengthMin+drF(Nr))       &       GGL90mixingLengthMin+drF(Nr))
300          ENDDO          ENDDO
301         ENDDO         ENDDO
302         DO k=Nr-1,2,-1         DO k=Nr-1,2,-1
303          DO J=jMin,jMax          DO j=jMin,jMax
304           DO I=iMin,iMax           DO i=iMin,iMax
305            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
306       &        GGL90mixingLength(I,J,K+1)+drF(k))       &        GGL90mixingLength(i,j,k+1)+drF(k))
307           ENDDO           ENDDO
308          ENDDO          ENDDO
309         ENDDO         ENDDO
310    
311         DO k=2,Nr         DO k=2,Nr
312          DO J=jMin,jMax          DO j=jMin,jMax
313           DO I=iMin,iMax           DO i=iMin,iMax
314            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
315       &                                  mxLength_Dn(I,J,K))       &                                  mxLength_Dn(i,j,k))
316            tmpmlx = SQRT( GGL90mixingLength(I,J,K)*mxLength_Dn(I,J,K) )            tmpmlx = SQRT( GGL90mixingLength(i,j,k)*mxLength_Dn(i,j,k) )
317            tmpmlx = MAX( tmpmlx, GGL90mixingLengthMin)            tmpmlx = MAX( tmpmlx, GGL90mixingLengthMin)
318            rMixingLength(I,J,K) = 1. _d 0 / tmpmlx            rMixingLength(i,j,k) = 1. _d 0 / tmpmlx
319           ENDDO           ENDDO
320          ENDDO          ENDDO
321         ENDDO         ENDDO
# Line 316  C- Line 326  C-
326    
327  C- Impose minimum mixing length (to avoid division by zero)  C- Impose minimum mixing length (to avoid division by zero)
328  c      DO k=2,Nr  c      DO k=2,Nr
329  c      DO J=jMin,jMax  c      DO j=jMin,jMax
330  c       DO I=iMin,iMax  c       DO i=iMin,iMax
331  c        GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),  c        GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
332  c    &        GGL90mixingLengthMin)  c    &        GGL90mixingLengthMin)
333  c        rMixingLength(I,J,K) = 1. _d 0 /GGL90mixingLength(I,J,K)  c        rMixingLength(i,j,k) = 1. _d 0 /GGL90mixingLength(i,j,k)
334  c       ENDDO  c       ENDDO
335  c      ENDDO  c      ENDDO
336  c     ENDDO  c     ENDDO
337    
338        DO k=2,Nr        DO k=2,Nr
339         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)  
          KappaH = KappaM/TKEPrandtlNumber(I,J,K)  
   
 C     Set a minium (= background) and maximum value  
          KappaM = MAX(KappaM,viscArNr(k))  
          KappaH = MAX(KappaH,diffKrNrT(k))  
          KappaM = MIN(KappaM,GGL90viscMax)  
          KappaH = MIN(KappaH,GGL90diffMax)  
   
 C     Mask land points and storage  
          GGL90viscAr(I,J,K,bi,bj) = KappaM * maskC(I,J,K,bi,bj)  
          GGL90diffKr(I,J,K,bi,bj) = KappaH * maskC(I,J,K,bi,bj)  
          KappaE(I,J,K) = GGL90alpha * GGL90viscAr(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  
340    
341    #ifdef ALLOW_GGL90_HORIZDIFF
342           IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
343  C     horizontal diffusion of TKE (requires an exchange in  C     horizontal diffusion of TKE (requires an exchange in
344  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  
345  C     common factors  C     common factors
346          DO j=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
347           DO i=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
# Line 416  C     ... across y-faces Line 379  C     ... across y-faces
379  C     Compute divergence of fluxes  C     Compute divergence of fluxes
380          DO j=1-Oly,sNy+Oly-1          DO j=1-Oly,sNy+Oly-1
381           DO i=1-Olx,sNx+Olx-1           DO i=1-Olx,sNx+Olx-1
382            gTKE(i,j,k)=gTKE(i,j,k)            gTKE(i,j) =
383       &   -_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)
384       &         *( (dfx(i+1,j)-dfx(i,j))       &         *( (dfx(i+1,j)-dfx(i,j))
385       &           +(dfy(i,j+1)-dfy(i,j))       &           +(dfy(i,j+1)-dfy(i,j))
386       &           )*deltaTggl90       &          )
387           ENDDO           ENDDO
388          ENDDO          ENDDO
389  C       end of k-loop  C      end if GGL90diffTKEh .eq. 0.
390           ENDIF
391    #endif /* ALLOW_GGL90_HORIZDIFF */
392    
393           DO j=jMin,jMax
394            DO i=iMin,iMax
395    C     vertical shear term (dU/dz)^2+(dV/dz)^2
396             tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)
397         &                 -( uVel(i,j,k  ,bi,bj)+uVel(i+1,j,k  ,bi,bj)) )
398         &        *recip_drC(k)
399             tempV= .5 _d 0*( vVel(i,j,km1,bi,bj)+vVel(i,j+1,km1,bi,bj)
400         &                 -( vVel(i,j,k  ,bi,bj)+vVel(i,j+1,k  ,bi,bj)) )
401         &        *recip_drC(k)
402             verticalShear = tempU*tempU + tempV*tempV
403             RiNumber = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)
404    C     compute Prandtl number (always greater than 0)
405             prTemp = 1. _d 0
406             IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
407             TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
408    c         TKEPrandtlNumber(i,j,k) = 1. _d 0
409    
410    C     viscosity and diffusivity
411             KappaM = GGL90ck*GGL90mixingLength(i,j,k)*SQRTTKE(i,j,k)
412             GGL90visctmp(i,j,k) = MAX(KappaM,diffKrNrT(k))
413         &                            * maskC(i,j,k,bi,bj)
414    c        note: storing GGL90visctmp like this, and using it later to compute
415    c              GGL9rdiffKr etc. is robust in case of smoothing (e.g. see OPA)
416             KappaM = MAX(KappaM,viscArNr(k)) * maskC(i,j,k,bi,bj)
417             KappaH = KappaM/TKEPrandtlNumber(i,j,k)
418             KappaE(i,j,k) = GGL90alpha * KappaM * maskC(i,j,k,bi,bj)
419    
420    C     dissipation term
421             TKEdissipation = explDissFac*GGL90ceps
422         &        *SQRTTKE(i,j,k)*rMixingLength(i,j,k)
423         &        *GGL90TKE(i,j,k,bi,bj)
424    C     partial update with sum of explicit contributions
425             GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
426         &        + deltaTggl90*(
427         &        + KappaM*verticalShear
428         &        - KappaH*Nsquare(i,j,k)
429         &        - TKEdissipation
430         &        )
431            ENDDO
432         ENDDO         ENDDO
433  C     end if GGL90diffTKEh .eq. 0.  
434        ENDIF  #ifdef ALLOW_GGL90_HORIZDIFF
435           IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
436    C--    Add horiz. diffusion tendency
437            DO j=jMin,jMax
438             DO i=iMin,iMax
439              GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
440         &                          + gTKE(i,j)*deltaTggl90
441             ENDDO
442            ENDDO
443           ENDIF
444  #endif /* ALLOW_GGL90_HORIZDIFF */  #endif /* ALLOW_GGL90_HORIZDIFF */
445    
446    C--   end of k loop
447          ENDDO
448    
449  C     ============================================  C     ============================================
450  C     Implicit time step to update TKE for k=1,Nr;  C     Implicit time step to update TKE for k=1,Nr;
451  C     TKE(Nr+1)=0 by default  C     TKE(Nr+1)=0 by default
# Line 480  C--   Center diagonal Line 497  C--   Center diagonal
497         DO j=jMin,jMax         DO j=jMin,jMax
498          DO i=iMin,iMax          DO i=iMin,iMax
499            b(i,j,k) = 1. _d 0 - c(i,j,k) - a(i,j,k)            b(i,j,k) = 1. _d 0 - c(i,j,k) - a(i,j,k)
500       &        + ab15*deltaTggl90*GGL90ceps*SQRTTKE(I,J,K)       &        + implDissFac*deltaTggl90*GGL90ceps*SQRTTKE(i,j,k)
501       &        * rMixingLength(I,J,K)       &        * rMixingLength(i,j,k)
502       &        * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)       &        * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)
503           ENDDO           ENDDO
504         ENDDO         ENDDO
# Line 490  C     end set up matrix Line 507  C     end set up matrix
507    
508  C     Apply boundary condition  C     Apply boundary condition
509        kp1 = MIN(Nr,kSurf+1)        kp1 = MIN(Nr,kSurf+1)
510        DO J=jMin,jMax        DO j=jMin,jMax
511         DO I=iMin,iMax         DO i=iMin,iMax
512  C     estimate friction velocity uStar from surface forcing  C     estimate friction velocity uStar from surface forcing
513          uStarSquare = SQRT(          uStarSquare = SQRT(
514       &    ( .5 _d 0*( surfaceForcingU(I,  J,  bi,bj)       &    ( .5 _d 0*( surfaceForcingU(i,  j,  bi,bj)
515       &              + surfaceForcingU(I+1,J,  bi,bj) ) )**2       &              + surfaceForcingU(i+1,j,  bi,bj) ) )**2
516       &  + ( .5 _d 0*( surfaceForcingV(I,  J,  bi,bj)       &  + ( .5 _d 0*( surfaceForcingV(i,  j,  bi,bj)
517       &              + surfaceForcingV(I,  J+1,bi,bj) ) )**2       &              + surfaceForcingV(i,  j+1,bi,bj) ) )**2
518       &                     )       &                     )
519  C     Dirichlet surface boundary condition for TKE  C     Dirichlet surface boundary condition for TKE
520          gTKE(I,J,kSurf) = MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)          GGL90TKE(i,j,kSurf,bi,bj) = maskC(i,j,kSurf,bi,bj)
521       &                     *maskC(I,J,kSurf,bi,bj)       &           *MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)
522          gTKE(i,j,kp1) = gTKE(i,j,kp1)          GGL90TKE(i,j,kp1,bi,bj) = GGL90TKE(i,j,kp1,bi,bj)
523       &                - a(i,j,kp1)*gTKE(i,j,kSurf)       &               - a(i,j,kp1)*GGL90TKE(i,j,kSurf,bi,bj)
524          a(i,j,kp1) = 0. _d 0          a(i,j,kp1) = 0. _d 0
525  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom
526          kBottom   = MAX(kLowC(I,J,bi,bj),1)          kBottom   = MAX(kLowC(i,j,bi,bj),1)
527          gTKE(I,J,kBottom) = gTKE(I,J,kBottom)          GGL90TKE(i,j,kBottom,bi,bj) = GGL90TKE(i,j,kBottom,bi,bj)
528       &                    - GGL90TKEbottom*c(I,J,kBottom)       &                              - GGL90TKEbottom*c(i,j,kBottom)
529          c(I,J,kBottom) = 0. _d 0          c(i,j,kBottom) = 0. _d 0
530         ENDDO         ENDDO
531        ENDDO        ENDDO
532    
533  C     solve tri-diagonal system, and store solution on gTKE (previously rhs)  C     solve tri-diagonal system
534        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
535       I                        a, b, c,       I                        a, b, c,
536       U                        gTKE,       U                        GGL90TKE,
537       O                        errCode,       O                        errCode,
538       I                        bi, bj, myThid )       I                        bi, bj, myThid )
539  c     CALL GGL90_SOLVE( bi, bj, iMin, iMax, jMin, jMax,  
540  c    I     a, b, c,        DO k=1,Nr
541  c    U     gTKE,         DO j=jMin,jMax
542  c    I     myThid )          DO i=iMin,iMax
   
 C     now update TKE  
       DO K=1,Nr  
        DO J=jMin,jMax  
         DO I=iMin,iMax  
543  C     impose minimum TKE to avoid numerical undershoots below zero  C     impose minimum TKE to avoid numerical undershoots below zero
544           GGL90TKE(I,J,K,bi,bj) = MAX( gTKE(I,J,K), GGL90TKEmin )           GGL90TKE(i,j,k,bi,bj) = maskC(i,j,k,bi,bj)
545       &        * maskC(I,J,K,bi,bj)       &                  *MAX( GGL90TKE(i,j,k,bi,bj), GGL90TKEmin )
546          ENDDO          ENDDO
547         ENDDO         ENDDO
548        ENDDO        ENDDO
# Line 538  C     impose minimum TKE to avoid numeri Line 550  C     impose minimum TKE to avoid numeri
550  C     end of time step  C     end of time step
551  C     ===============================  C     ===============================
552    
553          DO k=2,Nr
554           DO j=1,sNy
555            DO i=1,sNx
556  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
557        DO K=1,Nr           tmpVisc=
        DO J=jMin,jMax  
         DO I=iMin,iMax  
          tmpdiffKrS=  
558       &  (       &  (
559       &   p4 *  GGL90viscAr(i  ,j  ,k,bi,bj) * mskCor(i  ,j  ,bi,bj)       &   p4 *  GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
560       &  +p8 *( GGL90viscAr(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)       &  +p8 *( GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj)
561       &       + GGL90viscAr(i  ,j-1,k,bi,bj) * mskCor(i  ,j-1,bi,bj)       &       + GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj)
562       &       + GGL90viscAr(i+1,j  ,k,bi,bj) * mskCor(i+1,j  ,bi,bj)       &       + GGL90visctmp(i+1,j  ,k) * mskCor(i+1,j  ,bi,bj)
563       &       + GGL90viscAr(i  ,j+1,k,bi,bj) * mskCor(i  ,j+1,bi,bj))       &       + GGL90visctmp(i  ,j+1,k) * mskCor(i  ,j+1,bi,bj))
564       &  +p16*( GGL90viscAr(i+1,j+1,k,bi,bj) * mskCor(i+1,j+1,bi,bj)       &  +p16*( GGL90visctmp(i+1,j+1,k) * mskCor(i+1,j+1,bi,bj)
565       &       + GGL90viscAr(i+1,j-1,k,bi,bj) * mskCor(i+1,j-1,bi,bj)       &       + GGL90visctmp(i+1,j-1,k) * mskCor(i+1,j-1,bi,bj)
566       &       + GGL90viscAr(i-1,j+1,k,bi,bj) * mskCor(i-1,j+1,bi,bj)       &       + GGL90visctmp(i-1,j+1,k) * mskCor(i-1,j+1,bi,bj)
567       &       + GGL90viscAr(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj))       &       + GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj))
568       &  )       &  )
569       & /(p4       & /(p4
570       &  +p8 *(       maskC(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)       &  +p8 *(       maskC(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)
# Line 564  C     =============================== Line 576  C     ===============================
576       &       +       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)
577       &       +       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))
578       &  )*maskC(i,j,k,bi,bj)*mskCor(i,j,bi,bj)       &  )*maskC(i,j,k,bi,bj)*mskCor(i,j,bi,bj)
579       &   /TKEPrandtlNumber(i,j,k)  #else
580           GGL90diffKrS(I,J,K,bi,bj)= MAX( tmpdiffKrS , diffKrNrT(k) )           tmpVisc = GGL90visctmp(i,j,k)
581    #endif
582             tmpVisc = MIN(tmpVisc/TKEPrandtlNumber(i,j,k),GGL90diffMax)
583             GGL90diffKr(i,j,k,bi,bj)= MAX( tmpVisc , diffKrNrT(k) )
584            ENDDO
585           ENDDO
586          ENDDO
587    
588          DO k=2,Nr
589           DO j=1,sNy
590            DO i=1,sNx+1
591    #ifdef ALLOW_GGL90_SMOOTH
592            tmpVisc =
593         & (
594         &   p4 *(GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
595         &       +GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj))
596         &  +p8 *(GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj)
597         &       +GGL90visctmp(i-1,j+1,k) * mskCor(i-1,j+1,bi,bj)
598         &       +GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj)
599         &       +GGL90visctmp(i  ,j+1,k) * mskCor(i  ,j+1,bi,bj))
600         &  )
601         & /(p4 * 2. _d 0
602         &  +p8 *(      maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj)
603         &       +      maskC(i-1,j+1,k,bi,bj) * mskCor(i-1,j+1,bi,bj)
604         &       +      maskC(i  ,j-1,k,bi,bj) * mskCor(i  ,j-1,bi,bj)
605         &       +      maskC(i  ,j+1,k,bi,bj) * mskCor(i  ,j+1,bi,bj))
606         &  )
607         &  *maskC(i  ,j,k,bi,bj)*mskCor(i  ,j,bi,bj)
608         &  *maskC(i-1,j,k,bi,bj)*mskCor(i-1,j,bi,bj)
609    #else
610            tmpVisc = _maskW(i,j,k,bi,bj) *
611         &                   (.5 _d 0*(GGL90visctmp(i,j,k)
612         &                            +GGL90visctmp(i-1,j,k))
613         &                   )
614    #endif
615            tmpVisc = MIN( tmpVisc , GGL90viscMax )
616            GGL90viscArU(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
617          ENDDO          ENDDO
618         ENDDO         ENDDO
619        ENDDO        ENDDO
620    
621          DO k=2,Nr
622           DO j=1,sNy+1
623            DO i=1,sNx
624    #ifdef ALLOW_GGL90_SMOOTH
625            tmpVisc =
626         & (
627         &   p4 *(GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
628         &       +GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj))
629         &  +p8 *(GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj)
630         &       +GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj)
631         &       +GGL90visctmp(i+1,j  ,k) * mskCor(i+1,j  ,bi,bj)
632         &       +GGL90visctmp(i+1,j-1,k) * mskCor(i+1,j-1,bi,bj))
633         &  )
634         & /(p4 * 2. _d 0
635         &  +p8 *(      maskC(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)
636         &       +      maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj)
637         &       +      maskC(i+1,j  ,k,bi,bj) * mskCor(i+1,j  ,bi,bj)
638         &       +      maskC(i+1,j-1,k,bi,bj) * mskCor(i+1,j-1,bi,bj))
639         &  )
640         &   *maskC(i,j  ,k,bi,bj)*mskCor(i,j  ,bi,bj)
641         &   *maskC(i,j-1,k,bi,bj)*mskCor(i,j-1,bi,bj)
642    #else
643            tmpVisc = _maskS(i,j,k,bi,bj) *
644         &                   (.5 _d 0*(GGL90visctmp(i,j,k)
645         &                            +GGL90visctmp(i,j-1,k))
646         &                   )
647    
648  #endif  #endif
649            tmpVisc = MIN( tmpVisc , GGL90viscMax )
650            GGL90viscArV(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
651            ENDDO
652           ENDDO
653          ENDDO
654    
655  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
656        IF ( useDiagnostics ) THEN        IF ( useDiagnostics ) THEN
657           CALL DIAGNOSTICS_FILL( GGL90TKE   ,'GGL90TKE',           CALL DIAGNOSTICS_FILL( GGL90TKE   ,'GGL90TKE',
658       &                          0,Nr, 1, bi, bj, myThid )       &                          0,Nr, 1, bi, bj, myThid )
659           CALL DIAGNOSTICS_FILL( GGL90viscAr,'GGL90Ar ',           CALL DIAGNOSTICS_FILL( GGL90viscArU,'GGL90ArU',
660         &                          0,Nr, 1, bi, bj, myThid )
661             CALL DIAGNOSTICS_FILL( GGL90viscArV,'GGL90ArV',
662       &                          0,Nr, 1, bi, bj, myThid )       &                          0,Nr, 1, bi, bj, myThid )
663           CALL DIAGNOSTICS_FILL( GGL90diffKr,'GGL90Kr ',           CALL DIAGNOSTICS_FILL( GGL90diffKr,'GGL90Kr ',
664       &                          0,Nr, 1, bi, bj, myThid )       &                          0,Nr, 1, bi, bj, myThid )
Legend:
Removed from v.1.15  
changed lines
  Added in v.1.19
  ViewVC Help
Powered by ViewVC 1.1.22