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