/[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.15 by dfer, Sat Nov 14 14:27:56 2009 UTC revision 1.24 by mlosch, Mon Nov 5 08:11:45 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        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
# Line 84  c     _RL     SQRTTKE Line 89  c     _RL     SQRTTKE
89        _RL     RiNumber        _RL     RiNumber
90        _RL     TKEdissipation        _RL     TKEdissipation
91        _RL     tempU, tempV, prTemp        _RL     tempU, tempV, prTemp
92        _RL     MaxLength, tmpmlx        _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)        _RL     rMixingLength    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
96        _RL     mxLength_Dn      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _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        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  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
118        _RL p4, p8, p16, tmpdiffKrS        _RL p4, p8, p16
119        p4=0.25 _d 0        p4=0.25 _d 0
120        p8=0.125 _d 0        p8=0.125 _d 0
121        p16=0.0625 _d 0        p16=0.0625 _d 0
# Line 122  C     set separate time step (should be Line 129  C     set separate time step (should be
129        deltaTggl90 = dTtracerLev(1)        deltaTggl90 = dTtracerLev(1)
130    
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)  = 1. _d 0           GGL90mixingLength(i,j,k) = GGL90mixingLengthMin
143           GGL90mixingLength(I,J,K) = GGL90mixingLengthMin           GGL90visctmp(i,j,k)      = 0. _d 0
144    #ifndef SOLVE_DIAGONAL_LOWMEMORY
145             a3d(i,j,k) = 0. _d 0
146             b3d(i,j,k) = 1. _d 0
147             c3d(i,j,k) = 0. _d 0
148    #endif
149          ENDDO          ENDDO
150         ENDDO         ENDDO
151        ENDDO        ENDDO
152        DO J=1-Oly,sNy+Oly        DO j=1-OLy,sNy+OLy
153         DO I=1-Olx,sNx+Olx         DO i=1-OLx,sNx+OLx
154          rhoK(I,J)          = 0. _d 0          totalDepth(i,j)    = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
         rhoKm1(I,J)        = 0. _d 0  
         totalDepth(I,J)    = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)  
155          rMixingLength(i,j,1) = 0. _d 0          rMixingLength(i,j,1) = 0. _d 0
156          mxLength_Dn(I,J,1) = GGL90mixingLengthMin          mxLength_Dn(i,j,1) = GGL90mixingLengthMin
157          SQRTTKE(i,j,1) = SQRT( GGL90TKE(i,j,1,bi,bj) )          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  c      Kp1 = MIN(Nr,K+1)  c      kp1 = MIN(Nr,k+1)
165         CALL FIND_RHO_2D(         DO j=jMin,jMax
166       I      iMin, iMax, jMin, jMax, K,          DO i=iMin,iMax
167       I      theta(1-OLx,1-OLy,Km1,bi,bj), salt(1-OLx,1-OLy,Km1,bi,bj),           SQRTTKE(i,j,k)=SQRT( GGL90TKE(i,j,k,bi,bj) )
      O      rhoKm1,  
      I      Km1, bi, bj, myThid )  
   
        CALL FIND_RHO_2D(  
      I      iMin, iMax, jMin, jMax, K,  
      I      theta(1-OLx,1-OLy,K,bi,bj), salt(1-OLx,1-OLy,K,bi,bj),  
      O      rhoK,  
      I      K, bi, bj, myThid )  
        DO J=jMin,jMax  
         DO I=iMin,iMax  
          SQRTTKE(i,j,k)=SQRT( GGL90TKE(I,J,K,bi,bj) )  
