/[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.14 by jmc, Sat Oct 31 03:18:50 2009 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     *==========================================================*
# Line 41  C !USES: =============================== Line 42  C !USES: ===============================
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     sigmaR :: Vertical gradient of iso-neutral density
47  C     myTime :: Current time in simulation  C     myTime :: Current time in simulation
48    C     myIter :: Current time-step number
49  C     myThid :: My Thread Id number  C     myThid :: My Thread Id number
50        INTEGER bi, bj        INTEGER bi, bj
51          _RL     sigmaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
52        _RL     myTime        _RL     myTime
53          INTEGER myIter
54        INTEGER myThid        INTEGER myThid
55  CEOP  CEOP
56    
# Line 53  CEOP Line 58  CEOP
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     TKEdissipation   - dissipation of TKE  C     KappaH           :: (local) diffusivity parameter for temperature (eq.11)
72  C     GGL90mixingLength- mixing length of scheme following Banke+Delecuse  C     KappaE           :: (local) diffusivity parameter for TKE (eq.15)
73  C         rMixingLength- inverse of mixing length  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
# Line 84  c     _RL     SQRTTKE Line 90  c     _RL     SQRTTKE
90        _RL     RiNumber        _RL     RiNumber
91        _RL     TKEdissipation        _RL     TKEdissipation
92        _RL     tempU, tempV, prTemp        _RL     tempU, tempV, prTemp
93        _RL     MaxLength, tmpmlx        _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)        _RL     rMixingLength    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
# Line 93  c     _RL     SQRTTKE Line 99  c     _RL     SQRTTKE
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  c     INTEGER errCode        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  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
119        _RL p4, p8, p16, tmpdiffKrS        _RL p4, p8, p16
120        p4=0.25 _d 0        p4=0.25 _d 0
121        p8=0.125 _d 0        p8=0.125 _d 0
122        p16=0.0625 _d 0        p16=0.0625 _d 0
# Line 122  C     set separate time step (should be Line 130  C     set separate time step (should be
130        deltaTggl90 = dTtracerLev(1)        deltaTggl90 = dTtracerLev(1)
131    
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)  = 1. _d 0           GGL90mixingLength(i,j,k) = GGL90mixingLengthMin
144           GGL90mixingLength(I,J,K) = GGL90mixingLengthMin           GGL90visctmp(i,j,k)      = 0. _d 0
145    #ifndef SOLVE_DIAGONAL_LOWMEMORY
146             a3d(i,j,k) = 0. _d 0
147             b3d(i,j,k) = 1. _d 0
148             c3d(i,j,k) = 0. _d 0
149    #endif
150          ENDDO          ENDDO
151         ENDDO         ENDDO
152        ENDDO        ENDDO
153        DO J=1-Oly,sNy+Oly        DO j=1-OLy,sNy+OLy
154         DO I=1-Olx,sNx+Olx         DO i=1-OLx,sNx+OLx
155          rhoK(I,J)          = 0. _d 0          rhoK(i,j)          = 0. _d 0
156          rhoKm1(I,J)        = 0. _d 0          rhoKm1(i,j)        = 0. _d 0
157          totalDepth(I,J)    = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)          totalDepth(i,j)    = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
158          rMixingLength(i,j,1) = 0. _d 0          rMixingLength(i,j,1) = 0. _d 0
159          mxLength_Dn(I,J,1) = GGL90mixingLengthMin          mxLength_Dn(i,j,1) = GGL90mixingLengthMin
160          SQRTTKE(i,j,1) = SQRT( GGL90TKE(i,j,1,bi,bj) )          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  c      Kp1 = MIN(Nr,K+1)  c      kp1 = MIN(Nr,k+1)
168         CALL FIND_RHO_2D(         CALL FIND_RHO_2D(
169       I      iMin, iMax, jMin, jMax, K,       I      iMin, iMax, jMin, jMax, k,
170       I      theta(1-OLx,1-OLy,Km1,bi,bj), salt(1-OLx,1-OLy,Km1,bi,bj),       I      theta(1-OLx,1-OLy,km1,bi,bj), salt(1-OLx,1-OLy,km1,bi,bj),
171       O      rhoKm1,       O      rhoKm1,
172       I      Km1, bi, bj, myThid )       I      km1, bi, bj, myThid )
173    
174         CALL FIND_RHO_2D(         CALL FIND_RHO_2D(
175       I      iMin, iMax, jMin, jMax, K,       I      iMin, iMax, jMin, jMax, k,
176       I      theta(1-OLx,1-OLy,K,bi,bj), salt(1-OLx,1-OLy,K,bi,bj),       I      theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
177       O      rhoK,       O      rhoK,
178       I      K, bi, bj, 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(i,j,k)=SQRT( GGL90TKE(I,J,K,bi,bj) )           SQRTTKE(i,j,k)=SQRT( GGL90TKE(i,j,k,bi,bj) )
182    
183  C     buoyancy frequency  C     buoyancy frequency
184           Nsquare(i,j,k) = - gravity*recip_rhoConst*recip_drC(K)           Nsquare(i,j,k) = gravity*gravitySign*recip_rhoConst
185       &        * ( rhoKm1(I,J) - rhoK(I,J) )*maskC(I,J,K,bi,bj)       &                  * sigmaR(i,j,k)
186             Nsquare(i,j,k) = - gravity*recip_rhoConst*recip_drC(k)
187         &        * ( rhoKm1(i,j) - rhoK(i,j) )*maskC(i,j,k,bi,bj)
188  cC     vertical shear term (dU/dz)^2+(dV/dz)^2  cC     vertical shear term (dU/dz)^2+(dV/dz)^2
189  c         tempU= .5 _d 0*( uVel(I,J,Km1,bi,bj)+uVel(I+1,J,Km1,bi,bj)  c         tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)
190  c     &                 -( uVel(I,J,K  ,bi,bj)+uVel(I+1,J,K  ,bi,bj)) )  c     &                 -( uVel(i,j,k  ,bi,bj)+uVel(i+1,j,k  ,bi,bj)) )
191  c     &        *recip_drC(K)  c     &        *recip_drC(k)
192  c         tempV= .5 _d 0*( vVel(I,J,Km1,bi,bj)+vVel(I,J+1,Km1,bi,bj)  c         tempV= .5 _d 0*( vVel(i,j,km1,bi,bj)+vVel(i,j+1,km1,bi,bj)
193  c     &                 -( vVel(I,J,K  ,bi,bj)+vVel(I,J+1,K  ,bi,bj)) )  c     &                 -( vVel(i,j,k  ,bi,bj)+vVel(i,j+1,k  ,bi,bj)) )
194  c     &        *recip_drC(K)  c     &        *recip_drC(k)
195  c         verticalShear = tempU*tempU + tempV*tempV  c         verticalShear = tempU*tempU + tempV*tempV
196  c         RiNumber   = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)  c         RiNumber   = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)
197  cC     compute Prandtl number (always greater than 0)  cC     compute Prandtl number (always greater than 0)
198  c         prTemp = 1. _d 0  c         prTemp = 1. _d 0
199  c         IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber  c         IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
200  c         TKEPrandtlNumber(I,J,K) = MIN(10. _d 0,prTemp)  c         TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
201  C     mixing length  C     mixing length
202           GGL90mixingLength(I,J,K) = SQRTTWO *           GGL90mixingLength(i,j,k) = SQRTTWO *
203       &        SQRTTKE(i,j,k)/SQRT( MAX(Nsquare(i,j,k),GGL90eps) )       &        SQRTTKE(i,j,k)/SQRT( MAX(Nsquare(i,j,k),GGL90eps) )
204          ENDDO          ENDDO
205         ENDDO         ENDDO
# Line 194  C     mixing length Line 207  C     mixing length
207    
208  C- Impose upper and lower bound for mixing length  C- Impose upper and lower bound for mixing length
209        IF ( mxlMaxFlag .EQ. 0 ) THEN        IF ( mxlMaxFlag .EQ. 0 ) THEN
210  C-  
211         DO k=2,Nr         DO k=2,Nr
212          DO J=jMin,jMax          DO j=jMin,jMax
213           DO I=iMin,iMax           DO i=iMin,iMax
214            MaxLength=totalDepth(I,J)            MaxLength=totalDepth(i,j)
215            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
216       &                                   MaxLength)       &                                   MaxLength)
217           ENDDO           ENDDO
218          ENDDO          ENDDO
219         ENDDO         ENDDO
220    
221         DO k=2,Nr         DO k=2,Nr
222          DO J=jMin,jMax          DO j=jMin,jMax
223           DO I=iMin,iMax           DO i=iMin,iMax
224            GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
225       &                                   GGL90mixingLengthMin)       &                                   GGL90mixingLengthMin)
226            rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K)            rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
227           ENDDO           ENDDO
228          ENDDO          ENDDO
229         ENDDO         ENDDO
230    
231        ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN
232  C-  
233         DO k=2,Nr         DO k=2,Nr
234          DO J=jMin,jMax          DO j=jMin,jMax
235           DO I=iMin,iMax           DO i=iMin,iMax
236            MaxLength=MIN(Ro_surf(I,J,bi,bj)-rF(k),rF(k)-R_low(I,J,bi,bj))            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)  c         MaxLength=MAX(MaxLength,20. _d 0)
238            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
239       &                                   MaxLength)       &                                   MaxLength)
240           ENDDO           ENDDO
241          ENDDO          ENDDO
242         ENDDO         ENDDO
243    
244         DO k=2,Nr         DO k=2,Nr
245          DO J=jMin,jMax          DO j=jMin,jMax
246           DO I=iMin,iMax           DO i=iMin,iMax
247            GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
248       &                                   GGL90mixingLengthMin)       &                                   GGL90mixingLengthMin)
249            rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K)            rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
250           ENDDO           ENDDO
251          ENDDO          ENDDO
252         ENDDO         ENDDO
253    
254        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN
255  C-  
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         DO k=2,Nr
264          DO J=jMin,jMax          DO j=jMin,jMax
265           DO I=iMin,iMax           DO i=iMin,iMax
266            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
267       &        GGL90mixingLength(I,J,K-1)+drF(k-1))       &        GGL90mixingLength(i,j,k-1)+drF(k-1))
268           ENDDO           ENDDO
269          ENDDO          ENDDO
270         ENDDO         ENDDO
271         DO J=jMin,jMax         DO j=jMin,jMax
272          DO I=iMin,iMax          DO i=iMin,iMax
273            GGL90mixingLength(I,J,Nr) = MIN(GGL90mixingLength(I,J,Nr),            GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
274       &       GGL90mixingLengthMin+drF(Nr))       &       GGL90mixingLengthMin+drF(Nr))
275          ENDDO          ENDDO
276         ENDDO         ENDDO
277         DO k=Nr-1,2,-1         DO k=Nr-1,2,-1
278          DO J=jMin,jMax          DO j=jMin,jMax
279           DO I=iMin,iMax           DO i=iMin,iMax
280            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
281       &        GGL90mixingLength(I,J,K+1)+drF(k))       &        GGL90mixingLength(i,j,k+1)+drF(k))
282           ENDDO           ENDDO
283          ENDDO          ENDDO
284         ENDDO         ENDDO
285    
286         DO k=2,Nr         DO k=2,Nr
287          DO J=jMin,jMax          DO j=jMin,jMax
288           DO I=iMin,iMax           DO i=iMin,iMax
289            GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
290       &                                   GGL90mixingLengthMin)       &                                   GGL90mixingLengthMin)
291            rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K)            rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
292           ENDDO           ENDDO
293          ENDDO          ENDDO
294         ENDDO         ENDDO
295    
296        ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN
297  C-  
298         DO k=2,Nr         DO k=2,Nr
299          DO J=jMin,jMax          DO j=jMin,jMax
300           DO I=iMin,iMax           DO i=iMin,iMax
301            mxLength_Dn(I,J,K) = MIN(GGL90mixingLength(I,J,K),            mxLength_Dn(i,j,k) = MIN(GGL90mixingLength(i,j,k),
302       &        mxLength_Dn(I,J,K-1)+drF(k-1))       &        mxLength_Dn(i,j,k-1)+drF(k-1))
303           ENDDO           ENDDO
304          ENDDO          ENDDO
305         ENDDO         ENDDO
306         DO J=jMin,jMax         DO j=jMin,jMax
307          DO I=iMin,iMax          DO i=iMin,iMax
308            GGL90mixingLength(I,J,Nr) = MIN(GGL90mixingLength(I,J,Nr),            GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
309       &       GGL90mixingLengthMin+drF(Nr))       &       GGL90mixingLengthMin+drF(Nr))
310          ENDDO          ENDDO
311         ENDDO         ENDDO
312         DO k=Nr-1,2,-1         DO k=Nr-1,2,-1
313          DO J=jMin,jMax          DO j=jMin,jMax
314           DO I=iMin,iMax           DO i=iMin,iMax
315            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
316       &        GGL90mixingLength(I,J,K+1)+drF(k))       &        GGL90mixingLength(i,j,k+1)+drF(k))
317           ENDDO           ENDDO
318          ENDDO          ENDDO
319         ENDDO         ENDDO
320    
321         DO k=2,Nr         DO k=2,Nr
322          DO J=jMin,jMax          DO j=jMin,jMax
323           DO I=iMin,iMax           DO i=iMin,iMax
324            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
325       &                                  mxLength_Dn(I,J,K))       &                                  mxLength_Dn(i,j,k))
326            tmpmlx = SQRT( GGL90mixingLength(I,J,K)*mxLength_Dn(I,J,K) )            tmpmlx = SQRT( GGL90mixingLength(i,j,k)*mxLength_Dn(i,j,k) )
327            tmpmlx = MAX( tmpmlx, GGL90mixingLengthMin)            tmpmlx = MAX( tmpmlx, GGL90mixingLengthMin)
328            rMixingLength(I,J,K) = 1. _d 0 / tmpmlx            rMixingLength(i,j,k) = 1. _d 0 / tmpmlx
329           ENDDO           ENDDO
330          ENDDO          ENDDO
331         ENDDO         ENDDO
# Line 316  C- Line 336  C-
336    
337  C- Impose minimum mixing length (to avoid division by zero)  C- Impose minimum mixing length (to avoid division by zero)
338  c      DO k=2,Nr  c      DO k=2,Nr
339  c      DO J=jMin,jMax  c      DO j=jMin,jMax
340  c       DO I=iMin,iMax  c       DO i=iMin,iMax
341  c        GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),  c        GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
342  c    &        GGL90mixingLengthMin)  c    &        GGL90mixingLengthMin)
343  c        rMixingLength(I,J,K) = 1. _d 0 /GGL90mixingLength(I,J,K)  c        rMixingLength(i,j,k) = 1. _d 0 /GGL90mixingLength(i,j,k)
344  c       ENDDO  c       ENDDO
345  c      ENDDO  c      ENDDO
346  c     ENDDO  c     ENDDO
347    
348        DO k=2,Nr        DO k=2,Nr
349         Km1 = K-1         km1 = k-1
        DO J=jMin,jMax  
         DO I=iMin,iMax  
 C     vertical shear term (dU/dz)^2+(dV/dz)^2  
          tempU= .5 _d 0*( uVel(I,J,Km1,bi,bj)+uVel(I+1,J,Km1,bi,bj)  
      &                 -( uVel(I,J,K  ,bi,bj)+uVel(I+1,J,K  ,bi,bj)) )  
      &        *recip_drC(K)  
          tempV= .5 _d 0*( vVel(I,J,Km1,bi,bj)+vVel(I,J+1,Km1,bi,bj)  
      &                 -( vVel(I,J,K  ,bi,bj)+vVel(I,J+1,K  ,bi,bj)) )  
      &        *recip_drC(K)  
          verticalShear = tempU*tempU + tempV*tempV  
          RiNumber = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)  
 C     compute Prandtl number (always greater than 0)  
          prTemp = 1. _d 0  
          IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber  
          TKEPrandtlNumber(I,J,K) = MIN(10. _d 0,prTemp)  
 c         TKEPrandtlNumber(I,J,K) = 1. _d 0  
   
 C     viscosity and diffusivity  
          KappaM = GGL90ck*GGL90mixingLength(I,J,K)*SQRTTKE(i,j,k)  
          KappaH = KappaM/TKEPrandtlNumber(I,J,K)  
   
 C     Set a minium (= background) and maximum value  
          KappaM = MAX(KappaM,viscArNr(k))  
          KappaH = MAX(KappaH,diffKrNrT(k))  
          KappaM = MIN(KappaM,GGL90viscMax)  
          KappaH = MIN(KappaH,GGL90diffMax)  
   
 C     Mask land points and storage  
          GGL90viscAr(I,J,K,bi,bj) = KappaM * maskC(I,J,K,bi,bj)  
          GGL90diffKr(I,J,K,bi,bj) = KappaH * maskC(I,J,K,bi,bj)  
          KappaE(I,J,K) = GGL90alpha * GGL90viscAr(I,J,K,bi,bj)  
   
 C     dissipation term  
          TKEdissipation = ab05*GGL90ceps  
      &        *SQRTTKE(i,j,k)*rMixingLength(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  
      &        - KappaH*Nsquare(i,j,k)  
      &        - TKEdissipation  
      &        )  
         ENDDO  
        ENDDO  
       ENDDO  
