/[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.4 by jmc, Sat Dec 4 00:16:49 2004 UTC revision 1.23 by gforget, Wed Aug 8 22:22:42 2012 UTC
# Line 7  CBOP Line 7  CBOP
7  C !ROUTINE: GGL90_CALC  C !ROUTINE: GGL90_CALC
8    
9  C !INTERFACE: ======================================================  C !INTERFACE: ======================================================
10        subroutine GGL90_CALC(        SUBROUTINE GGL90_CALC(
11       I     bi, bj, myTime, myThid )       I                 bi, bj, sigmaR, myTime, myIter, myThid )
12    
13    
14  C !DESCRIPTION: \bv  C !DESCRIPTION: \bv
15  C     /==========================================================\  C     *==========================================================*
16  C     | SUBROUTINE GGL90_CALC                                    |  C     | SUBROUTINE GGL90_CALC                                    |
17  C     | o Compute all GGL90 fields defined in GGL90.h            |  C     | o Compute all GGL90 fields defined in GGL90.h            |
18  C     |==========================================================|  C     *==========================================================*
19  C     | Equation numbers refer to                                |  C     | Equation numbers refer to                                |
20  C     | Gaspar et al. (1990), JGR 95 (C9), pp 16,179             |  C     | Gaspar et al. (1990), JGR 95 (C9), pp 16,179             |
21  C     | Some parts of the implementation follow Blanke and       |  C     | Some parts of the implementation follow Blanke and       |
22  C     | Delecuse (1993), JPO, and OPA code, in particular the    |  C     | Delecuse (1993), JPO, and OPA code, in particular the    |
23  C     | computation of the                                       |  C     | computation of the                                       |
24  C     | mixing length = max(min(lk,depth),lkmin)                 |  C     | mixing length = max(min(lk,depth),lkmin)                 |
25  C     \==========================================================/  C     *==========================================================*
       IMPLICIT NONE  
 C  
 C--------------------------------------------------------------------  
26    
27  C global parameters updated by ggl90_calc  C global parameters updated by ggl90_calc
28  C     GGL90TKE     - sub-grid turbulent kinetic energy           (m^2/s^2)  C     GGL90TKE     :: sub-grid turbulent kinetic energy          (m^2/s^2)
29  C     GGL90viscAz  - GGL90 eddy viscosity coefficient              (m^2/s)  C     GGL90viscAz  :: GGL90 eddy viscosity coefficient             (m^2/s)
30  C     GGL90diffKzT - GGL90 diffusion coefficient for temperature   (m^2/s)  C     GGL90diffKzT :: GGL90 diffusion coefficient for temperature  (m^2/s)
 C  
31  C \ev  C \ev
32    
33  C !USES: ============================================================  C !USES: ============================================================
34          IMPLICIT NONE
35  #include "SIZE.h"  #include "SIZE.h"
36  #include "EEPARAMS.h"  #include "EEPARAMS.h"
37  #include "PARAMS.h"  #include "PARAMS.h"
# Line 43  C !USES: =============================== Line 41  C !USES: ===============================
41  #include "GRID.h"  #include "GRID.h"
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     myTime - Current time in simulation  C     sigmaR :: Vertical gradient of iso-neutral density
47    C     myTime :: Current time in simulation
48    C     myIter :: Current time-step number
49    C     myThid :: My Thread Id number
50        INTEGER bi, bj        INTEGER bi, bj
51        INTEGER myThid        _RL     sigmaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
52        _RL     myTime        _RL     myTime
53          INTEGER myIter
54          INTEGER myThid
55    CEOP
56    
57  #ifdef ALLOW_GGL90  #ifdef ALLOW_GGL90
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     buoyFreq         - buoyancy freqency  C     KappaH           :: (local) diffusivity parameter for temperature (eq.11)
72  C     TKEdissipation   - dissipation of TKE  C     KappaE           :: (local) diffusivity parameter for TKE (eq.15)
73  C     GGL90mixingLength- mixing length of scheme following Banke+Delecuse  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
84        _RL     Nsquare  c     _RL     Nsquare
85          _RL     Nsquare(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
86        _RL     deltaTggl90        _RL     deltaTggl90
87        _RL     SQRTTKE  c     _RL     SQRTTKE
88          _RL     SQRTTKE(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
89        _RL     RiNumber        _RL     RiNumber
90        _RL     TKEdissipation        _RL     TKEdissipation
91        _RL     tempU, tempV, prTemp        _RL     tempU, tempV, prTemp
92          _RL     MaxLength, tmpmlx, tmpVisc
93        _RL     TKEPrandtlNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     TKEPrandtlNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
94        _RL     GGL90mixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     GGL90mixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
95          _RL     rMixingLength    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
96          _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)
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     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
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  CEOP  #ifdef ALLOW_GGL90_SMOOTH
116          _RL p4, p8, p16
117          p4=0.25 _d 0
118          p8=0.125 _d 0
119          p16=0.0625 _d 0
120    #endif
121        iMin = 2-OLx        iMin = 2-OLx
122        iMax = sNx+OLx-1        iMax = sNx+OLx-1
123        jMin = 2-OLy        jMin = 2-OLy
# Line 111  CEOP Line 125  CEOP
125    
126  C     set separate time step (should be deltaTtracer)  C     set separate time step (should be deltaTtracer)
127        deltaTggl90 = dTtracerLev(1)        deltaTggl90 = dTtracerLev(1)
128  C      
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)  = 0. _d 0           GGL90mixingLength(i,j,k) = GGL90mixingLengthMin
141           GGL90mixingLength(I,J,K) = 0. _d 0           GGL90visctmp(i,j,k)      = 0. _d 0
142          ENDDO  #ifndef SOLVE_DIAGONAL_LOWMEMORY
143         ENDDO               a3d(i,j,k) = 0. _d 0
144        ENDDO           b3d(i,j,k) = 1. _d 0
145        DO J=1-Oly,sNy+Oly           c3d(i,j,k) = 0. _d 0
146         DO I=1-Olx,sNx+Olx  #endif
147          rhoK    (I,J) = 0. _d 0          ENDDO
148          rhoKm1  (I,J) = 0. _d 0         ENDDO
149          totalDepth(I,J)   = 0. _d 0        ENDDO
150          IF ( recip_Rcol(I,J,bi,bj) .NE. 0. )        DO j=1-OLy,sNy+OLy
151       &       totalDepth(I,J) = 1./recip_Rcol(I,J,bi,bj)         DO i=1-OLx,sNx+OLx
152            totalDepth(i,j)    = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
153            rMixingLength(i,j,1) = 0. _d 0
154            mxLength_Dn(i,j,1) = GGL90mixingLengthMin
155            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         Kp1 = MIN(Nr,K+1)  c      kp1 = MIN(Nr,k+1)
163         CALL FIND_RHO(         DO j=jMin,jMax
164       I      bi, bj, iMin, iMax, jMin, jMax, Km1, K,          DO i=iMin,iMax
165       I      theta, salt,           SQRTTKE(i,j,k)=SQRT( GGL90TKE(i,j,k,bi,bj) )
166       O      rhoKm1,  
      I      myThid )  
        CALL FIND_RHO(  
      I      bi, bj, iMin, iMax, jMin, jMax, K, K,  
      I      theta, salt,  
      O      rhoK,  
      I      myThid )  
        DO J=jMin,jMax  
         DO I=iMin,iMax  
          SQRTTKE=SQRT( GGL90TKE(I,J,K,bi,bj) )  
 C  
167  C     buoyancy frequency  C     buoyancy frequency
168  C           Nsquare(i,j,k) = gravity*gravitySign*recip_rhoConst
169           Nsquare = - gravity*recip_rhoConst*recip_drC(K)       &                  * sigmaR(i,j,k)
170       &        * ( rhoKm1(I,J) - rhoK(I,J) )*maskC(I,J,K,bi,bj)  cC     vertical shear term (dU/dz)^2+(dV/dz)^2
171  C     vertical shear term (dU/dz)^2+(dV/dz)^2  c         tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)
172           tempu= .5*(  uVel(I,J,Km1,bi,bj)+uVel(I+1,J,Km1,bi,bj)  c     &                 -( uVel(i,j,k  ,bi,bj)+uVel(i+1,j,k  ,bi,bj)) )
173       &             - (uVel(I,J,K  ,bi,bj)+uVel(I+1,J,K  ,bi,bj)) )  c     &        *recip_drC(k)
174       &        *recip_drC(K)  c         tempV= .5 _d 0*( vVel(i,j,km1,bi,bj)+vVel(i,j+1,km1,bi,bj)
175           tempv= .5*(  vVel(I,J,Km1,bi,bj)+vVel(I,J+1,Km1,bi,bj)  c     &                 -( vVel(i,j,k  ,bi,bj)+vVel(i,j+1,k  ,bi,bj)) )
176       &             - (vVel(I,J,K  ,bi,bj)+vVel(I,J+1,K  ,bi,bj)) )  c     &        *recip_drC(k)
177       &        *recip_drC(K)  c         verticalShear = tempU*tempU + tempV*tempV
178           verticalShear = tempU*tempU + tempV*tempV  c         RiNumber   = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)
179           RiNumber   = MAX(Nsquare,0. _d 0)/(verticalShear+GGL90eps)  cC     compute Prandtl number (always greater than 0)
180  C     compute Prandtl number (always greater than 0)  c         prTemp = 1. _d 0
181           prTemp = 1. _d 0  c         IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
182           IF ( RiNumber .GE. 0.2 ) prTemp = 5.0 * RiNumber  c         TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
          TKEPrandtlNumber(I,J,K) = MIN(10.,prTemp)  