168    
169  C     buoyancy frequency  C     buoyancy frequency
170           Nsquare(i,j,k) = - gravity*recip_rhoConst*recip_drC(K)           Nsquare(i,j,k) = gravity*gravitySign*recip_rhoConst
171       &        * ( rhoKm1(I,J) - rhoK(I,J) )*maskC(I,J,K,bi,bj)       &                  * sigmaR(i,j,k)
172  cC     vertical shear term (dU/dz)^2+(dV/dz)^2  cC     vertical shear term (dU/dz)^2+(dV/dz)^2
173  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)
174  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)) )
175  c     &        *recip_drC(K)  c     &        *recip_drC(k)
176  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)
177  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)) )
178  c     &        *recip_drC(K)  c     &        *recip_drC(k)
179  c         verticalShear = tempU*tempU + tempV*tempV  c         verticalShear = tempU*tempU + tempV*tempV
180  c         RiNumber   = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)  c         RiNumber   = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)
181  cC     compute Prandtl number (always greater than 0)  cC     compute Prandtl number (always greater than 0)
182  c         prTemp = 1. _d 0  c         prTemp = 1. _d 0
183  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
184  c         TKEPrandtlNumber(I,J,K) = MIN(10. _d 0,prTemp)  c         TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
185  C     mixing length  C     mixing length
186           GGL90mixingLength(I,J,K) = SQRTTWO *           GGL90mixingLength(i,j,k) = SQRTTWO *
187       &        SQRTTKE(i,j,k)/SQRT( MAX(Nsquare(i,j,k),GGL90eps) )       &        SQRTTKE(i,j,k)/SQRT( MAX(Nsquare(i,j,k),GGL90eps) )
188          ENDDO          ENDDO
189         ENDDO         ENDDO
190        ENDDO        ENDDO
191    
192    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  C- Impose upper and lower bound for mixing length
202        IF ( mxlMaxFlag .EQ. 0 ) THEN        IF ( mxlMaxFlag .EQ. 0 ) THEN
203  C-  
204         DO k=2,Nr         DO k=2,Nr
205          DO J=jMin,jMax          DO j=jMin,jMax
206           DO I=iMin,iMax           DO i=iMin,iMax
207            MaxLength=totalDepth(I,J)            MaxLength=totalDepth(i,j)
208            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
209       &                                   MaxLength)       &                                   MaxLength)
210           ENDDO           ENDDO
211          ENDDO          ENDDO
212         ENDDO         ENDDO
213    
214         DO k=2,Nr         DO k=2,Nr
215          DO J=jMin,jMax          DO j=jMin,jMax
216           DO I=iMin,iMax           DO i=iMin,iMax
217            GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
218       &                                   GGL90mixingLengthMin)       &                                   GGL90mixingLengthMin)
219            rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K)            rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
220           ENDDO           ENDDO
221          ENDDO          ENDDO
222         ENDDO         ENDDO
223    
224        ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN
225  C-  
226         DO k=2,Nr         DO k=2,Nr
227          DO J=jMin,jMax          DO j=jMin,jMax
228           DO I=iMin,iMax           DO i=iMin,iMax
229            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))
230  c         MaxLength=MAX(MaxLength,20. _d 0)  c         MaxLength=MAX(MaxLength,20. _d 0)
231            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
232       &                                   MaxLength)       &                                   MaxLength)
233           ENDDO           ENDDO
234          ENDDO          ENDDO
235         ENDDO         ENDDO
236    
237         DO k=2,Nr         DO k=2,Nr
238          DO J=jMin,jMax          DO j=jMin,jMax
239           DO I=iMin,iMax           DO i=iMin,iMax
240            GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
241       &                                   GGL90mixingLengthMin)       &                                   GGL90mixingLengthMin)
242            rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K)            rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
243           ENDDO           ENDDO
244          ENDDO          ENDDO
245         ENDDO         ENDDO
246    
247        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN
248  C-  
249         DO k=2,Nr         DO k=2,Nr
250          DO J=jMin,jMax          DO j=jMin,jMax
251           DO I=iMin,iMax           DO i=iMin,iMax
252            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
253       &        GGL90mixingLength(I,J,K-1)+drF(k-1))       &        GGL90mixingLength(i,j,k-1)+drF(k-1))
254           ENDDO           ENDDO
255          ENDDO          ENDDO
256         ENDDO         ENDDO
257         DO J=jMin,jMax         DO j=jMin,jMax
258          DO I=iMin,iMax          DO i=iMin,iMax
259            GGL90mixingLength(I,J,Nr) = MIN(GGL90mixingLength(I,J,Nr),            GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
260       &       GGL90mixingLengthMin+drF(Nr))       &       GGL90mixingLengthMin+drF(Nr))
261          ENDDO          ENDDO
262         ENDDO         ENDDO
263         DO k=Nr-1,2,-1         DO k=Nr-1,2,-1
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))       &        GGL90mixingLength(i,j,k+1)+drF(k))
268           ENDDO           ENDDO
269          ENDDO          ENDDO
270         ENDDO         ENDDO
271    
272         DO k=2,Nr         DO k=2,Nr
273          DO J=jMin,jMax          DO j=jMin,jMax
274           DO I=iMin,iMax           DO i=iMin,iMax
275            GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
276       &                                   GGL90mixingLengthMin)       &                                   GGL90mixingLengthMin)
277            rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K)            rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
278           ENDDO           ENDDO
279          ENDDO          ENDDO
280         ENDDO         ENDDO
281    
282        ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN
283  C-  
284         DO k=2,Nr         DO k=2,Nr
285          DO J=jMin,jMax          DO j=jMin,jMax
286           DO I=iMin,iMax           DO i=iMin,iMax
287            mxLength_Dn(I,J,K) = MIN(GGL90mixingLength(I,J,K),            mxLength_Dn(i,j,k) = MIN(GGL90mixingLength(i,j,k),
288       &        mxLength_Dn(I,J,K-1)+drF(k-1))       &        mxLength_Dn(i,j,k-1)+drF(k-1))
289           ENDDO           ENDDO
290          ENDDO          ENDDO
291         ENDDO         ENDDO
292         DO J=jMin,jMax         DO j=jMin,jMax
293          DO I=iMin,iMax          DO i=iMin,iMax
294            GGL90mixingLength(I,J,Nr) = MIN(GGL90mixingLength(I,J,Nr),            GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
295       &       GGL90mixingLengthMin+drF(Nr))       &       GGL90mixingLengthMin+drF(Nr))
296          ENDDO          ENDDO
297         ENDDO         ENDDO
298         DO k=Nr-1,2,-1         DO k=Nr-1,2,-1
299          DO J=jMin,jMax          DO j=jMin,jMax
300           DO I=iMin,iMax           DO i=iMin,iMax
301            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
302       &        GGL90mixingLength(I,J,K+1)+drF(k))       &        GGL90mixingLength(i,j,k+1)+drF(k))
303           ENDDO           ENDDO
304          ENDDO          ENDDO
305         ENDDO         ENDDO
306    
307         DO k=2,Nr         DO k=2,Nr
308          DO J=jMin,jMax          DO j=jMin,jMax
309           DO I=iMin,iMax           DO i=iMin,iMax
310            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
311       &                                  mxLength_Dn(I,J,K))       &                                  mxLength_Dn(i,j,k))
312            tmpmlx = SQRT( GGL90mixingLength(I,J,K)*mxLength_Dn(I,J,K) )            tmpmlx = SQRT( GGL90mixingLength(i,j,k)*mxLength_Dn(i,j,k) )
313            tmpmlx = MAX( tmpmlx, GGL90mixingLengthMin)            tmpmlx = MAX( tmpmlx, GGL90mixingLengthMin)
314            rMixingLength(I,J,K) = 1. _d 0 / tmpmlx            rMixingLength(i,j,k) = 1. _d 0 / tmpmlx
315           ENDDO           ENDDO
316          ENDDO          ENDDO
317         ENDDO         ENDDO
# Line 316  C- Line 322  C-
322    
323  C- Impose minimum mixing length (to avoid division by zero)  C- Impose minimum mixing length (to avoid division by zero)
324  c      DO k=2,Nr  c      DO k=2,Nr
325  c      DO J=jMin,jMax  c      DO j=jMin,jMax
326  c       DO I=iMin,iMax  c       DO i=iMin,iMax
327  c        GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),  c        GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
328  c    &        GGL90mixingLengthMin)  c    &        GGL90mixingLengthMin)
329  c        rMixingLength(I,J,K) = 1. _d 0 /GGL90mixingLength(I,J,K)  c        rMixingLength(i,j,k) = 1. _d 0 /GGL90mixingLength(i,j,k)
330  c       ENDDO  c       ENDDO
331  c      ENDDO  c      ENDDO
332  c     ENDDO  c     ENDDO
333    
334        DO k=2,Nr        DO k=2,Nr
335         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  
336    
 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  
   
 C     horizontal diffusion of TKE (requires an exchange in  
 C      do_fields_blocking_exchanges)  
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. _d 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         &         *maskW(i,j,k,bi,bj) ! paranoia setting
364           ENDDO           ENDDO
365          ENDDO          ENDDO
366  C     ... across y-faces  C     ... across y-faces
367          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
368           dfy(i,1-Oly)=0. _d 0           dfy(i,1-OLy)=0. _d 0
369          ENDDO          ENDDO
370          DO j=1-Oly+1,sNy+Oly          DO j=1-OLy+1,sNy+OLy
371           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
372            dfy(i,j) = -GGL90diffTKEh*yA(i,j)            dfy(i,j) = -GGL90diffTKEh*yA(i,j)
373       &      *_recip_dyC(i,j,bi,bj)       &      *_recip_dyC(i,j,bi,bj)
374       &      *(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))
375  #ifdef ISOTROPIC_COS_SCALING  #ifdef ISOTROPIC_COS_SCALING
376       &      *CosFacV(j,bi,bj)       &      *CosFacV(j,bi,bj)
377  #endif /* ISOTROPIC_COS_SCALING */  #endif /* ISOTROPIC_COS_SCALING */
378         &         *maskS(i,j,k,bi,bj) ! paranoia setting
379           ENDDO           ENDDO
380          ENDDO          ENDDO
381  C     Compute divergence of fluxes  C     Compute divergence of fluxes
382          DO j=1-Oly,sNy+Oly-1          DO j=1-OLy,sNy+OLy-1
383           DO i=1-Olx,sNx+Olx-1           DO i=1-OLx,sNx+OLx-1
384            gTKE(i,j,k)=gTKE(i,j,k)            hFac = min(.5 _d 0,_hFacC(i,j,k-1,bi,bj) ) +
385       &   -_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) )
386       &         *( (dfx(i+1,j)-dfx(i,j))            gTKE(i,j) = 0.0
387       &           +(dfy(i,j+1)-dfy(i,j))            if ( hFac .ne. 0.0 )
388       &           )*deltaTggl90       &      gTKE(i,j) = -recip_drC(k)*recip_rA(i,j,bi,bj)/hFac
389         &         *((dfx(i+1,j)-dfx(i,j))
390         &          +(dfy(i,j+1)-dfy(i,j)) )
391         &         *maskC(i,j,k,bi,bj) ! paranoia setting
392           ENDDO           ENDDO
393          ENDDO          ENDDO
394  C       end of k-loop  C      end if GGL90diffTKEh .eq. 0.
395           ENDIF
396    #endif /* ALLOW_GGL90_HORIZDIFF */
397    
398           DO j=jMin,jMax
399            DO i=iMin,iMax
400    C     vertical shear term (dU/dz)^2+(dV/dz)^2
401             tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)
402         &                 -( uVel(i,j,k  ,bi,bj)+uVel(i+1,j,k  ,bi,bj)) )
403         &        *recip_drC(k)
404             tempV= .5 _d 0*( vVel(i,j,km1,bi,bj)+vVel(i,j+1,km1,bi,bj)
405         &                 -( vVel(i,j,k  ,bi,bj)+vVel(i,j+1,k  ,bi,bj)) )
406         &        *recip_drC(k)
407             verticalShear = tempU*tempU + tempV*tempV
408             RiNumber = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)
409    C     compute Prandtl number (always greater than 0)
410             prTemp = 1. _d 0
411             IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
412             TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
413    c         TKEPrandtlNumber(i,j,k) = 1. _d 0
414    
415    C     viscosity and diffusivity
416             KappaM = GGL90ck*GGL90mixingLength(i,j,k)*SQRTTKE(i,j,k)
417             GGL90visctmp(i,j,k) = MAX(KappaM,diffKrNrT(k))
418         &                            * maskC(i,j,k,bi,bj)
419    c        note: storing GGL90visctmp like this, and using it later to compute
420    c              GGL9rdiffKr etc. is robust in case of smoothing (e.g. see OPA)
421             KappaM = MAX(KappaM,viscArNr(k)) * maskC(i,j,k,bi,bj)
422             KappaH = KappaM/TKEPrandtlNumber(i,j,k)
423             KappaE(i,j,k) = GGL90alpha * KappaM * maskC(i,j,k,bi,bj)
424    
425    C     dissipation term
426             TKEdissipation = explDissFac*GGL90ceps
427         &        *SQRTTKE(i,j,k)*rMixingLength(i,j,k)
428         &        *GGL90TKE(i,j,k,bi,bj)
429    C     partial update with sum of explicit contributions
430             GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
431         &        + deltaTggl90*(
432         &        + KappaM*verticalShear
433         &        - KappaH*Nsquare(i,j,k)
434         &        - TKEdissipation
435         &        )
436            ENDDO
437         ENDDO         ENDDO
438  C     end if GGL90diffTKEh .eq. 0.  
439        ENDIF  #ifdef ALLOW_GGL90_HORIZDIFF
440           IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
441    C--    Add horiz. diffusion tendency
442            DO j=jMin,jMax
443             DO i=iMin,iMax
444              GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
445         &                          + gTKE(i,j)*deltaTggl90
446             ENDDO
447            ENDDO
448           ENDIF
449  #endif /* ALLOW_GGL90_HORIZDIFF */  #endif /* ALLOW_GGL90_HORIZDIFF */
450    
451    C--   end of k loop
452          ENDDO
453    
454  C     ============================================  C     ============================================
455  C     Implicit time step to update TKE for k=1,Nr;  C     Implicit time step to update TKE for k=1,Nr;
456  C     TKE(Nr+1)=0 by default  C     TKE(Nr+1)=0 by default
# Line 437  C     set up matrix Line 459  C     set up matrix
459  C--   Lower diagonal  C--   Lower diagonal
460        DO j=jMin,jMax        DO j=jMin,jMax
461         DO i=iMin,iMax         DO i=iMin,iMax
462           a(i,j,1) = 0. _d 0           a3d(i,j,1) = 0. _d 0
463         ENDDO         ENDDO
464        ENDDO        ENDDO
465        DO k=2,Nr        DO k=2,Nr
# Line 447  C--   Lower diagonal Line 469  C--   Lower diagonal
469  C-    We keep recip_hFacC in the diffusive flux calculation,  C-    We keep recip_hFacC in the diffusive flux calculation,
470  C-    but no hFacC in TKE volume control  C-    but no hFacC in TKE volume control
471  C-    No need for maskC(k-1) with recip_hFacC(k-1)  C-    No need for maskC(k-1) with recip_hFacC(k-1)
472           a(i,j,k) = -deltaTggl90           a3d(i,j,k) = -deltaTggl90
473       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)
474       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))
475       &        *recip_drC(k)*maskC(i,j,k,bi,bj)       &        *recip_drC(k)*maskC(i,j,k,bi,bj)
# Line 457  C-    No need for maskC(k-1) with recip_ Line 479  C-    No need for maskC(k-1) with recip_
479  C--   Upper diagonal  C--   Upper diagonal
480        DO j=jMin,jMax        DO j=jMin,jMax
481         DO i=iMin,iMax         DO i=iMin,iMax
482           c(i,j,1)  = 0. _d 0           c3d(i,j,1)  = 0. _d 0
483         ENDDO         ENDDO
484        ENDDO        ENDDO
485        DO k=2,Nr        DO k=2,Nr
# Line 467  C--   Upper diagonal Line 489  C--   Upper diagonal
489  C-    We keep recip_hFacC in the diffusive flux calculation,  C-    We keep recip_hFacC in the diffusive flux calculation,
490  C-    but no hFacC in TKE volume control  C-    but no hFacC in TKE volume control
491  C-    No need for maskC(k) with recip_hFacC(k)  C-    No need for maskC(k) with recip_hFacC(k)
492            c(i,j,k) = -deltaTggl90            c3d(i,j,k) = -deltaTggl90
493       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)
494       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))
495       &        *recip_drC(k)*maskC(i,j,k-1,bi,bj)       &        *recip_drC(k)*maskC(i,j,k-1,bi,bj)
# Line 479  C--   Center diagonal Line 501  C--   Center diagonal
501         km1 = MAX(k-1,1)         km1 = MAX(k-1,1)
502         DO j=jMin,jMax         DO j=jMin,jMax
503          DO i=iMin,iMax          DO i=iMin,iMax
504            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)
505       &        + ab15*deltaTggl90*GGL90ceps*SQRTTKE(I,J,K)       &        + implDissFac*deltaTggl90*GGL90ceps*SQRTTKE(i,j,k)
506       &        * rMixingLength(I,J,K)       &        * rMixingLength(i,j,k)
507       &        * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)       &        * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)
508           ENDDO           ENDDO
509         ENDDO         ENDDO
# Line 490  C     end set up matrix Line 512  C     end set up matrix
512    
513  C     Apply boundary condition  C     Apply boundary condition
514        kp1 = MIN(Nr,kSurf+1)        kp1 = MIN(Nr,kSurf+1)
515        DO J=jMin,jMax        DO j=jMin,jMax
516         DO I=iMin,iMax         DO i=iMin,iMax
517  C     estimate friction velocity uStar from surface forcing  C     estimate friction velocity uStar from surface forcing
518          uStarSquare = SQRT(          uStarSquare = SQRT(
519       &    ( .5 _d 0*( surfaceForcingU(I,  J,  bi,bj)       &    ( .5 _d 0*( surfaceForcingU(i,  j,  bi,bj)
520       &              + surfaceForcingU(I+1,J,  bi,bj) ) )**2       &              + surfaceForcingU(i+1,j,  bi,bj) ) )**2
521       &  + ( .5 _d 0*( surfaceForcingV(I,  J,  bi,bj)       &  + ( .5 _d 0*( surfaceForcingV(i,  j,  bi,bj)
522       &              + surfaceForcingV(I,  J+1,bi,bj) ) )**2       &              + surfaceForcingV(i,  j+1,bi,bj) ) )**2
523       &                     )       &                     )
524  C     Dirichlet surface boundary condition for TKE  C     Dirichlet surface boundary condition for TKE
525          gTKE(I,J,kSurf) = MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)          GGL90TKE(i,j,kSurf,bi,bj) = maskC(i,j,kSurf,bi,bj)
526       &                     *maskC(I,J,kSurf,bi,bj)       &           *MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)
527          gTKE(i,j,kp1) = gTKE(i,j,kp1)          GGL90TKE(i,j,kp1,bi,bj) = GGL90TKE(i,j,kp1,bi,bj)
528       &                - a(i,j,kp1)*gTKE(i,j,kSurf)       &               - a3d(i,j,kp1)*GGL90TKE(i,j,kSurf,bi,bj)
529          a(i,j,kp1) = 0. _d 0          a3d(i,j,kp1) = 0. _d 0
530  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom
531          kBottom   = MAX(kLowC(I,J,bi,bj),1)          kBottom   = MAX(kLowC(i,j,bi,bj),1)
532          gTKE(I,J,kBottom) = gTKE(I,J,kBottom)          GGL90TKE(i,j,kBottom,bi,bj) = GGL90TKE(i,j,kBottom,bi,bj)
533       &                    - GGL90TKEbottom*c(I,J,kBottom)       &                              - GGL90TKEbottom*c3d(i,j,kBottom)
534          c(I,J,kBottom) = 0. _d 0          c3d(i,j,kBottom) = 0. _d 0
535         ENDDO         ENDDO
536        ENDDO        ENDDO
537    
538  C     solve tri-diagonal system, and store solution on gTKE (previously rhs)  C     solve tri-diagonal system
539        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
540       I                        a, b, c,       I                        a3d, b3d, c3d,
541       U                        gTKE,       U                        GGL90TKE,
542       O                        errCode,       O                        errCode,
543       I                        bi, bj, myThid )       I                        bi, bj, myThid )
544  c     CALL GGL90_SOLVE( bi, bj, iMin, iMax, jMin, jMax,  
545  c    I     a, b, c,        DO k=1,Nr
546  c    U     gTKE,         DO j=jMin,jMax
547  c    I     myThid )          DO i=iMin,iMax
   
 C     now update TKE  
       DO K=1,Nr  
        DO J=jMin,jMax  
         DO I=iMin,iMax  
