/[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.2 by mlosch, Mon Sep 27 08:02:04 2004 UTC revision 1.22 by jmc, Thu Jun 28 15:35:52 2012 UTC
# Line 7  CBOP Line 7  CBOP
7  C !ROUTINE: GGL90_CALC  C !ROUTINE: GGL90_CALC
8    
9  C !INTERFACE: ======================================================  C !INTERFACE: ======================================================
10        subroutine GGL90_CALC(        SUBROUTINE GGL90_CALC(
11       I     bi, bj, myTime, myThid )       I                 bi, bj, sigmaR, myTime, myIter, myThid )
12    
13    
14  C !DESCRIPTION: \bv  C !DESCRIPTION: \bv
15  C     /==========================================================\  C     *==========================================================*
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
124        jMax = sNy+OLy-1        jMax = sNy+OLy-1
125    
126  C     set separate time step (should be deltaTtracer)  C     set separate time step (should be deltaTtracer)
127        deltaTggl90 = deltaTtracer        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- Impose upper and lower bound for mixing length
191          IF ( mxlMaxFlag .EQ. 0 ) THEN
192    
193           DO k=2,Nr
194            DO j=jMin,jMax
195             DO i=iMin,iMax
196              MaxLength=totalDepth(i,j)
197              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
198         &                                   MaxLength)
199             ENDDO
200            ENDDO
201           ENDDO
202    
203           DO k=2,Nr
204            DO j=jMin,jMax
205             DO i=iMin,iMax
206              GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
207         &                                   GGL90mixingLengthMin)
208              rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
209             ENDDO
210            ENDDO
211           ENDDO
212    
213          ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN
214    
215           DO k=2,Nr
216            DO j=jMin,jMax
217             DO i=iMin,iMax
218              MaxLength=MIN(Ro_surf(i,j,bi,bj)-rF(k),rF(k)-R_low(i,j,bi,bj))
219    c         MaxLength=MAX(MaxLength,20. _d 0)
220              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
221         &                                   MaxLength)
222             ENDDO
223            ENDDO
224           ENDDO
225    
226           DO k=2,Nr
227            DO j=jMin,jMax
228             DO i=iMin,iMax
229              GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
230         &                                   GGL90mixingLengthMin)
231              rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
232             ENDDO
233            ENDDO
234           ENDDO
235    
236          ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN
237    
238    cgf ensure mixing between first and second level
239    c      DO j=jMin,jMax
240    c        DO i=iMin,iMax
241    c         GGL90mixingLength(i,j,2)=drF(1)
242    c        ENDDO
243    c      ENDDO
244    cgf
245           DO k=2,Nr
246            DO j=jMin,jMax
247             DO i=iMin,iMax
248              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
249         &        GGL90mixingLength(i,j,k-1)+drF(k-1))
250             ENDDO
251            ENDDO
252           ENDDO
253           DO j=jMin,jMax
254            DO i=iMin,iMax
255              GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
256         &       GGL90mixingLengthMin+drF(Nr))
257            ENDDO
258           ENDDO
259           DO k=Nr-1,2,-1
260            DO j=jMin,jMax
261             DO i=iMin,iMax
262              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
263         &        GGL90mixingLength(i,j,k+1)+drF(k))
264             ENDDO
265            ENDDO
266           ENDDO
267    
268           DO k=2,Nr
269            DO j=jMin,jMax
270             DO i=iMin,iMax
271              GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
272         &                                   GGL90mixingLengthMin)
273              rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
274             ENDDO
275            ENDDO
276           ENDDO
277    
278          ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN
279    
280           DO k=2,Nr
281            DO j=jMin,jMax
282             DO i=iMin,iMax
283              mxLength_Dn(i,j,k) = MIN(GGL90mixingLength(i,j,k),
284         &        mxLength_Dn(i,j,k-1)+drF(k-1))
285             ENDDO
286            ENDDO
287           ENDDO
288           DO j=jMin,jMax
289            DO i=iMin,iMax
290              GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
291         &       GGL90mixingLengthMin+drF(Nr))
292            ENDDO
293           ENDDO
294           DO k=Nr-1,2,-1
295            DO j=jMin,jMax
296             DO i=iMin,iMax
297              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
298         &        GGL90mixingLength(i,j,k+1)+drF(k))
299             ENDDO
300            ENDDO
301           ENDDO
302    
303           DO k=2,Nr
304            DO j=jMin,jMax
305             DO i=iMin,iMax
306              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
307         &                                  mxLength_Dn(i,j,k))
308              tmpmlx = SQRT( GGL90mixingLength(i,j,k)*mxLength_Dn(i,j,k) )
309              tmpmlx = MAX( tmpmlx, GGL90mixingLengthMin)
310              rMixingLength(i,j,k) = 1. _d 0 / tmpmlx
311             ENDDO
312            ENDDO
313           ENDDO
314    
315          ELSE
316           STOP 'GGL90_CALC: Wrong mxlMaxFlag (mixing length limit)'
317          ENDIF
318    
319    C- Impose minimum mixing length (to avoid division by zero)
320    c      DO k=2,Nr
321    c      DO j=jMin,jMax
322    c       DO i=iMin,iMax
323    c        GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
324    c    &        GGL90mixingLengthMin)
325    c        rMixingLength(i,j,k) = 1. _d 0 /GGL90mixingLength(i,j,k)
326    c       ENDDO
327    c      ENDDO
328    c     ENDDO
329    
330          DO k=2,Nr
331           km1 = k-1
332    
333  #ifdef ALLOW_GGL90_HORIZDIFF  #ifdef ALLOW_GGL90_HORIZDIFF
334        IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN         IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
335         DO K=2,Nr  C     horizontal diffusion of TKE (requires an exchange in
336    C      do_fields_blocking_exchanges)
337  C     common factors  C     common factors
338          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
339           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
340            xA(i,j) = _dyG(i,j,bi,bj)            xA(i,j) = _dyG(i,j,bi,bj)
341       &         *drF(k)*_hFacW(i,j,k,bi,bj)       &         *drF(k)*_hFacW(i,j,k,bi,bj)
342            yA(i,j) = _dxG(i,j,bi,bj)            yA(i,j) = _dxG(i,j,bi,bj)
343       &         *drF(k)*_hFacS(i,j,k,bi,bj)       &         *drF(k)*_hFacS(i,j,k,bi,bj)
344           ENDDO           ENDDO
345          ENDDO            ENDDO
346  C     Compute diffusive fluxes  C     Compute diffusive fluxes
347  C     ... across x-faces  C     ... across x-faces
348          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
349           dfx(1-Olx,j)=0.           dfx(1-OLx,j)=0. _d 0
350           DO i=1-Olx+1,sNx+Olx           DO i=1-OLx+1,sNx+OLx
351            dfx(i,j) = -GGL90diffTKEh*xA(i,j)            dfx(i,j) = -GGL90diffTKEh*xA(i,j)
352       &      *_recip_dxC(i,j,bi,bj)       &      *_recip_dxC(i,j,bi,bj)
353       &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i-1,j,k,bi,bj))       &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i-1,j,k,bi,bj))
# Line 227  C     ... across x-faces Line 355  C     ... across x-faces
355           ENDDO           ENDDO
356          ENDDO          ENDDO
357  C     ... across y-faces  C     ... across y-faces
358          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
359           dfy(i,1-Oly)=0.           dfy(i,1-OLy)=0. _d 0
360          ENDDO          ENDDO
361          DO j=1-Oly+1,sNy+Oly          DO j=1-OLy+1,sNy+OLy
362           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
363            dfy(i,j) = -GGL90diffTKEh*yA(i,j)            dfy(i,j) = -GGL90diffTKEh*yA(i,j)
364       &      *_recip_dyC(i,j,bi,bj)       &      *_recip_dyC(i,j,bi,bj)
365       &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i,j-1,k,bi,bj))       &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i,j-1,k,bi,bj))
# Line 239  C     ... across y-faces Line 367  C     ... across y-faces
367       &      *CosFacV(j,bi,bj)       &      *CosFacV(j,bi,bj)
368  #endif /* ISOTROPIC_COS_SCALING */  #endif /* ISOTROPIC_COS_SCALING */
369           ENDDO           ENDDO
370          ENDDO            ENDDO
371  C     Compute divergence of fluxes  C     Compute divergence of fluxes
372          DO j=1-Oly,sNy+Oly-1          DO j=1-OLy,sNy+OLy-1
373           DO i=1-Olx,sNx+Olx-1           DO i=1-OLx,sNx+OLx-1
374            gTKE(i,j,k)=gTKE(i,j,k)            gTKE(i,j) =
375       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj)       &    -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj)
376       &         *( (dfx(i+1,j)-dfx(i,j))       &         *( (dfx(i+1,j)-dfx(i,j))
377       &           +(dfy(i,j+1)-dfy(i,j))       &           +(dfy(i,j+1)-dfy(i,j))
378       &           )       &          )
379           ENDDO             ENDDO
380            ENDDO
381    C      end if GGL90diffTKEh .eq. 0.
382           ENDIF
383    #endif /* ALLOW_GGL90_HORIZDIFF */
384    
385           DO j=jMin,jMax
386            DO i=iMin,iMax
387    C     vertical shear term (dU/dz)^2+(dV/dz)^2
388             tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)
389         &                 -( uVel(i,j,k  ,bi,bj)+uVel(i+1,j,k  ,bi,bj)) )
390         &        *recip_drC(k)
391             tempV= .5 _d 0*( vVel(i,j,km1,bi,bj)+vVel(i,j+1,km1,bi,bj)
392         &                 -( vVel(i,j,k  ,bi,bj)+vVel(i,j+1,k  ,bi,bj)) )
393         &        *recip_drC(k)
394             verticalShear = tempU*tempU + tempV*tempV
395             RiNumber = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)
396    C     compute Prandtl number (always greater than 0)
397             prTemp = 1. _d 0
398             IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
399             TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
400    c         TKEPrandtlNumber(i,j,k) = 1. _d 0
401    
402    C     viscosity and diffusivity
403             KappaM = GGL90ck*GGL90mixingLength(i,j,k)*SQRTTKE(i,j,k)
404             GGL90visctmp(i,j,k) = MAX(KappaM,diffKrNrT(k))
405         &                            * maskC(i,j,k,bi,bj)
406    c        note: storing GGL90visctmp like this, and using it later to compute
407    c              GGL9rdiffKr etc. is robust in case of smoothing (e.g. see OPA)
408             KappaM = MAX(KappaM,viscArNr(k)) * maskC(i,j,k,bi,bj)
409             KappaH = KappaM/TKEPrandtlNumber(i,j,k)
410             KappaE(i,j,k) = GGL90alpha * KappaM * maskC(i,j,k,bi,bj)
411    
412    C     dissipation term
413             TKEdissipation = explDissFac*GGL90ceps
414         &        *SQRTTKE(i,j,k)*rMixingLength(i,j,k)
415         &        *GGL90TKE(i,j,k,bi,bj)
416    C     partial update with sum of explicit contributions
417             GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
418         &        + deltaTggl90*(
419         &        + KappaM*verticalShear
420         &        - KappaH*Nsquare(i,j,k)
421         &        - TKEdissipation
422         &        )
423          ENDDO          ENDDO
 C       end of k-loop  