183  C     mixing length  C     mixing length
184           GGL90mixingLength(I,J,K) =           GGL90mixingLength(i,j,k) = SQRTTWO *
185       &        SQRTTKE/SQRT( MAX(Nsquare,GGL90eps) )       &        SQRTTKE(i,j,k)/SQRT( MAX(Nsquare(i,j,k),GGL90eps) )
186  C     impose upper bound for mixing length (total depth)          ENDDO
          GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),  
      &        totalDepth(I,J))  
 C     impose minimum mixing length (to avoid division by zero)  
          GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),  
      &        GGL90mixingLengthMin)  
 C     viscosity of last timestep  
          KappaM = GGL90ck*GGL90mixingLength(I,J,K)*SQRTTKE  
          KappaE(I,J,K) = KappaM*GGL90alpha  
 C     dissipation term  
          TKEdissipation = ab05*GGL90ceps  
      &        *SQRTTKE/GGL90mixingLength(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  
      &        - KappaM*Nsquare/TKEPrandtlNumber(I,J,K)  
      &        - TKEdissipation  
      &        )  
         ENDDO    
187         ENDDO         ENDDO
188        ENDDO        ENDDO
189  C     horizontal diffusion of TKE (requires an exchange in  
190  C       do_fields_blocking_exchanges)  C- ensure mixing between first and second level
191          IF (mxlSurfFlag) THEN
192           DO j=jMin,jMax
193            DO i=iMin,iMax
194             GGL90mixingLength(i,j,2)=drF(1)
195            ENDDO
196           ENDDO
197          ENDIF
198    
199    C- Impose upper and lower bound for mixing length
200          IF ( mxlMaxFlag .EQ. 0 ) THEN
201    
202           DO k=2,Nr
203            DO j=jMin,jMax
204             DO i=iMin,iMax
205              MaxLength=totalDepth(i,j)
206              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
207         &                                   MaxLength)
208             ENDDO
209            ENDDO
210           ENDDO
211    
212           DO k=2,Nr
213            DO j=jMin,jMax
214             DO i=iMin,iMax
215              GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
216         &                                   GGL90mixingLengthMin)
217              rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
218             ENDDO
219            ENDDO
220           ENDDO
221    
222          ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN
223    
224           DO k=2,Nr
225            DO j=jMin,jMax
226             DO i=iMin,iMax
227              MaxLength=MIN(Ro_surf(i,j,bi,bj)-rF(k),rF(k)-R_low(i,j,bi,bj))
228    c         MaxLength=MAX(MaxLength,20. _d 0)
229              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
230         &                                   MaxLength)
231             ENDDO
232            ENDDO
233           ENDDO
234    
235           DO k=2,Nr
236            DO j=jMin,jMax
237             DO i=iMin,iMax
238              GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
239         &                                   GGL90mixingLengthMin)
240              rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
241             ENDDO
242            ENDDO
243           ENDDO
244    
245          ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN
246    
247           DO k=2,Nr
248            DO j=jMin,jMax
249             DO i=iMin,iMax
250              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
251         &        GGL90mixingLength(i,j,k-1)+drF(k-1))
252             ENDDO
253            ENDDO
254           ENDDO
255           DO j=jMin,jMax
256            DO i=iMin,iMax
257              GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
258         &       GGL90mixingLengthMin+drF(Nr))
259            ENDDO
260           ENDDO
261           DO k=Nr-1,2,-1
262            DO j=jMin,jMax
263             DO i=iMin,iMax
264              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
265         &        GGL90mixingLength(i,j,k+1)+drF(k))
266             ENDDO
267            ENDDO
268           ENDDO
269    
270           DO k=2,Nr
271            DO j=jMin,jMax
272             DO i=iMin,iMax
273              GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
274         &                                   GGL90mixingLengthMin)
275              rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
276             ENDDO
277            ENDDO
278           ENDDO
279    
280          ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN
281    
282           DO k=2,Nr
283            DO j=jMin,jMax
284             DO i=iMin,iMax
285              mxLength_Dn(i,j,k) = MIN(GGL90mixingLength(i,j,k),
286         &        mxLength_Dn(i,j,k-1)+drF(k-1))
287             ENDDO
288            ENDDO
289           ENDDO
290           DO j=jMin,jMax
291            DO i=iMin,iMax
292              GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
293         &       GGL90mixingLengthMin+drF(Nr))
294            ENDDO
295           ENDDO
296           DO k=Nr-1,2,-1
297            DO j=jMin,jMax
298             DO i=iMin,iMax
299              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
300         &        GGL90mixingLength(i,j,k+1)+drF(k))
301             ENDDO
302            ENDDO
303           ENDDO
304    
305           DO k=2,Nr
306            DO j=jMin,jMax
307             DO i=iMin,iMax
308              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
309         &                                  mxLength_Dn(i,j,k))
310              tmpmlx = SQRT( GGL90mixingLength(i,j,k)*mxLength_Dn(i,j,k) )
311              tmpmlx = MAX( tmpmlx, GGL90mixingLengthMin)
312              rMixingLength(i,j,k) = 1. _d 0 / tmpmlx
313             ENDDO
314            ENDDO
315           ENDDO
316    
317          ELSE
318           STOP 'GGL90_CALC: Wrong mxlMaxFlag (mixing length limit)'
319          ENDIF
320    
321    C- Impose minimum mixing length (to avoid division by zero)
322    c      DO k=2,Nr
323    c      DO j=jMin,jMax
324    c       DO i=iMin,iMax
325    c        GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
326    c    &        GGL90mixingLengthMin)
327    c        rMixingLength(i,j,k) = 1. _d 0 /GGL90mixingLength(i,j,k)
328    c       ENDDO
329    c      ENDDO
330    c     ENDDO
331    
332          DO k=2,Nr
333           km1 = k-1
334    
335  #ifdef ALLOW_GGL90_HORIZDIFF  #ifdef ALLOW_GGL90_HORIZDIFF
336        IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN         IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
337         DO K=2,Nr  C     horizontal diffusion of TKE (requires an exchange in
338    C      do_fields_blocking_exchanges)
339  C     common factors  C     common factors
340          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
341           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
342            xA(i,j) = _dyG(i,j,bi,bj)            xA(i,j) = _dyG(i,j,bi,bj)
343       &         *drF(k)*_hFacW(i,j,k,bi,bj)       &         *drF(k)*_hFacW(i,j,k,bi,bj)
344            yA(i,j) = _dxG(i,j,bi,bj)            yA(i,j) = _dxG(i,j,bi,bj)
345       &         *drF(k)*_hFacS(i,j,k,bi,bj)       &         *drF(k)*_hFacS(i,j,k,bi,bj)
346           ENDDO           ENDDO
347          ENDDO            ENDDO
348  C     Compute diffusive fluxes  C     Compute diffusive fluxes
349  C     ... across x-faces  C     ... across x-faces
350          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
351           dfx(1-Olx,j)=0.           dfx(1-OLx,j)=0. _d 0
352           DO i=1-Olx+1,sNx+Olx           DO i=1-OLx+1,sNx+OLx
353            dfx(i,j) = -GGL90diffTKEh*xA(i,j)            dfx(i,j) = -GGL90diffTKEh*xA(i,j)
354       &      *_recip_dxC(i,j,bi,bj)       &      *_recip_dxC(i,j,bi,bj)
355       &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i-1,j,k,bi,bj))       &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i-1,j,k,bi,bj))
# Line 227  C     ... across x-faces Line 357  C     ... across x-faces
357           ENDDO           ENDDO
358          ENDDO          ENDDO
359  C     ... across y-faces  C     ... across y-faces
360          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
361           dfy(i,1-Oly)=0.           dfy(i,1-OLy)=0. _d 0
362          ENDDO          ENDDO
363          DO j=1-Oly+1,sNy+Oly          DO j=1-OLy+1,sNy+OLy
364           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
365            dfy(i,j) = -GGL90diffTKEh*yA(i,j)            dfy(i,j) = -GGL90diffTKEh*yA(i,j)
366       &      *_recip_dyC(i,j,bi,bj)       &      *_recip_dyC(i,j,bi,bj)
367       &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i,j-1,k,bi,bj))       &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i,j-1,k,bi,bj))
# Line 239  C     ... across y-faces Line 369  C     ... across y-faces
369       &      *CosFacV(j,bi,bj)       &      *CosFacV(j,bi,bj)
370  #endif /* ISOTROPIC_COS_SCALING */  #endif /* ISOTROPIC_COS_SCALING */
371           ENDDO           ENDDO
372          ENDDO            ENDDO
373  C     Compute divergence of fluxes  C     Compute divergence of fluxes
374          DO j=1-Oly,sNy+Oly-1          DO j=1-OLy,sNy+OLy-1
375           DO i=1-Olx,sNx+Olx-1           DO i=1-OLx,sNx+OLx-1
376            gTKE(i,j,k)=gTKE(i,j,k)            gTKE(i,j) =
377       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj)       &    -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj)
378       &         *( (dfx(i+1,j)-dfx(i,j))       &         *( (dfx(i+1,j)-dfx(i,j))
379       &           +(dfy(i,j+1)-dfy(i,j))       &           +(dfy(i,j+1)-dfy(i,j))
380       &           )       &          )
381           ENDDO             ENDDO
382            ENDDO
383    C      end if GGL90diffTKEh .eq. 0.
384           ENDIF
385    #endif /* ALLOW_GGL90_HORIZDIFF */
386    
387           DO j=jMin,jMax
388            DO i=iMin,iMax
389    C     vertical shear term (dU/dz)^2+(dV/dz)^2
390             tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)
391         &                 -( uVel(i,j,k  ,bi,bj)+uVel(i+1,j,k  ,bi,bj)) )
392         &        *recip_drC(k)
393             tempV= .5 _d 0*( vVel(i,j,km1,bi,bj)+vVel(i,j+1,km1,bi,bj)
394         &                 -( vVel(i,j,k  ,bi,bj)+vVel(i,j+1,k  ,bi,bj)) )
395         &        *recip_drC(k)
396             verticalShear = tempU*tempU + tempV*tempV
397             RiNumber = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)
398    C     compute Prandtl number (always greater than 0)
399             prTemp = 1. _d 0
400             IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
401             TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
402    c         TKEPrandtlNumber(i,j,k) = 1. _d 0
403    
404    C     viscosity and diffusivity
405             KappaM = GGL90ck*GGL90mixingLength(i,j,k)*SQRTTKE(i,j,k)
406             GGL90visctmp(i,j,k) = MAX(KappaM,diffKrNrT(k))
407         &                            * maskC(i,j,k,bi,bj)
408    c        note: storing GGL90visctmp like this, and using it later to compute
409    c              GGL9rdiffKr etc. is robust in case of smoothing (e.g. see OPA)
410             KappaM = MAX(KappaM,viscArNr(k)) * maskC(i,j,k,bi,bj)
411             KappaH = KappaM/TKEPrandtlNumber(i,j,k)
412             KappaE(i,j,k) = GGL90alpha * KappaM * maskC(i,j,k,bi,bj)
413    
414    C     dissipation term
415             TKEdissipation = explDissFac*GGL90ceps
416         &        *SQRTTKE(i,j,k)*rMixingLength(i,j,k)
417         &        *GGL90TKE(i,j,k,bi,bj)
418    C     partial update with sum of explicit contributions
419             GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
420         &        + deltaTggl90*(
421         &        + KappaM*verticalShear
422         &        - KappaH*Nsquare(i,j,k)
423         &        - TKEdissipation
424         &        )
425          ENDDO          ENDDO
 C       end of k-loop  