548  C     impose minimum TKE to avoid numerical undershoots below zero  C     impose minimum TKE to avoid numerical undershoots below zero
549           GGL90TKE(I,J,K,bi,bj) = MAX( gTKE(I,J,K), GGL90TKEmin )           GGL90TKE(i,j,k,bi,bj) = maskC(i,j,k,bi,bj)
550       &        * maskC(I,J,K,bi,bj)       &                  *MAX( GGL90TKE(i,j,k,bi,bj), GGL90TKEmin )
551          ENDDO          ENDDO
552         ENDDO         ENDDO
553        ENDDO        ENDDO
# Line 538  C     impose minimum TKE to avoid numeri Line 555  C     impose minimum TKE to avoid numeri
555  C     end of time step  C     end of time step
556  C     ===============================  C     ===============================
557    
558          DO k=2,Nr
559           DO j=1,sNy
560            DO i=1,sNx
561  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
562        DO K=1,Nr           tmpVisc=
        DO J=jMin,jMax  
         DO I=iMin,iMax  
          tmpdiffKrS=  
563       &  (       &  (
564       &   p4 *  GGL90viscAr(i  ,j  ,k,bi,bj) * mskCor(i  ,j  ,bi,bj)       &   p4 *  GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
565       &  +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)
566       &       + GGL90viscAr(i  ,j-1,k,bi,bj) * mskCor(i  ,j-1,bi,bj)       &       + GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj)
567       &       + GGL90viscAr(i+1,j  ,k,bi,bj) * mskCor(i+1,j  ,bi,bj)       &       + GGL90visctmp(i+1,j  ,k) * mskCor(i+1,j  ,bi,bj)
568       &       + GGL90viscAr(i  ,j+1,k,bi,bj) * mskCor(i  ,j+1,bi,bj))       &       + GGL90visctmp(i  ,j+1,k) * mskCor(i  ,j+1,bi,bj))
569       &  +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)
570       &       + 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)
571       &       + 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)
572       &       + 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))
573       &  )       &  )
574       & /(p4       & /(p4
575       &  +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 564  C     =============================== Line 581  C     ===============================
581       &       +       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)
582       &       +       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))
583       &  )*maskC(i,j,k,bi,bj)*mskCor(i,j,bi,bj)       &  )*maskC(i,j,k,bi,bj)*mskCor(i,j,bi,bj)
584       &   /TKEPrandtlNumber(i,j,k)  #else
585           GGL90diffKrS(I,J,K,bi,bj)= MAX( tmpdiffKrS , diffKrNrT(k) )           tmpVisc = GGL90visctmp(i,j,k)
586    #endif
587             tmpVisc = MIN(tmpVisc/TKEPrandtlNumber(i,j,k),GGL90diffMax)
588             GGL90diffKr(i,j,k,bi,bj)= MAX( tmpVisc , diffKrNrT(k) )
589            ENDDO
590           ENDDO
591          ENDDO
592    
593          DO k=2,Nr
594           DO j=1,sNy
595            DO i=1,sNx+1
596    #ifdef ALLOW_GGL90_SMOOTH
597            tmpVisc =
598         & (
599         &   p4 *(GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
600         &       +GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj))
601         &  +p8 *(GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj)
602         &       +GGL90visctmp(i-1,j+1,k) * mskCor(i-1,j+1,bi,bj)
603         &       +GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj)
604         &       +GGL90visctmp(i  ,j+1,k) * mskCor(i  ,j+1,bi,bj))
605         &  )
606         & /(p4 * 2. _d 0
607         &  +p8 *(      maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj)
608         &       +      maskC(i-1,j+1,k,bi,bj) * mskCor(i-1,j+1,bi,bj)
609         &       +      maskC(i  ,j-1,k,bi,bj) * mskCor(i  ,j-1,bi,bj)
610         &       +      maskC(i  ,j+1,k,bi,bj) * mskCor(i  ,j+1,bi,bj))
611         &  )
612         &  *maskC(i  ,j,k,bi,bj)*mskCor(i  ,j,bi,bj)
613         &  *maskC(i-1,j,k,bi,bj)*mskCor(i-1,j,bi,bj)
614    #else
615            tmpVisc = _maskW(i,j,k,bi,bj) *
616         &                   (.5 _d 0*(GGL90visctmp(i,j,k)
617         &                            +GGL90visctmp(i-1,j,k))
618         &                   )
619    #endif
620            tmpVisc = MIN( tmpVisc , GGL90viscMax )
621            GGL90viscArU(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
622          ENDDO          ENDDO
623         ENDDO         ENDDO
624        ENDDO        ENDDO
625    
626          DO k=2,Nr
627           DO j=1,sNy+1
628            DO i=1,sNx
629    #ifdef ALLOW_GGL90_SMOOTH
630            tmpVisc =
631         & (
632         &   p4 *(GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
633         &       +GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj))
634         &  +p8 *(GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj)
635         &       +GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj)
636         &       +GGL90visctmp(i+1,j  ,k) * mskCor(i+1,j  ,bi,bj)
637         &       +GGL90visctmp(i+1,j-1,k) * mskCor(i+1,j-1,bi,bj))
638         &  )
639         & /(p4 * 2. _d 0
640         &  +p8 *(      maskC(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)
641         &       +      maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj)
642         &       +      maskC(i+1,j  ,k,bi,bj) * mskCor(i+1,j  ,bi,bj)
643         &       +      maskC(i+1,j-1,k,bi,bj) * mskCor(i+1,j-1,bi,bj))
644         &  )
645         &   *maskC(i,j  ,k,bi,bj)*mskCor(i,j  ,bi,bj)
646         &   *maskC(i,j-1,k,bi,bj)*mskCor(i,j-1,bi,bj)
647    #else
648            tmpVisc = _maskS(i,j,k,bi,bj) *
649         &                   (.5 _d 0*(GGL90visctmp(i,j,k)
650         &                            +GGL90visctmp(i,j-1,k))
651         &                   )
652    
653  #endif  #endif
654            tmpVisc = MIN( tmpVisc , GGL90viscMax )
655            GGL90viscArV(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
656            ENDDO
657           ENDDO
658          ENDDO
659    
660  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
661        IF ( useDiagnostics ) THEN        IF ( useDiagnostics ) THEN
662           CALL DIAGNOSTICS_FILL( GGL90TKE   ,'GGL90TKE',           CALL DIAGNOSTICS_FILL( GGL90TKE   ,'GGL90TKE',
663       &                          0,Nr, 1, bi, bj, myThid )       &                          0,Nr, 1, bi, bj, myThid )
664           CALL DIAGNOSTICS_FILL( GGL90viscAr,'GGL90Ar ',           CALL DIAGNOSTICS_FILL( GGL90viscArU,'GGL90ArU',
665         &                          0,Nr, 1, bi, bj, myThid )
666             CALL DIAGNOSTICS_FILL( GGL90viscArV,'GGL90ArV',
667       &                          0,Nr, 1, bi, bj, myThid )       &                          0,Nr, 1, bi, bj, myThid )
668           CALL DIAGNOSTICS_FILL( GGL90diffKr,'GGL90Kr ',           CALL DIAGNOSTICS_FILL( GGL90diffKr,'GGL90Kr ',
669       &                          0,Nr, 1, bi, bj, myThid )       &                          0,Nr, 1, bi, bj, myThid )
Legend:
Removed from v.1.15  
changed lines
  Added in v.1.24
  ViewVC Help
Powered by ViewVC 1.1.22