424         ENDDO         ENDDO
425  C     end if GGL90diffTKEh .eq. 0.  
426        ENDIF  #ifdef ALLOW_GGL90_HORIZDIFF
427           IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
428    C--    Add horiz. diffusion tendency
429            DO j=jMin,jMax
430             DO i=iMin,iMax
431              GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
432         &                          + gTKE(i,j)*deltaTggl90
433             ENDDO
434            ENDDO
435           ENDIF
436  #endif /* ALLOW_GGL90_HORIZDIFF */  #endif /* ALLOW_GGL90_HORIZDIFF */
437    
438    C--   end of k loop
439          ENDDO
440    
441  C     ============================================  C     ============================================
442  C     Implicit time step to update TKE for k=1,Nr;  C     Implicit time step to update TKE for k=1,Nr;
443  C     TKE(Nr+1)=0 by default  C     TKE(Nr+1)=0 by default
# Line 264  C     set up matrix Line 446  C     set up matrix
446  C--   Lower diagonal  C--   Lower diagonal
447        DO j=jMin,jMax        DO j=jMin,jMax
448         DO i=iMin,iMax         DO i=iMin,iMax
449           a(i,j,1) = 0. _d 0           a3d(i,j,1) = 0. _d 0
450         ENDDO         ENDDO
451        ENDDO        ENDDO
452        DO k=2,Nr        DO k=2,Nr
453         km1=MAX(1,k-1)         km1=MAX(2,k-1)
454         DO j=jMin,jMax         DO j=jMin,jMax
455          DO i=iMin,iMax          DO i=iMin,iMax
456           a(i,j,k) = -deltaTggl90  C-    We keep recip_hFacC in the diffusive flux calculation,
457       &        *recip_drF(km1)*recip_hFacI(i,j,k,bi,bj)  C-    but no hFacC in TKE volume control
458       &        *.5*(KappaE(i,j, k )+KappaE(i,j,km1))  C-    No need for maskC(k-1) with recip_hFacC(k-1)
459       &        *recip_drC(k)           a3d(i,j,k) = -deltaTggl90
460            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)
461         &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))
462         &        *recip_drC(k)*maskC(i,j,k,bi,bj)
463          ENDDO          ENDDO
464         ENDDO         ENDDO
465        ENDDO        ENDDO
466  C--   Upper diagonal  C--   Upper diagonal
467        DO j=jMin,jMax        DO j=jMin,jMax
468         DO i=iMin,iMax         DO i=iMin,iMax
469           c(i,j,1)  = 0. _d 0           c3d(i,j,1)  = 0. _d 0
          c(i,j,Nr) = 0. _d 0  