426         ENDDO         ENDDO
427  C     end if GGL90diffTKEh .eq. 0.  
428        ENDIF  #ifdef ALLOW_GGL90_HORIZDIFF
429           IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
430    C--    Add horiz. diffusion tendency
431            DO j=jMin,jMax
432             DO i=iMin,iMax
433              GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
434         &                          + gTKE(i,j)*deltaTggl90
435             ENDDO
436            ENDDO
437           ENDIF
438  #endif /* ALLOW_GGL90_HORIZDIFF */  #endif /* ALLOW_GGL90_HORIZDIFF */
439    
440    C--   end of k loop
441          ENDDO
442    
443  C     ============================================  C     ============================================
444  C     Implicit time step to update TKE for k=1,Nr;  C     Implicit time step to update TKE for k=1,Nr;
445  C     TKE(Nr+1)=0 by default  C     TKE(Nr+1)=0 by default
# Line 264  C     set up matrix Line 448  C     set up matrix
448  C--   Lower diagonal  C--   Lower diagonal
449        DO j=jMin,jMax        DO j=jMin,jMax
450         DO i=iMin,iMax         DO i=iMin,iMax
451           a(i,j,1) = 0. _d 0           a3d(i,j,1) = 0. _d 0
452         ENDDO         ENDDO
453        ENDDO        ENDDO
454        DO k=2,Nr        DO k=2,Nr
455         km1=MAX(1,k-1)         km1=MAX(2,k-1)
456         DO j=jMin,jMax         DO j=jMin,jMax
457          DO i=iMin,iMax          DO i=iMin,iMax
458           a(i,j,k) = -deltaTggl90  C-    We keep recip_hFacC in the diffusive flux calculation,
459       &        *recip_drF(km1)*recip_hFacI(i,j,k,bi,bj)  C-    but no hFacC in TKE volume control
460       &        *.5*(KappaE(i,j, k )+KappaE(i,j,km1))  C-    No need for maskC(k-1) with recip_hFacC(k-1)
461       &        *recip_drC(k)           a3d(i,j,k) = -deltaTggl90
462            IF (recip_hFacI(i,j,km1,bi,bj).EQ.0.) a(i,j,k)=0.       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)
463         &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))
464         &        *recip_drC(k)*maskC(i,j,k,bi,bj)
465          ENDDO          ENDDO
466         ENDDO         ENDDO
467        ENDDO        ENDDO
468  C--   Upper diagonal  C--   Upper diagonal
469        DO j=jMin,jMax        DO j=jMin,jMax
470         DO i=iMin,iMax         DO i=iMin,iMax
471           c(i,j,1)  = 0. _d 0           c3d(i,j,1)  = 0. _d 0
          c(i,j,Nr) = 0. _d 0  