350    
351    #ifdef ALLOW_GGL90_HORIZDIFF
352           IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
353  C     horizontal diffusion of TKE (requires an exchange in  C     horizontal diffusion of TKE (requires an exchange in
354  C      do_fields_blocking_exchanges)  C      do_fields_blocking_exchanges)
 #ifdef ALLOW_GGL90_HORIZDIFF  
       IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN  
        DO K=2,Nr  
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)
# Line 390  C     common factors Line 363  C     common factors
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. _d 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 400  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. _d 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 414  C     ... across y-faces Line 387  C     ... across y-faces
387           ENDDO           ENDDO
388          ENDDO          ENDDO
389  C     Compute divergence of fluxes  C     Compute divergence of fluxes
390  C- jmc: concerned about missing a deltaT since gTKE is already the future TKE.          DO j=1-OLy,sNy+OLy-1
391          DO j=1-Oly,sNy+Oly-1           DO i=1-OLx,sNx+OLx-1
392           DO i=1-Olx,sNx+Olx-1            gTKE(i,j) =
393            gTKE(i,j,k)=gTKE(i,j,k)       &    -_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          ENDDO
399  C       end of k-loop  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
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 438  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  C- jmc: concerned that a(k=2) should always be zero         km1=MAX(2,k-1)
 C       and would be better without recip_hFacC  
        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  c     &        *recip_drF(km1)*recip_hFacI(i,j,k,bi,bj)  C-    but no hFacC in TKE volume control
