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