472         ENDDO         ENDDO
473        ENDDO        ENDDO
474  CML      DO k=1,Nr-1        DO k=2,Nr
       DO k=2,Nr-1  
        kp1=min(Nr,k+1)  
475         DO j=jMin,jMax         DO j=jMin,jMax
476          DO i=iMin,iMax          DO i=iMin,iMax
477            c(i,j,k) = -deltaTggl90            kp1=MAX(1,MIN(klowC(i,j,bi,bj),k+1))
478       &        *recip_drF( k )*recip_hFacI(i,j,k,bi,bj)  C-    We keep recip_hFacC in the diffusive flux calculation,
479       &               *.5*(KappaE(i,j,k)+KappaE(i,j,kp1))  C-    but no hFacC in TKE volume control
480       &        *recip_drC(k)  C-    No need for maskC(k) with recip_hFacC(k)
481            IF (recip_hFacI(i,j,kp1,bi,bj).EQ.0.) c(i,j,k)=0.            c3d(i,j,k) = -deltaTggl90
482         &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)
483         &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))
484         &        *recip_drC(k)*maskC(i,j,k-1,bi,bj)
485          ENDDO          ENDDO
486         ENDDO         ENDDO
487        ENDDO        ENDDO
488  C--   Center diagonal  C--   Center diagonal
489        DO k=1,Nr        DO k=1,Nr
490           km1 = MAX(k-1,1)
491         DO j=jMin,jMax         DO j=jMin,jMax
492          DO i=iMin,iMax          DO i=iMin,iMax
493            b(i,j,k) = 1. _d 0 - c(i,j,k) - a(i,j,k)            b3d(i,j,k) = 1. _d 0 - c3d(i,j,k) - a3d(i,j,k)
494       &        + ab15*deltaTggl90*GGL90ceps*SQRT(GGL90TKE(I,J,K,bi,bj))       &        + implDissFac*deltaTggl90*GGL90ceps*SQRTTKE(i,j,k)
495       &        /GGL90mixingLength(I,J,K)       &        * rMixingLength(i,j,k)
496         &        * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)
497           ENDDO           ENDDO
498         ENDDO         ENDDO
499        ENDDO        ENDDO
500  C     end set up matrix  C     end set up matrix
501    
 C  