470         ENDDO         ENDDO
471        ENDDO        ENDDO
472  CML      DO k=1,Nr-1        DO k=2,Nr
       DO k=2,Nr-1  
        kp1=min(Nr,k+1)  
473         DO j=jMin,jMax         DO j=jMin,jMax
474          DO i=iMin,iMax          DO i=iMin,iMax
475            c(i,j,k) = -deltaTggl90            kp1=MAX(1,MIN(klowC(i,j,bi,bj),k+1))
476       &        *recip_drF( k )*recip_hFacI(i,j,k,bi,bj)  C-    We keep recip_hFacC in the diffusive flux calculation,
477       &               *.5*(KappaE(i,j,k)+KappaE(i,j,kp1))  C-    but no hFacC in TKE volume control
478       &        *recip_drC(k)  C-    No need for maskC(k) with recip_hFacC(k)
479            IF (recip_hFacI(i,j,kp1,bi,bj).EQ.0.) c(i,j,k)=0.            c3d(i,j,k) = -deltaTggl90
480         &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)
481         &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))
482         &        *recip_drC(k)*maskC(i,j,k-1,bi,bj)
483          ENDDO          ENDDO
484         ENDDO         ENDDO
485        ENDDO        ENDDO
486  C--   Center diagonal  C--   Center diagonal
487        DO k=1,Nr        DO k=1,Nr
488           km1 = MAX(k-1,1)
489         DO j=jMin,jMax         DO j=jMin,jMax
490          DO i=iMin,iMax          DO i=iMin,iMax
491            b(i,j,k) = 1. _d 0 - c(i,j,k) - a(i,j,k)            b3d(i,j,k) = 1. _d 0 - c3d(i,j,k) - a3d(i,j,k)
492       &        + ab15*deltaTggl90*GGL90ceps*SQRT(GGL90TKE(I,J,K,bi,bj))       &        + implDissFac*deltaTggl90*GGL90ceps*SQRTTKE(i,j,k)
493       &        /GGL90mixingLength(I,J,K)       &        * rMixingLength(i,j,k)
494         &        * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)
495           ENDDO           ENDDO
496         ENDDO         ENDDO
497        ENDDO        ENDDO
498  C     end set up matrix  C     end set up matrix
499    
 C  