476    C-    No need for maskC(k-1) with recip_hFacC(k-1)
477             a3d(i,j,k) = -deltaTggl90
478       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)
479       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))
480       &        *recip_drC(k)*maskC(i,j,k,bi,bj)*maskC(i,j,k-1,bi,bj)       &        *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
488         ENDDO         ENDDO
489        ENDDO        ENDDO
490        DO k=2,Nr        DO k=2,Nr
 C- jmc: concerned that c(k) from k=kLow to k=Nr should always be zero  
491         DO j=jMin,jMax         DO j=jMin,jMax
492          DO i=iMin,iMax          DO i=iMin,iMax
493  C- jmc: this is dangerous since klowC could be zero if land column            kp1=MAX(1,MIN(klowC(i,j,bi,bj),k+1))
494  C       and would be better without recip_hFacC  C-    We keep recip_hFacC in the diffusive flux calculation,
495            kp1=min(klowC(i,j,bi,bj),k+1)  C-    but no hFacC in TKE volume control
496            c(i,j,k) = -deltaTggl90  C-    No need for maskC(k) with recip_hFacC(k)
497  c     &        *recip_drF( k )*recip_hFacI(i,j,k,bi,bj)            c3d(i,j,k) = -deltaTggl90
498       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)
499       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))
500       &        *recip_drC(k)*maskC(i,j,k,bi,bj)*maskC(i,j,k-1,bi,bj)       &        *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*SQRTTKE(I,J,K)       &        + implDissFac*deltaTggl90*GGL90ceps*SQRTTKE(i,j,k)
511       &        *rMixingLength(I,J,K)*maskC(i,j,k,bi,bj)       &        * 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    
518  C     Apply boundary condition  C     Apply boundary condition
519  C- jmc: concerned about conservation when a or c are changed after computing b        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 _d 0*( 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 _d 0*( 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  C- jmc: would be much better to update the provisional TKE (i.e. gTKE) at k=2          GGL90TKE(i,j,kSurf,bi,bj) = maskC(i,j,kSurf,bi,bj)
531          gTKE(I,J,kSurf) = MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)       &           *MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)
532       &                     *maskC(I,J,kSurf,bi,bj)          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   = MAX(kLowC(I,J,bi,bj),1)          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          c(I,J,kBottom) = 0. _d 0          c3d(i,j,kBottom) = 0. _d 0
540         ENDDO         ENDDO
541        ENDDO        ENDDO
542    
543  C     solve tri-diagonal system, and store solution on gTKE (previously rhs)  C     solve tri-diagonal system
544  c     CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
545  c    I                        a, b, c,       I                        a3d, b3d, c3d,
546  c    U                        gTKE,       U                        GGL90TKE,
547  c    O                        errCode,       O                        errCode,
548  c    I                        bi, bj, myThid )       I                        bi, bj, myThid )
549        CALL GGL90_SOLVE( bi, bj, iMin, iMax, jMin, jMax,  
550       I     a, b, c,        DO k=1,Nr
551       U     gTKE,         DO j=jMin,jMax
552       I     myThid )          DO i=iMin,iMax
   
 C     now update TKE  
       DO K=1,Nr  
        DO J=jMin,jMax  
         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