502  C     Apply boundary condition  C     Apply boundary condition
503  C            kp1 = MIN(Nr,kSurf+1)
504        DO J=jMin,jMax        DO j=jMin,jMax
505         DO I=iMin,iMax         DO i=iMin,iMax
506  C     estimate friction velocity uStar from surface forcing  C     estimate friction velocity uStar from surface forcing
507          uStarSquare = SQRT(          uStarSquare = SQRT(
508       &         ( .5*( surfaceForcingU(I,  J,  bi,bj)       &    ( .5 _d 0*( surfaceForcingU(i,  j,  bi,bj)
509       &              + surfaceForcingU(I+1,J,  bi,bj) ) )**2       &              + surfaceForcingU(i+1,j,  bi,bj) ) )**2
510       &       + ( .5*( surfaceForcingV(I,  J,  bi,bj)       &  + ( .5 _d 0*( surfaceForcingV(i,  j,  bi,bj)
511       &              + surfaceForcingV(I,  J+1,bi,bj) ) )**2       &              + surfaceForcingV(i,  j+1,bi,bj) ) )**2
512       &                     )       &                     )
513  C     Dirichlet surface boundary condition for TKE  C     Dirichlet surface boundary condition for TKE
514          gTKE(I,J,kSurf) = MAX(GGL90TKEmin,GGL90m2*uStarSquare)          GGL90TKE(i,j,kSurf,bi,bj) = maskC(i,j,kSurf,bi,bj)
515       &                     *maskC(I,J,kSurf,bi,bj)       &           *MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)
516            GGL90TKE(i,j,kp1,bi,bj) = GGL90TKE(i,j,kp1,bi,bj)
517         &               - a3d(i,j,kp1)*GGL90TKE(i,j,kSurf,bi,bj)
518            a3d(i,j,kp1) = 0. _d 0
519  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom
520          kBottom   = MIN(MAX(kLowC(I,J,bi,bj),1),Nr)          kBottom   = MAX(kLowC(i,j,bi,bj),1)
521          gTKE(I,J,kBottom) = gTKE(I,J,kBottom)          GGL90TKE(i,j,kBottom,bi,bj) = GGL90TKE(i,j,kBottom,bi,bj)
522       &       - GGL90TKEbottom*c(I,J,kBottom)       &                              - GGL90TKEbottom*c3d(i,j,kBottom)
523         ENDDO          c3d(i,j,kBottom) = 0. _d 0
524        ENDDO             ENDDO
525  C        ENDDO
526  C     solve tri-diagonal system, and store solution on gTKE (previously rhs)  
527  C  C     solve tri-diagonal system
528        CALL GGL90_SOLVE( bi, bj, iMin, iMax, jMin, jMax,        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
529       I     a, b, c,       I                        a3d, b3d, c3d,
530       U     gTKE,       U                        GGL90TKE,
531       I     myThid )       O                        errCode,
532  C       I                        bi, bj, myThid )
533  C     now update TKE  
534  C            DO k=1,Nr
535        DO K=1,Nr         DO j=jMin,jMax
536         DO J=jMin,jMax          DO i=iMin,iMax
         DO I=iMin,iMax  