500  C     Apply boundary condition  C     Apply boundary condition
501  C            kp1 = MIN(Nr,kSurf+1)
502        DO J=jMin,jMax        DO j=jMin,jMax
503         DO I=iMin,iMax         DO i=iMin,iMax
504  C     estimate friction velocity uStar from surface forcing  C     estimate friction velocity uStar from surface forcing
505          uStarSquare = SQRT(          uStarSquare = SQRT(
506       &         ( .5*( surfaceForcingU(I,  J,  bi,bj)       &    ( .5 _d 0*( surfaceForcingU(i,  j,  bi,bj)
507       &              + surfaceForcingU(I+1,J,  bi,bj) ) )**2       &              + surfaceForcingU(i+1,j,  bi,bj) ) )**2
508       &       + ( .5*( surfaceForcingV(I,  J,  bi,bj)       &  + ( .5 _d 0*( surfaceForcingV(i,  j,  bi,bj)
509       &              + surfaceForcingV(I,  J+1,bi,bj) ) )**2       &              + surfaceForcingV(i,  j+1,bi,bj) ) )**2
510       &                     )       &                     )
511  C     Dirichlet surface boundary condition for TKE  C     Dirichlet surface boundary condition for TKE
512          gTKE(I,J,kSurf) = MAX(GGL90TKEmin,GGL90m2*uStarSquare)          GGL90TKE(i,j,kSurf,bi,bj) = maskC(i,j,kSurf,bi,bj)
513       &                     *maskC(I,J,kSurf,bi,bj)       &           *MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)
514            GGL90TKE(i,j,kp1,bi,bj) = GGL90TKE(i,j,kp1,bi,bj)
515         &               - a3d(i,j,kp1)*GGL90TKE(i,j,kSurf,bi,bj)
516            a3d(i,j,kp1) = 0. _d 0
517  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom
518          kBottom   = MIN(MAX(kLowC(I,J,bi,bj),1),Nr)          kBottom   = MAX(kLowC(i,j,bi,bj),1)
519          gTKE(I,J,kBottom) = gTKE(I,J,kBottom)          GGL90TKE(i,j,kBottom,bi,bj) = GGL90TKE(i,j,kBottom,bi,bj)
520       &       - GGL90TKEbottom*c(I,J,kBottom)       &                              - GGL90TKEbottom*c3d(i,j,kBottom)
521         ENDDO          c3d(i,j,kBottom) = 0. _d 0
522        ENDDO             ENDDO
523  C        ENDDO
524  C     solve tri-diagonal system, and store solution on gTKE (previously rhs)  
525  C  C     solve tri-diagonal system
526        CALL GGL90_SOLVE( bi, bj, iMin, iMax, jMin, jMax,        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
527       I     a, b, c,       I                        a3d, b3d, c3d,
528       U     gTKE,       U                        GGL90TKE,
529       I     myThid )       O                        errCode,
530  C       I                        bi, bj, myThid )
531  C     now update TKE  
532  C            DO k=1,Nr
533        DO K=1,Nr         DO j=jMin,jMax
534         DO J=jMin,jMax          DO i=iMin,iMax
         DO I=iMin,iMax  
