/[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.3 by jmc, Sun Oct 17 23:11:59 2004 UTC revision 1.10 by dfer, Tue Oct 7 19:35:10 2008 UTC
# Line 60  C     K, Kp1, km1, kSurf, kBottom  - ver Line 60  C     K, Kp1, km1, kSurf, kBottom  - ver
60  C     ab15, ab05       - weights for implicit timestepping  C     ab15, ab05       - weights for implicit timestepping
61  C     uStarSquare      - square of friction velocity  C     uStarSquare      - square of friction velocity
62  C     verticalShear    - (squared) vertical shear of horizontal velocity  C     verticalShear    - (squared) vertical shear of horizontal velocity
63  C     Nsquare          - squared buoyancy freqency  C     Nsquare          - squared buoyancy freqency
64  C     RiNumber         - local Richardson number  C     RiNumber         - local Richardson number
65  C     KappaM           - (local) viscosity parameter (eq.10)  C     KappaM           - (local) viscosity parameter (eq.10)
66  C     KappaH           - (local) diffusivity parameter for temperature (eq.11)  C     KappaH           - (local) diffusivity parameter for temperature (eq.11)
67  C     KappaE           - (local) diffusivity parameter for TKE (eq.15)  C     KappaE           - (local) diffusivity parameter for TKE (eq.15)
 C     buoyFreq         - buoyancy freqency  
68  C     TKEdissipation   - dissipation of TKE  C     TKEdissipation   - dissipation of TKE
69  C     GGL90mixingLength- mixing length of scheme following Banke+Delecuse  C     GGL90mixingLength- mixing length of scheme following Banke+Delecuse
70    C         rMixingLength- inverse of mixing length
71  C     totalDepth       - thickness of water column (inverse of recip_Rcol)  C     totalDepth       - thickness of water column (inverse of recip_Rcol)
72  C     TKEPrandtlNumber - here, an empirical function of the Richardson number  C     TKEPrandtlNumber - here, an empirical function of the Richardson number
73  C     rhoK, rhoKm1     - density at layer K and Km1 (relative to K)  C     rhoK, rhoKm1     - density at layer K and Km1 (relative to K)
# Line 86  C     gTKE             - right hand side Line 86  C     gTKE             - right hand side
86        _RL     tempU, tempV, prTemp        _RL     tempU, tempV, prTemp
87        _RL     TKEPrandtlNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     TKEPrandtlNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
88        _RL     GGL90mixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     GGL90mixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
89          _RL         rMixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
90        _RL     KappaE           (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     KappaE           (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
91        _RL     rhoK             (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     rhoK             (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92        _RL     rhoKm1           (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     rhoKm1           (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 110  CEOP Line 111  CEOP
111        jMax = sNy+OLy-1        jMax = sNy+OLy-1
112    
113  C     set separate time step (should be deltaTtracer)  C     set separate time step (should be deltaTtracer)
114        deltaTggl90 = deltaTtracer        deltaTggl90 = dTtracerLev(1)
115  C      C
116        kSurf = 1        kSurf = 1
117  C     implicit timestepping weights for dissipation  C     implicit timestepping weights for dissipation
118        ab15 =  1.5 _d 0        ab15 =  1.5 _d 0
# Line 126  C     Initialize local fields Line 127  C     Initialize local fields
127           gTKE(I,J,K)              = 0. _d 0           gTKE(I,J,K)              = 0. _d 0
128           KappaE(I,J,K)            = 0. _d 0           KappaE(I,J,K)            = 0. _d 0
129           TKEPrandtlNumber(I,J,K)  = 0. _d 0           TKEPrandtlNumber(I,J,K)  = 0. _d 0
130           GGL90mixingLength(I,J,K) = 0. _d 0           GGL90mixingLength(I,J,K) = GGL90mixingLengthMin
131                 rMixingLength(I,J,K) = 0. _d 0
132          ENDDO          ENDDO
133         ENDDO             ENDDO
134        ENDDO        ENDDO
135        DO J=1-Oly,sNy+Oly        DO J=1-Oly,sNy+Oly
136         DO I=1-Olx,sNx+Olx         DO I=1-Olx,sNx+Olx
137          rhoK    (I,J) = 0. _d 0          rhoK    (I,J) = 0. _d 0
138          rhoKm1  (I,J) = 0. _d 0          rhoKm1  (I,J) = 0. _d 0
139          totalDepth(I,J)   = 0. _d 0          totalDepth(I,J)   = 0. _d 0
140          IF ( recip_Rcol(I,J,bi,bj) .NE. 0. )          IF ( recip_Rcol(I,J,bi,bj) .NE. 0. _d 0 )
141       &       totalDepth(I,J) = 1./recip_Rcol(I,J,bi,bj)       &       totalDepth(I,J) = 1./recip_Rcol(I,J,bi,bj)
142         ENDDO         ENDDO
143        ENDDO        ENDDO
# Line 144  C     start k-loop Line 146  C     start k-loop
146        DO K = 2, Nr        DO K = 2, Nr
147         Km1 = K-1         Km1 = K-1
148         Kp1 = MIN(Nr,K+1)         Kp1 = MIN(Nr,K+1)
149         CALL FIND_RHO(         CALL FIND_RHO_2D(
150       I      bi, bj, iMin, iMax, jMin, jMax, Km1, K,       I      iMin, iMax, jMin, jMax, K,
151       I      theta, salt,       I      theta(1-OLx,1-OLy,Km1,bi,bj), salt(1-OLx,1-OLy,Km1,bi,bj),
152       O      rhoKm1,       O      rhoKm1,
153       I      myThid )       I      Km1, bi, bj, myThid )
154         CALL FIND_RHO(  
155       I      bi, bj, iMin, iMax, jMin, jMax, K, K,         CALL FIND_RHO_2D(
156       I      theta, salt,       I      iMin, iMax, jMin, jMax, K,
157         I      theta(1-OLx,1-OLy,K,bi,bj), salt(1-OLx,1-OLy,K,bi,bj),
158       O      rhoK,       O      rhoK,
159       I      myThid )       I      K, bi, bj, myThid )
160         DO J=jMin,jMax         DO J=jMin,jMax
161          DO I=iMin,iMax          DO I=iMin,iMax
162           SQRTTKE=SQRT( GGL90TKE(I,J,K,bi,bj) )           SQRTTKE=SQRT( GGL90TKE(I,J,K,bi,bj) )
# Line 163  C Line 166  C
166           Nsquare = - gravity*recip_rhoConst*recip_drC(K)           Nsquare = - gravity*recip_rhoConst*recip_drC(K)
167       &        * ( rhoKm1(I,J) - rhoK(I,J) )*maskC(I,J,K,bi,bj)       &        * ( rhoKm1(I,J) - rhoK(I,J) )*maskC(I,J,K,bi,bj)
168  C     vertical shear term (dU/dz)^2+(dV/dz)^2  C     vertical shear term (dU/dz)^2+(dV/dz)^2
169           tempu= .5*(  uVel(I,J,Km1,bi,bj)+uVel(I+1,J,Km1,bi,bj)           tempu= .5 _d 0*( uVel(I,J,Km1,bi,bj)+uVel(I+1,J,Km1,bi,bj)
170       &             - (uVel(I,J,K  ,bi,bj)+uVel(I+1,J,K  ,bi,bj)) )       &                 -( uVel(I,J,K  ,bi,bj)+uVel(I+1,J,K  ,bi,bj)) )
171       &        *recip_drC(K)       &        *recip_drC(K)
172           tempv= .5*(  vVel(I,J,Km1,bi,bj)+vVel(I,J+1,Km1,bi,bj)           tempv= .5 _d 0*( vVel(I,J,Km1,bi,bj)+vVel(I,J+1,Km1,bi,bj)
173       &             - (vVel(I,J,K  ,bi,bj)+vVel(I,J+1,K  ,bi,bj)) )       &                 -( vVel(I,J,K  ,bi,bj)+vVel(I,J+1,K  ,bi,bj)) )
174       &        *recip_drC(K)       &        *recip_drC(K)
175           verticalShear = tempU*tempU + tempV*tempV           verticalShear = tempU*tempU + tempV*tempV
176           RiNumber   = MAX(Nsquare,0. _d 0)/(verticalShear+GGL90eps)           RiNumber   = MAX(Nsquare,0. _d 0)/(verticalShear+GGL90eps)
177  C     compute Prandtl number (always greater than 0)  C     compute Prandtl number (always greater than 0)
178           prTemp = 1. _d 0           prTemp = 1. _d 0
179           IF ( RiNumber .GE. 0.2 ) prTemp = 5.0 * RiNumber           IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
180           TKEPrandtlNumber(I,J,K) = MIN(10.,prTemp)           TKEPrandtlNumber(I,J,K) = MIN(10. _d 0,prTemp)
181  C     mixing length  C     mixing length
182           GGL90mixingLength(I,J,K) =           GGL90mixingLength(I,J,K) = SQRTTWO *
183       &        SQRTTKE/SQRT( MAX(Nsquare,GGL90eps) )       &        SQRTTKE/SQRT( MAX(Nsquare,GGL90eps) )
184  C     impose upper bound for mixing length (total depth)  C     impose upper bound for mixing length (total depth)
185           GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),           GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),
186       &        totalDepth(I,J))       &        totalDepth(I,J))
187  C     impose minimum mixing length (to avoid division by zero)  C     impose minimum mixing length (to avoid division by zero)
188           GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),           GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),
189       &        GGL90mixingLengthMin)       &        GGL90mixingLengthMin)
190             rMixingLength(I,J,K) = 1. _d 0 /GGL90mixingLength(I,J,K)
191  C     viscosity of last timestep  C     viscosity of last timestep
192           KappaM = GGL90ck*GGL90mixingLength(I,J,K)*SQRTTKE           KappaM = GGL90ck*GGL90mixingLength(I,J,K)*SQRTTKE
193             KappaH = KappaM/TKEPrandtlNumber(I,J,K)
194    
195    C     Set a minium (= background) and maximum value
196             KappaM = MAX(KappaM,viscAr)
197             KappaH = MAX(KappaH,diffKrNrT(k))
198             KappaM = MIN(KappaM,GGL90viscMax)
199             KappaH = MIN(KappaH,GGL90diffMax)
200    
201    C     Mask land points and save
202           KappaE(I,J,K) = KappaM*GGL90alpha           KappaE(I,J,K) = KappaM*GGL90alpha
203             GGL90viscAr(I,J,K,bi,bj) = KappaM * maskC(I,J,K,bi,bj)
204             GGL90diffKr(I,J,K,bi,bj) = KappaH * maskC(I,J,K,bi,bj)
205    
206  C     dissipation term  C     dissipation term
207           TKEdissipation = ab05*GGL90ceps           TKEdissipation = ab05*GGL90ceps
208       &        *SQRTTKE/GGL90mixingLength(I,J,K)       &        *SQRTTKE*rMixingLength(I,J,K)
209       &        *GGL90TKE(I,J,K,bi,bj)             &        *GGL90TKE(I,J,K,bi,bj)
210  C     sum up contributions to form the right hand side  C     sum up contributions to form the right hand side
211           gTKE(I,J,K) = GGL90TKE(I,J,K,bi,bj)           gTKE(I,J,K) = GGL90TKE(I,J,K,bi,bj)
212       &        + deltaTggl90*(       &        + deltaTggl90*(
213       &        + KappaM*verticalShear       &        + KappaM*verticalShear
214       &        - KappaM*Nsquare/TKEPrandtlNumber(I,J,K)       &        - KappaH*Nsquare
215       &        - TKEdissipation  c    &        - KappaM*Nsquare/TKEPrandtlNumber(I,J,K)
216         &        - TKEdissipation
217       &        )       &        )
218          ENDDO            ENDDO
219         ENDDO         ENDDO
220        ENDDO        ENDDO
221  C     horizontal diffusion of TKE (requires an exchange in  C     horizontal diffusion of TKE (requires an exchange in
222  C       do_fields_blocking_exchanges)  C      do_fields_blocking_exchanges)
223  #ifdef ALLOW_GGL90_HORIZDIFF  #ifdef ALLOW_GGL90_HORIZDIFF
224        IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN        IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
225         DO K=2,Nr         DO K=2,Nr
# Line 214  C     common factors Line 231  C     common factors
231            yA(i,j) = _dxG(i,j,bi,bj)            yA(i,j) = _dxG(i,j,bi,bj)
232       &         *drF(k)*_hFacS(i,j,k,bi,bj)       &         *drF(k)*_hFacS(i,j,k,bi,bj)
233           ENDDO           ENDDO
234          ENDDO            ENDDO
235  C     Compute diffusive fluxes  C     Compute diffusive fluxes
236  C     ... across x-faces  C     ... across x-faces
237          DO j=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
238           dfx(1-Olx,j)=0.           dfx(1-Olx,j)=0. _d 0
239           DO i=1-Olx+1,sNx+Olx           DO i=1-Olx+1,sNx+Olx
240            dfx(i,j) = -GGL90diffTKEh*xA(i,j)            dfx(i,j) = -GGL90diffTKEh*xA(i,j)
241       &      *_recip_dxC(i,j,bi,bj)       &      *_recip_dxC(i,j,bi,bj)
# Line 228  C     ... across x-faces Line 245  C     ... across x-faces
245          ENDDO          ENDDO
246  C     ... across y-faces  C     ... across y-faces
247          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
248           dfy(i,1-Oly)=0.           dfy(i,1-Oly)=0. _d 0
249          ENDDO          ENDDO
250          DO j=1-Oly+1,sNy+Oly          DO j=1-Oly+1,sNy+Oly
251           DO i=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
# Line 239  C     ... across y-faces Line 256  C     ... across y-faces
256       &      *CosFacV(j,bi,bj)       &      *CosFacV(j,bi,bj)
257  #endif /* ISOTROPIC_COS_SCALING */  #endif /* ISOTROPIC_COS_SCALING */
258           ENDDO           ENDDO
259          ENDDO            ENDDO
260  C     Compute divergence of fluxes  C     Compute divergence of fluxes
261          DO j=1-Oly,sNy+Oly-1          DO j=1-Oly,sNy+Oly-1
262           DO i=1-Olx,sNx+Olx-1           DO i=1-Olx,sNx+Olx-1
263            gTKE(i,j,k)=gTKE(i,j,k)            gTKE(i,j,k)=gTKE(i,j,k)
264       &   -_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)
265       &         *( (dfx(i+1,j)-dfx(i,j))       &         *( (dfx(i+1,j)-dfx(i,j))
266       &           +(dfy(i,j+1)-dfy(i,j))       &           +(dfy(i,j+1)-dfy(i,j))
267       &           )       &           )
268           ENDDO             ENDDO
269          ENDDO          ENDDO
270  C       end of k-loop  C       end of k-loop
271         ENDDO         ENDDO
272  C     end if GGL90diffTKEh .eq. 0.  C     end if GGL90diffTKEh .eq. 0.
273        ENDIF        ENDIF
# Line 264  C     set up matrix Line 281  C     set up matrix
281  C--   Lower diagonal  C--   Lower diagonal
282        DO j=jMin,jMax        DO j=jMin,jMax
283         DO i=iMin,iMax         DO i=iMin,iMax
284           a(i,j,1) = 0. _d 0           a(i,j,1) = 0. _d 0
285         ENDDO         ENDDO
286        ENDDO        ENDDO
287        DO k=2,Nr        DO k=2,Nr
# Line 273  C--   Lower diagonal Line 290  C--   Lower diagonal
290          DO i=iMin,iMax          DO i=iMin,iMax
291           a(i,j,k) = -deltaTggl90           a(i,j,k) = -deltaTggl90
292       &        *recip_drF(km1)*recip_hFacI(i,j,k,bi,bj)       &        *recip_drF(km1)*recip_hFacI(i,j,k,bi,bj)
293       &        *.5*(KappaE(i,j, k )+KappaE(i,j,km1))       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))
294       &        *recip_drC(k)       &        *recip_drC(k)
295            IF (recip_hFacI(i,j,km1,bi,bj).EQ.0.) a(i,j,k)=0.            IF (recip_hFacI(i,j,km1,bi,bj).EQ.0. _d 0) a(i,j,k)=0. _d 0
296          ENDDO          ENDDO
297         ENDDO         ENDDO
298        ENDDO        ENDDO
# Line 293  CML      DO k=1,Nr-1 Line 310  CML      DO k=1,Nr-1
310          DO i=iMin,iMax          DO i=iMin,iMax
311            c(i,j,k) = -deltaTggl90            c(i,j,k) = -deltaTggl90
312       &        *recip_drF( k )*recip_hFacI(i,j,k,bi,bj)       &        *recip_drF( k )*recip_hFacI(i,j,k,bi,bj)
313       &               *.5*(KappaE(i,j,k)+KappaE(i,j,kp1))       &               *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))
314       &        *recip_drC(k)       &        *recip_drC(k)
315            IF (recip_hFacI(i,j,kp1,bi,bj).EQ.0.) c(i,j,k)=0.            IF (recip_hFacI(i,j,kp1,bi,bj).EQ.0. _d 0) c(i,j,k)=0. _d 0
316          ENDDO          ENDDO
317         ENDDO         ENDDO
318        ENDDO        ENDDO
# Line 305  C--   Center diagonal Line 322  C--   Center diagonal
322          DO i=iMin,iMax          DO i=iMin,iMax
323            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)
324       &        + ab15*deltaTggl90*GGL90ceps*SQRT(GGL90TKE(I,J,K,bi,bj))       &        + ab15*deltaTggl90*GGL90ceps*SQRT(GGL90TKE(I,J,K,bi,bj))
325       &        /GGL90mixingLength(I,J,K)       &        *rMixingLength(I,J,K)
326           ENDDO           ENDDO
327         ENDDO         ENDDO
328        ENDDO        ENDDO
# Line 313  C     end set up matrix Line 330  C     end set up matrix
330    
331  C  C
332  C     Apply boundary condition  C     Apply boundary condition
333  C      C
334        DO J=jMin,jMax        DO J=jMin,jMax
335         DO I=iMin,iMax         DO I=iMin,iMax
336  C     estimate friction velocity uStar from surface forcing  C     estimate friction velocity uStar from surface forcing
337          uStarSquare = SQRT(          uStarSquare = SQRT(
338       &         ( .5*( surfaceForcingU(I,  J,  bi,bj)       &    ( .5 _d 0*( surfaceForcingU(I,  J,  bi,bj)
339       &              + surfaceForcingU(I+1,J,  bi,bj) ) )**2       &              + surfaceForcingU(I+1,J,  bi,bj) ) )**2
340       &       + ( .5*( surfaceForcingV(I,  J,  bi,bj)       &  + ( .5 _d 0*( surfaceForcingV(I,  J,  bi,bj)
341       &              + surfaceForcingV(I,  J+1,bi,bj) ) )**2       &              + surfaceForcingV(I,  J+1,bi,bj) ) )**2
342       &                     )       &                     )
343  C     Dirichlet surface boundary condition for TKE  C     Dirichlet surface boundary condition for TKE
344          gTKE(I,J,kSurf) = MAX(GGL90TKEmin,GGL90m2*uStarSquare)          gTKE(I,J,kSurf) = MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)
345       &                     *maskC(I,J,kSurf,bi,bj)       &                     *maskC(I,J,kSurf,bi,bj)
346  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom
347          kBottom   = MIN(MAX(kLowC(I,J,bi,bj),1),Nr)          kBottom   = MIN(MAX(kLowC(I,J,bi,bj),1),Nr)
348          gTKE(I,J,kBottom) = gTKE(I,J,kBottom)          gTKE(I,J,kBottom) = gTKE(I,J,kBottom)
349       &       - GGL90TKEbottom*c(I,J,kBottom)       &       - GGL90TKEbottom*c(I,J,kBottom)
350         ENDDO         ENDDO
351        ENDDO            ENDDO
352  C  C
353  C     solve tri-diagonal system, and store solution on gTKE (previously rhs)  C     solve tri-diagonal system, and store solution on gTKE (previously rhs)
354  C  C
# Line 341  C Line 358  C
358       I     myThid )       I     myThid )
359  C  C
360  C     now update TKE  C     now update TKE
361  C      C
362        DO K=1,Nr        DO K=1,Nr
363         DO J=jMin,jMax         DO J=jMin,jMax
364          DO I=iMin,iMax          DO I=iMin,iMax
365  C     impose minimum TKE to avoid numerical undershoots below zero  C     impose minimum TKE to avoid numerical undershoots below zero
366           GGL90TKE(I,J,K,bi,bj) = MAX( gTKE(I,J,K), GGL90TKEmin )           GGL90TKE(I,J,K,bi,bj) = MAX( gTKE(I,J,K), GGL90TKEmin )
367       &        * maskC(I,J,K,bi,bj)       &        * maskC(I,J,K,bi,bj)
368          ENDDO          ENDDO
369         ENDDO         ENDDO
370        ENDDO            ENDDO
371  C  C
372  C     end of time step  C     end of time step
373  C     ===============================  C     ===============================
374  C     compute viscosity coefficients  C     compute viscosity coefficients
375  C      C
376        DO K=2,Nr  c      DO K=2,Nr
377         DO J=jMin,jMax  c       DO J=jMin,jMax
378          DO I=iMin,iMax  c        DO I=iMin,iMax
379  C     Eq. (11), (18) and (21)  C     Eq. (11), (18) and (21)
380           KappaM = GGL90ck*GGL90mixingLength(I,J,K)*  c         KappaM = GGL90ck*GGL90mixingLength(I,J,K)*
381       &                  SQRT( GGL90TKE(I,J,K,bi,bj) )  c     &                  SQRT( GGL90TKE(I,J,K,bi,bj) )
382           KappaH = KappaM/TKEPrandtlNumber(I,J,K)  c         KappaH = KappaM/TKEPrandtlNumber(I,J,K)
383  C     Set a minium (= background) value  C     Set a minium (= background) value
384           KappaM = MAX(KappaM,viscAr)  c         KappaM = MAX(KappaM,viscAr)
385           KappaH = MAX(KappaH,diffKrNrT(k))  c         KappaH = MAX(KappaH,diffKrNrT(k))
386  C     Set a maximum and mask land point  C     Set a maximum and mask land point
387           GGL90viscAr(I,J,K,bi,bj) = MIN(KappaM,GGL90viscMax)  c        GGL90viscAr(I,J,K,bi,bj) = MIN(KappaM,GGL90viscMax)
388       &        * maskC(I,J,K,bi,bj)  c    &        * maskC(I,J,K,bi,bj)
389           GGL90diffKr(I,J,K,bi,bj) = MIN(KappaH,GGL90diffMax)  c        GGL90diffKr(I,J,K,bi,bj) = MIN(KappaH,GGL90diffMax)
390       &        * maskC(I,J,K,bi,bj)  c    &        * maskC(I,J,K,bi,bj)
391          ENDDO  c        ENDDO
392         ENDDO  c       ENDDO
393  C     end third k-loop  C     end third k-loop
394        ENDDO      c      ENDDO
395    
396    #ifdef ALLOW_DIAGNOSTICS
397          IF ( useDiagnostics ) THEN
398             CALL DIAGNOSTICS_FILL( GGL90TKE   ,'GGL90TKE',
399         &                          0,Nr, 1, bi, bj, myThid )
400             CALL DIAGNOSTICS_FILL( GGL90viscAr,'GGL90Ar ',
401         &                          0,Nr, 1, bi, bj, myThid )
402             CALL DIAGNOSTICS_FILL( GGL90diffKr,'GGL90Kr ',
403         &                          0,Nr, 1, bi, bj, myThid )
404             CALL DIAGNOSTICS_FILL( TKEPrandtlNumber ,'GGL90Prl',
405         &                          0,Nr, 2, bi, bj, myThid )
406             CALL DIAGNOSTICS_FILL( GGL90mixingLength,'GGL90Lmx',
407         &                          0,Nr, 2, bi, bj, myThid )
408          ENDIF
409    #endif
410    
411  #endif /* ALLOW_GGL90 */  #endif /* ALLOW_GGL90 */
412    
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.10
  ViewVC Help
Powered by ViewVC 1.1.22