537  C     impose minimum TKE to avoid numerical undershoots below zero  C     impose minimum TKE to avoid numerical undershoots below zero
538           GGL90TKE(I,J,K,bi,bj) = MAX( gTKE(I,J,K), GGL90TKEmin )           GGL90TKE(i,j,k,bi,bj) = maskC(i,j,k,bi,bj)
539       &        * maskC(I,J,K,bi,bj)       &                  *MAX( GGL90TKE(i,j,k,bi,bj), GGL90TKEmin )
540          ENDDO          ENDDO
541         ENDDO         ENDDO
542        ENDDO            ENDDO
543  C  
544  C     end of time step  C     end of time step
545  C     ===============================  C     ===============================
546  C     compute viscosity coefficients  
547  C            DO k=2,Nr
548        DO K=2,Nr         DO j=1,sNy
549         DO J=jMin,jMax          DO i=1,sNx
550          DO I=iMin,iMax  #ifdef ALLOW_GGL90_SMOOTH
551  C     Eq. (11), (18) and (21)           tmpVisc=
552           KappaM = GGL90ck*GGL90mixingLength(I,J,K)*       &  (
553       &                  SQRT( GGL90TKE(I,J,K,bi,bj) )       &   p4 *  GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
554           KappaH = KappaM/TKEPrandtlNumber(I,J,K)       &  +p8 *( GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj)
555  C     Set a minium (= background) value       &       + GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj)
556           KappaM = MAX(KappaM,viscAr)       &       + GGL90visctmp(i+1,j  ,k) * mskCor(i+1,j  ,bi,bj)
557           KappaH = MAX(KappaH,diffKrNrT(k))       &       + GGL90visctmp(i  ,j+1,k) * mskCor(i  ,j+1,bi,bj))
558  C     Set a maximum and mask land point       &  +p16*( GGL90visctmp(i+1,j+1,k) * mskCor(i+1,j+1,bi,bj)
559           GGL90viscAr(I,J,K,bi,bj) = MIN(KappaM,GGL90viscMax)       &       + GGL90visctmp(i+1,j-1,k) * mskCor(i+1,j-1,bi,bj)
560       &        * maskC(I,J,K,bi,bj)       &       + GGL90visctmp(i-1,j+1,k) * mskCor(i-1,j+1,bi,bj)
561           GGL90diffKr(I,J,K,bi,bj) = MIN(KappaH,GGL90diffMax)       &       + GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj))
562       &        * maskC(I,J,K,bi,bj)       &  )
563          ENDDO       & /(p4
564         ENDDO       &  +p8 *(       maskC(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)
565  C     end third k-loop       &       +       maskC(i  ,j-1,k,bi,bj) * mskCor(i  ,j-1,bi,bj)
566        ENDDO           &       +       maskC(i+1,j  ,k,bi,bj) * mskCor(i+1,j  ,bi,bj)
567         &       +       maskC(i  ,j+1,k,bi,bj) * mskCor(i  ,j+1,bi,bj))
568         &  +p16*(       maskC(i+1,j+1,k,bi,bj) * mskCor(i+1,j+1,bi,bj)
569         &       +       maskC(i+1,j-1,k,bi,bj) * mskCor(i+1,j-1,bi,bj)
570         &       +       maskC(i-1,j+1,k,bi,bj) * mskCor(i-1,j+1,bi,bj)
571         &       +       maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj))
572         &  )*maskC(i,j,k,bi,bj)*mskCor(i,j,bi,bj)
573    #else
574             tmpVisc = GGL90visctmp(i,j,k)
575    #endif
576             tmpVisc = MIN(tmpVisc/TKEPrandtlNumber(i,j,k),GGL90diffMax)
577             GGL90diffKr(i,j,k,bi,bj)= MAX( tmpVisc , diffKrNrT(k) )
578            ENDDO
579           ENDDO
580          ENDDO
581    
582          DO k=2,Nr
583           DO j=1,sNy
584            DO i=1,sNx+1
585    #ifdef ALLOW_GGL90_SMOOTH
586            tmpVisc =
587         & (
588         &   p4 *(GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
589         &       +GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj))
590         &  +p8 *(GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj)
591         &       +GGL90visctmp(i-1,j+1,k) * mskCor(i-1,j+1,bi,bj)
592         &       +GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj)
593         &       +GGL90visctmp(i  ,j+1,k) * mskCor(i  ,j+1,bi,bj))
594         &  )
595         & /(p4 * 2. _d 0
596         &  +p8 *(      maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj)
597         &       +      maskC(i-1,j+1,k,bi,bj) * mskCor(i-1,j+1,bi,bj)
598         &       +      maskC(i  ,j-1,k,bi,bj) * mskCor(i  ,j-1,bi,bj)
599         &       +      maskC(i  ,j+1,k,bi,bj) * mskCor(i  ,j+1,bi,bj))
600         &  )
601         &  *maskC(i  ,j,k,bi,bj)*mskCor(i  ,j,bi,bj)
602         &  *maskC(i-1,j,k,bi,bj)*mskCor(i-1,j,bi,bj)
603    #else
604            tmpVisc = _maskW(i,j,k,bi,bj) *
605         &                   (.5 _d 0*(GGL90visctmp(i,j,k)
606         &                            +GGL90visctmp(i-1,j,k))
607         &                   )
608    #endif
609            tmpVisc = MIN( tmpVisc , GGL90viscMax )
610            GGL90viscArU(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
611            ENDDO
612           ENDDO
613          ENDDO
614    
615          DO k=2,Nr
616           DO j=1,sNy+1
617            DO i=1,sNx
618    #ifdef ALLOW_GGL90_SMOOTH
619            tmpVisc =
620         & (
621         &   p4 *(GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
622         &       +GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj))
623         &  +p8 *(GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj)
624         &       +GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj)
625         &       +GGL90visctmp(i+1,j  ,k) * mskCor(i+1,j  ,bi,bj)
626         &       +GGL90visctmp(i+1,j-1,k) * mskCor(i+1,j-1,bi,bj))
627         &  )
628         & /(p4 * 2. _d 0
629         &  +p8 *(      maskC(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)
630         &       +      maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj)
631         &       +      maskC(i+1,j  ,k,bi,bj) * mskCor(i+1,j  ,bi,bj)
632         &       +      maskC(i+1,j-1,k,bi,bj) * mskCor(i+1,j-1,bi,bj))
633         &  )
634         &   *maskC(i,j  ,k,bi,bj)*mskCor(i,j  ,bi,bj)
635         &   *maskC(i,j-1,k,bi,bj)*mskCor(i,j-1,bi,bj)
636    #else
637            tmpVisc = _maskS(i,j,k,bi,bj) *
638         &                   (.5 _d 0*(GGL90visctmp(i,j,k)
639         &                            +GGL90visctmp(i,j-1,k))
640         &                   )
641    
642    #endif
643            tmpVisc = MIN( tmpVisc , GGL90viscMax )
644            GGL90viscArV(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
645            ENDDO
646           ENDDO
647          ENDDO
648    
649    #ifdef ALLOW_DIAGNOSTICS
650          IF ( useDiagnostics ) THEN
651             CALL DIAGNOSTICS_FILL( GGL90TKE   ,'GGL90TKE',
652         &                          0,Nr, 1, bi, bj, myThid )
653             CALL DIAGNOSTICS_FILL( GGL90viscArU,'GGL90ArU',
654         &                          0,Nr, 1, bi, bj, myThid )
655             CALL DIAGNOSTICS_FILL( GGL90viscArV,'GGL90ArV',
656         &                          0,Nr, 1, bi, bj, myThid )
657             CALL DIAGNOSTICS_FILL( GGL90diffKr,'GGL90Kr ',
658         &                          0,Nr, 1, bi, bj, myThid )
659             CALL DIAGNOSTICS_FILL( TKEPrandtlNumber ,'GGL90Prl',
660         &                          0,Nr, 2, bi, bj, myThid )
661             CALL DIAGNOSTICS_FILL( GGL90mixingLength,'GGL90Lmx',
662         &                          0,Nr, 2, bi, bj, myThid )
663          ENDIF
664    #endif
665    
666  #endif /* ALLOW_GGL90 */  #endif /* ALLOW_GGL90 */
667    
668        RETURN        RETURN
669        END        END
   
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.23
  ViewVC Help
Powered by ViewVC 1.1.22