# Line 536  C     impose minimum TKE to avoid numeri Line 560  C     impose minimum TKE to avoid numeri
560  C     end of time step  C     end of time step
561  C     ===============================  C     ===============================
562    
563          DO k=2,Nr
564           DO j=1,sNy
565            DO i=1,sNx
566  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
567        DO K=1,Nr           tmpVisc=
        DO J=jMin,jMax  
         DO I=iMin,iMax  
          tmpdiffKrS=  
568       &  (       &  (
569       &   p4 *  GGL90viscAr(i  ,j  ,k,bi,bj) * mskCor(i  ,j  ,bi,bj)       &   p4 *  GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
570       &  +p8 *( GGL90viscAr(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)       &  +p8 *( GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj)
571       &       + GGL90viscAr(i  ,j-1,k,bi,bj) * mskCor(i  ,j-1,bi,bj)       &       + GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj)
572       &       + GGL90viscAr(i+1,j  ,k,bi,bj) * mskCor(i+1,j  ,bi,bj)       &       + GGL90visctmp(i+1,j  ,k) * mskCor(i+1,j  ,bi,bj)
573       &       + GGL90viscAr(i  ,j+1,k,bi,bj) * mskCor(i  ,j+1,bi,bj))       &       + GGL90visctmp(i  ,j+1,k) * mskCor(i  ,j+1,bi,bj))
574       &  +p16*( GGL90viscAr(i+1,j+1,k,bi,bj) * mskCor(i+1,j+1,bi,bj)       &  +p16*( GGL90visctmp(i+1,j+1,k) * mskCor(i+1,j+1,bi,bj)
575       &       + GGL90viscAr(i+1,j-1,k,bi,bj) * mskCor(i+1,j-1,bi,bj)       &       + GGL90visctmp(i+1,j-1,k) * mskCor(i+1,j-1,bi,bj)
576       &       + GGL90viscAr(i-1,j+1,k,bi,bj) * mskCor(i-1,j+1,bi,bj)       &       + GGL90visctmp(i-1,j+1,k) * mskCor(i-1,j+1,bi,bj)
577       &       + GGL90viscAr(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj))       &       + GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj))
578       &  )       &  )
579       & /(p4       & /(p4
580       &  +p8 *(       maskC(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)       &  +p8 *(       maskC(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)
# Line 562  C     =============================== Line 586  C     ===============================
586       &       +       maskC(i-1,j+1,k,bi,bj) * mskCor(i-1,j+1,bi,bj)       &       +       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))       &       +       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)       &  )*maskC(i,j,k,bi,bj)*mskCor(i,j,bi,bj)
589       &   /TKEPrandtlNumber(i,j,k)  #else
590           GGL90diffKrS(I,J,K,bi,bj)= MAX( tmpdiffKrS , diffKrNrT(k) )           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          ENDDO
628         ENDDO         ENDDO
629        ENDDO        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  #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  #ifdef ALLOW_DIAGNOSTICS
666        IF ( useDiagnostics ) THEN        IF ( useDiagnostics ) THEN
667           CALL DIAGNOSTICS_FILL( GGL90TKE   ,'GGL90TKE',           CALL DIAGNOSTICS_FILL( GGL90TKE   ,'GGL90TKE',
668       &                          0,Nr, 1, bi, bj, myThid )       &                          0,Nr, 1, bi, bj, myThid )
669           CALL DIAGNOSTICS_FILL( GGL90viscAr,'GGL90Ar ',           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 )       &                          0,Nr, 1, bi, bj, myThid )
673           CALL DIAGNOSTICS_FILL( GGL90diffKr,'GGL90Kr ',           CALL DIAGNOSTICS_FILL( GGL90diffKr,'GGL90Kr ',
674       &                          0,Nr, 1, bi, bj, myThid )       &                          0,Nr, 1, bi, bj, myThid )
Legend:
Removed from v.1.14  
changed lines
  Added in v.1.21
  ViewVC Help
Powered by ViewVC 1.1.22