535  C     impose minimum TKE to avoid numerical undershoots below zero  C     impose minimum TKE to avoid numerical undershoots below zero
536           GGL90TKE(I,J,K,bi,bj) = MAX( gTKE(I,J,K), GGL90TKEmin )           GGL90TKE(i,j,k,bi,bj) = maskC(i,j,k,bi,bj)
537       &        * maskC(I,J,K,bi,bj)       &                  *MAX( GGL90TKE(i,j,k,bi,bj), GGL90TKEmin )
538          ENDDO          ENDDO
539         ENDDO         ENDDO
540        ENDDO            ENDDO
541  C  
542  C     end of time step  C     end of time step
543  C     ===============================  C     ===============================
544  C     compute viscosity coefficients  
545  C            DO k=2,Nr
546        DO K=2,Nr         DO j=1,sNy
547         DO J=jMin,jMax          DO i=1,sNx
548          DO I=iMin,iMax  #ifdef ALLOW_GGL90_SMOOTH
549  C     Eq. (11), (18) and (21)           tmpVisc=
550           KappaM = GGL90ck*GGL90mixingLength(I,J,K)*       &  (
551       &                  SQRT( GGL90TKE(I,J,K,bi,bj) )       &   p4 *  GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
552           KappaH = KappaM/TKEPrandtlNumber(I,J,K)       &  +p8 *( GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj)
553  C     Set a minium (= background) value       &       + GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj)
554           KappaM = MAX(KappaM,viscAr)       &       + GGL90visctmp(i+1,j  ,k) * mskCor(i+1,j  ,bi,bj)
555           KappaH = MAX(KappaH,diffKrT)       &       + GGL90visctmp(i  ,j+1,k) * mskCor(i  ,j+1,bi,bj))
556  C     Set a maximum and mask land point       &  +p16*( GGL90visctmp(i+1,j+1,k) * mskCor(i+1,j+1,bi,bj)
557           GGL90viscAr(I,J,K,bi,bj) = MIN(KappaM,GGL90viscMax)       &       + GGL90visctmp(i+1,j-1,k) * mskCor(i+1,j-1,bi,bj)
558       &        * maskC(I,J,K,bi,bj)       &       + GGL90visctmp(i-1,j+1,k) * mskCor(i-1,j+1,bi,bj)
559           GGL90diffKr(I,J,K,bi,bj) = MIN(KappaH,GGL90diffMax)       &       + GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj))
560       &        * maskC(I,J,K,bi,bj)       &  )
561          ENDDO       & /(p4
562         ENDDO       &  +p8 *(       maskC(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)
563  C     end third k-loop       &       +       maskC(i  ,j-1,k,bi,bj) * mskCor(i  ,j-1,bi,bj)
564        ENDDO           &       +       maskC(i+1,j  ,k,bi,bj) * mskCor(i+1,j  ,bi,bj)
565         &       +       maskC(i  ,j+1,k,bi,bj) * mskCor(i  ,j+1,bi,bj))
566         &  +p16*(       maskC(i+1,j+1,k,bi,bj) * mskCor(i+1,j+1,bi,bj)
567         &       +       maskC(i+1,j-1,k,bi,bj) * mskCor(i+1,j-1,bi,bj)
568         &       +       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,j,k,bi,bj)*mskCor(i,j,bi,bj)
571    #else
572             tmpVisc = GGL90visctmp(i,j,k)
573    #endif
574             tmpVisc = MIN(tmpVisc/TKEPrandtlNumber(i,j,k),GGL90diffMax)
575             GGL90diffKr(i,j,k,bi,bj)= MAX( tmpVisc , diffKrNrT(k) )
576            ENDDO
577           ENDDO
578          ENDDO
579    
580          DO k=2,Nr
581           DO j=1,sNy
582            DO i=1,sNx+1
583    #ifdef ALLOW_GGL90_SMOOTH
584            tmpVisc =
585         & (
586         &   p4 *(GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
587         &       +GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj))
588         &  +p8 *(GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj)
589         &       +GGL90visctmp(i-1,j+1,k) * mskCor(i-1,j+1,bi,bj)
590         &       +GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj)
591         &       +GGL90visctmp(i  ,j+1,k) * mskCor(i  ,j+1,bi,bj))
592         &  )
593         & /(p4 * 2. _d 0
594         &  +p8 *(      maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj)
595         &       +      maskC(i-1,j+1,k,bi,bj) * mskCor(i-1,j+1,bi,bj)
596         &       +      maskC(i  ,j-1,k,bi,bj) * mskCor(i  ,j-1,bi,bj)
597         &       +      maskC(i  ,j+1,k,bi,bj) * mskCor(i  ,j+1,bi,bj))
598         &  )
599         &  *maskC(i  ,j,k,bi,bj)*mskCor(i  ,j,bi,bj)
600         &  *maskC(i-1,j,k,bi,bj)*mskCor(i-1,j,bi,bj)
601    #else
602            tmpVisc = _maskW(i,j,k,bi,bj) *
603         &                   (.5 _d 0*(GGL90visctmp(i,j,k)
604         &                            +GGL90visctmp(i-1,j,k))
605         &                   )
606    #endif
607            tmpVisc = MIN( tmpVisc , GGL90viscMax )
608            GGL90viscArU(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
609            ENDDO
610           ENDDO
611          ENDDO
612    
613          DO k=2,Nr
614           DO j=1,sNy+1
615            DO i=1,sNx
616    #ifdef ALLOW_GGL90_SMOOTH
617            tmpVisc =
618         & (
619         &   p4 *(GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
620         &       +GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj))
621         &  +p8 *(GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj)
622         &       +GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj)
623         &       +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         &  )
626         & /(p4 * 2. _d 0
627         &  +p8 *(      maskC(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)
628         &       +      maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj)
629         &       +      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         &  )
632         &   *maskC(i,j  ,k,bi,bj)*mskCor(i,j  ,bi,bj)
633         &   *maskC(i,j-1,k,bi,bj)*mskCor(i,j-1,bi,bj)
634    #else
635            tmpVisc = _maskS(i,j,k,bi,bj) *
636         &                   (.5 _d 0*(GGL90visctmp(i,j,k)
637         &                            +GGL90visctmp(i,j-1,k))
638         &                   )
639    
640    #endif
641            tmpVisc = MIN( tmpVisc , GGL90viscMax )
642            GGL90viscArV(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
643            ENDDO
644           ENDDO
645          ENDDO
646    
647    #ifdef ALLOW_DIAGNOSTICS
648          IF ( useDiagnostics ) THEN
649             CALL DIAGNOSTICS_FILL( GGL90TKE   ,'GGL90TKE',
650         &                          0,Nr, 1, bi, bj, myThid )
651             CALL DIAGNOSTICS_FILL( GGL90viscArU,'GGL90ArU',
652         &                          0,Nr, 1, bi, bj, myThid )
653             CALL DIAGNOSTICS_FILL( GGL90viscArV,'GGL90ArV',
654         &                          0,Nr, 1, bi, bj, myThid )
655             CALL DIAGNOSTICS_FILL( GGL90diffKr,'GGL90Kr ',
656         &                          0,Nr, 1, bi, bj, myThid )
657             CALL DIAGNOSTICS_FILL( TKEPrandtlNumber ,'GGL90Prl',
658         &                          0,Nr, 2, bi, bj, myThid )
659             CALL DIAGNOSTICS_FILL( GGL90mixingLength,'GGL90Lmx',
660         &                          0,Nr, 2, bi, bj, myThid )
661          ENDIF
662    #endif
663    
664  #endif /* ALLOW_GGL90 */  #endif /* ALLOW_GGL90 */
665    
666        RETURN        RETURN
667        END        END
   
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.22
  ViewVC Help
Powered by ViewVC 1.1.22