/[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.30 by jmc, Sat Feb 21 20:11:20 2015 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  C !DESCRIPTION: \bv  C !DESCRIPTION: \bv
14  C     *==========================================================*  C     *==========================================================*
# Line 41  C !USES: =============================== Line 41  C !USES: ===============================
41    
42  C !INPUT PARAMETERS: ===================================================  C !INPUT PARAMETERS: ===================================================
43  C Routine arguments  C Routine arguments
44  C     bi, bj :: array indices on which to apply calculations  C     bi, bj :: Current tile indices
45    C     sigmaR :: Vertical gradient of iso-neutral density
46  C     myTime :: Current time in simulation  C     myTime :: Current time in simulation
47    C     myIter :: Current time-step number
48  C     myThid :: My Thread Id number  C     myThid :: My Thread Id number
49        INTEGER bi, bj        INTEGER bi, bj
50          _RL     sigmaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
51        _RL     myTime        _RL     myTime
52          INTEGER myIter
53        INTEGER myThid        INTEGER myThid
 CEOP  
54    
55  #ifdef ALLOW_GGL90  #ifdef ALLOW_GGL90
56    
57  C !LOCAL VARIABLES: ====================================================  C !LOCAL VARIABLES: ====================================================
58  C Local constants  C Local constants
59  C     iMin, iMax, jMin, jMax, I, J - array computation indices  C     iMin,iMax,jMin,jMax :: index boundaries of computation domain
60  C     K, Kp1, km1, kSurf, kBottom  - vertical loop indices  C     i, j, k, kp1,km1 :: array computation indices
61  C     ab15, ab05       - weights for implicit timestepping  C     kSurf, kBottom   :: vertical indices of domain boundaries
62  C     uStarSquare      - square of friction velocity  C     explDissFac      :: explicit Dissipation Factor (in [0-1])
63  C     verticalShear    - (squared) vertical shear of horizontal velocity  C     implDissFac      :: implicit Dissipation Factor (in [0-1])
64  C     Nsquare          - squared buoyancy freqency  C     uStarSquare      :: square of friction velocity
65  C     RiNumber         - local Richardson number  C     verticalShear    :: (squared) vertical shear of horizontal velocity
66  C     KappaM           - (local) viscosity parameter (eq.10)  C     Nsquare          :: squared buoyancy freqency
67  C     KappaH           - (local) diffusivity parameter for temperature (eq.11)  C     RiNumber         :: local Richardson number
68  C     KappaE           - (local) diffusivity parameter for TKE (eq.15)  C     KappaM           :: (local) viscosity parameter (eq.10)
69  C     TKEdissipation   - dissipation of TKE  C     KappaH           :: (local) diffusivity parameter for temperature (eq.11)
70  C     GGL90mixingLength- mixing length of scheme following Banke+Delecuse  C     KappaE           :: (local) diffusivity parameter for TKE (eq.15)
71  C         rMixingLength- inverse of mixing length  C     TKEdissipation   :: dissipation of TKE
72  C     totalDepth       - thickness of water column (inverse of recip_Rcol)  C     GGL90mixingLength:: mixing length of scheme following Banke+Delecuse
73  C     TKEPrandtlNumber - here, an empirical function of the Richardson number  C         rMixingLength:: inverse of mixing length
74  C     rhoK, rhoKm1     - density at layer K and Km1 (relative to K)  C     totalDepth       :: thickness of water column (inverse of recip_Rcol)
75  C     gTKE             - right hand side of implicit equation  C     TKEPrandtlNumber :: here, an empirical function of the Richardson number
76        INTEGER iMin ,iMax ,jMin ,jMax        INTEGER iMin ,iMax ,jMin ,jMax
77        INTEGER I, J, K, Kp1, Km1, kSurf, kBottom        INTEGER i, j, k, kp1, km1, kSurf, kBottom
78        _RL     ab15, ab05        _RL     explDissFac, implDissFac
79        _RL     uStarSquare        _RL     uStarSquare
80        _RL     verticalShear        _RL     verticalShear
81        _RL     KappaM, KappaH        _RL     KappaM, KappaH
# Line 82  c     _RL     Nsquare Line 85  c     _RL     Nsquare
85  c     _RL     SQRTTKE  c     _RL     SQRTTKE
86        _RL     SQRTTKE(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     SQRTTKE(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
87        _RL     RiNumber        _RL     RiNumber
88    #ifdef ALLOW_GGL90_IDEMIX
89          _RL     IDEMIX_RiNumber
90    #endif
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)
97        _RL     mxLength_Dn      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _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)
       _RL     rhoK             (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL     rhoKm1           (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
99        _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
100        _RL     gTKE             (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     GGL90visctmp     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
101    #ifdef ALLOW_DIAGNOSTICS
102          _RL     surf_flx_tke     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
103    #endif /* ALLOW_DIAGNOSTICS */
104  C-    tri-diagonal matrix  C-    tri-diagonal matrix
105        _RL     a(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     a3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
106        _RL     b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     b3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
107        _RL     c(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     c3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
108        INTEGER errCode        INTEGER errCode
109  #ifdef ALLOW_GGL90_HORIZDIFF  #ifdef ALLOW_GGL90_HORIZDIFF
110  C-    xA, yA   - area of lateral faces  C     hFac     :: fractional thickness of W-cell
111  C-    dfx, dfy - diffusive flux across lateral faces  C     xA, yA   :: area of lateral faces
112    C     dfx, dfy :: diffusive flux across lateral faces
113    C     gTKE     :: right hand side of diffusion equation
114          _RL     hFac
115        _RL     xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
116        _RL     yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
117        _RL     dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
118        _RL     dfy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     dfy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
119          _RL    gTKE(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
120  #endif /* ALLOW_GGL90_HORIZDIFF */  #endif /* ALLOW_GGL90_HORIZDIFF */
121  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
122        _RL p4, p8, p16, tmpdiffKrS        _RL p4, p8, p16
123    CEOP
124        p4=0.25 _d 0        p4=0.25 _d 0
125        p8=0.125 _d 0        p8=0.125 _d 0
126        p16=0.0625 _d 0        p16=0.0625 _d 0
# Line 122  C     set separate time step (should be Line 134  C     set separate time step (should be
134        deltaTggl90 = dTtracerLev(1)        deltaTggl90 = dTtracerLev(1)
135    
136        kSurf = 1        kSurf = 1
137  C     implicit timestepping weights for dissipation  C     explicit/implicit timestepping weights for dissipation
138        ab15 =  1.5 _d 0        explDissFac = 0. _d 0
139        ab05 = -0.5 _d 0        implDissFac = 1. _d 0 - explDissFac
       ab15 =  1. _d 0  
       ab05 =  0. _d 0  
140    
141  C     Initialize local fields  C     Initialize local fields
142        DO K = 1, Nr        DO k = 1, Nr
143         DO J=1-Oly,sNy+Oly         DO j=1-OLy,sNy+OLy
144          DO I=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
145           gTKE(I,J,K)              = 0. _d 0           rMixingLength(i,j,k)     = 0. _d 0
146           KappaE(I,J,K)            = 0. _d 0           mxLength_Dn(i,j,k)       = 0. _d 0
147           TKEPrandtlNumber(I,J,K)  = 1. _d 0           GGL90visctmp(i,j,k)      = 0. _d 0
148           GGL90mixingLength(I,J,K) = GGL90mixingLengthMin           KappaE(i,j,k)            = 0. _d 0
149             TKEPrandtlNumber(i,j,k)  = 1. _d 0
150             GGL90mixingLength(i,j,k) = GGL90mixingLengthMin
151             GGL90visctmp(i,j,k)      = 0. _d 0
152    #ifndef SOLVE_DIAGONAL_LOWMEMORY
153             a3d(i,j,k) = 0. _d 0
154             b3d(i,j,k) = 1. _d 0
155             c3d(i,j,k) = 0. _d 0
156    #endif
157             Nsquare(i,j,k) = 0. _d 0
158             SQRTTKE(i,j,k) = 0. _d 0
159          ENDDO          ENDDO
160         ENDDO         ENDDO
161        ENDDO        ENDDO
162        DO J=1-Oly,sNy+Oly        DO j=1-OLy,sNy+OLy
163         DO I=1-Olx,sNx+Olx         DO i=1-OLx,sNx+OLx
164          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)  
165          rMixingLength(i,j,1) = 0. _d 0          rMixingLength(i,j,1) = 0. _d 0
166          mxLength_Dn(I,J,1) = GGL90mixingLengthMin          mxLength_Dn(i,j,1) = GGL90mixingLengthMin
167          SQRTTKE(i,j,1) = SQRT( GGL90TKE(i,j,1,bi,bj) )          SQRTTKE(i,j,1) = SQRT( GGL90TKE(i,j,1,bi,bj) )
168    #ifdef ALLOW_GGL90_HORIZDIFF
169            xA(i,j)  = 0. _d 0
170            yA(i,j)  = 0. _d 0
171            dfx(i,j) = 0. _d 0
172            dfy(i,j) = 0. _d 0
173            gTKE(i,j) = 0. _d 0
174    #endif /* ALLOW_GGL90_HORIZDIFF */
175         ENDDO         ENDDO
176        ENDDO        ENDDO
177    
178    #ifdef ALLOW_GGL90_IDEMIX
179          IF ( useIDEMIX) CALL GGL90_IDEMIX(
180         &   bi, bj, sigmaR, myTime, myIter, myThid )
181    #endif /* ALLOW_GGL90_IDEMIX */
182    
183  C     start k-loop  C     start k-loop
184        DO K = 2, Nr        DO k = 2, Nr
185         Km1 = K-1  c      km1 = k-1
186  c      Kp1 = MIN(Nr,K+1)  c      kp1 = MIN(Nr,k+1)
187         CALL FIND_RHO_2D(         DO j=jMin,jMax
188       I      iMin, iMax, jMin, jMax, K,          DO i=iMin,iMax
189       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) )  
190    
191  C     buoyancy frequency  C     buoyancy frequency
192           Nsquare(i,j,k) = - gravity*recip_rhoConst*recip_drC(K)           Nsquare(i,j,k) = gravity*gravitySign*recip_rhoConst
193       &        * ( rhoKm1(I,J) - rhoK(I,J) )*maskC(I,J,K,bi,bj)       &                  * sigmaR(i,j,k)
194  cC     vertical shear term (dU/dz)^2+(dV/dz)^2  cC     vertical shear term (dU/dz)^2+(dV/dz)^2
195  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)
196  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)) )
197  c     &        *recip_drC(K)  c     &        *recip_drC(k)
198  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)
199  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)) )
200  c     &        *recip_drC(K)  c     &        *recip_drC(k)
201  c         verticalShear = tempU*tempU + tempV*tempV  c         verticalShear = tempU*tempU + tempV*tempV
202  c         RiNumber   = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)  c         RiNumber   = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)
203  cC     compute Prandtl number (always greater than 0)  cC     compute Prandtl number (always greater than 0)
204  c         prTemp = 1. _d 0  c         prTemp = 1. _d 0
205  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
206  c         TKEPrandtlNumber(I,J,K) = MIN(10. _d 0,prTemp)  c         TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
207  C     mixing length  C     mixing length
208           GGL90mixingLength(I,J,K) = SQRTTWO *           GGL90mixingLength(i,j,k) = SQRTTWO *
209       &        SQRTTKE(i,j,k)/SQRT( MAX(Nsquare(i,j,k),GGL90eps) )       &        SQRTTKE(i,j,k)/SQRT( MAX(Nsquare(i,j,k),GGL90eps) )
210          ENDDO          ENDDO
211         ENDDO         ENDDO
212        ENDDO        ENDDO
213    
214    C- ensure mixing between first and second level
215          IF (mxlSurfFlag) THEN
216           DO j=jMin,jMax
217            DO i=iMin,iMax
218             GGL90mixingLength(i,j,2)=drF(1)
219            ENDDO
220           ENDDO
221          ENDIF
222    
223  C- Impose upper and lower bound for mixing length  C- Impose upper and lower bound for mixing length
224        IF ( mxlMaxFlag .EQ. 0 ) THEN        IF ( mxlMaxFlag .EQ. 0 ) 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=totalDepth(I,J)            MaxLength=totalDepth(i,j)
230            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
231       &                                   MaxLength)       &                                   MaxLength)
232           ENDDO           ENDDO
233          ENDDO          ENDDO
234         ENDDO         ENDDO
235    
236         DO k=2,Nr         DO k=2,Nr
237          DO J=jMin,jMax          DO j=jMin,jMax
238           DO I=iMin,iMax           DO i=iMin,iMax
239            GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
240       &                                   GGL90mixingLengthMin)       &                                   GGL90mixingLengthMin)
241            rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K)            rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
242           ENDDO           ENDDO
243          ENDDO          ENDDO
244         ENDDO         ENDDO
245    
246        ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN
247  C-  
248         DO k=2,Nr         DO k=2,Nr
249          DO J=jMin,jMax          DO j=jMin,jMax
250           DO I=iMin,iMax           DO i=iMin,iMax
251            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))
252  c         MaxLength=MAX(MaxLength,20. _d 0)  c         MaxLength=MAX(MaxLength,20. _d 0)
253            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
254       &                                   MaxLength)       &                                   MaxLength)
255           ENDDO           ENDDO
256          ENDDO          ENDDO
257         ENDDO         ENDDO
258    
259         DO k=2,Nr         DO k=2,Nr
260          DO J=jMin,jMax          DO j=jMin,jMax
261           DO I=iMin,iMax           DO i=iMin,iMax
262            GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
263       &                                   GGL90mixingLengthMin)       &                                   GGL90mixingLengthMin)
264            rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K)            rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
265           ENDDO           ENDDO
266          ENDDO          ENDDO
267         ENDDO         ENDDO
268    
269        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN
270  C-  
271         DO k=2,Nr         DO k=2,Nr
272          DO J=jMin,jMax          DO j=jMin,jMax
273           DO I=iMin,iMax           DO i=iMin,iMax
274            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
275       &        GGL90mixingLength(I,J,K-1)+drF(k-1))       &        GGL90mixingLength(i,j,k-1)+drF(k-1))
276           ENDDO           ENDDO
277          ENDDO          ENDDO
278         ENDDO         ENDDO
279         DO J=jMin,jMax         DO j=jMin,jMax
280          DO I=iMin,iMax          DO i=iMin,iMax
281            GGL90mixingLength(I,J,Nr) = MIN(GGL90mixingLength(I,J,Nr),            GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
282       &       GGL90mixingLengthMin+drF(Nr))       &       GGL90mixingLengthMin+drF(Nr))
283          ENDDO          ENDDO
284         ENDDO         ENDDO
285         DO k=Nr-1,2,-1         DO k=Nr-1,2,-1
286          DO J=jMin,jMax          DO j=jMin,jMax
287           DO I=iMin,iMax           DO i=iMin,iMax
288            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
289       &        GGL90mixingLength(I,J,K+1)+drF(k))       &        GGL90mixingLength(i,j,k+1)+drF(k))
290           ENDDO           ENDDO
291          ENDDO          ENDDO
292         ENDDO         ENDDO
293    
294         DO k=2,Nr         DO k=2,Nr
295          DO J=jMin,jMax          DO j=jMin,jMax
296           DO I=iMin,iMax           DO i=iMin,iMax
297            GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
298       &                                   GGL90mixingLengthMin)       &                                   GGL90mixingLengthMin)
299            rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K)            rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
300           ENDDO           ENDDO
301          ENDDO          ENDDO
302         ENDDO         ENDDO
303    
304        ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN
305  C-  
306         DO k=2,Nr         DO k=2,Nr
307          DO J=jMin,jMax          DO j=jMin,jMax
308           DO I=iMin,iMax           DO i=iMin,iMax
309            mxLength_Dn(I,J,K) = MIN(GGL90mixingLength(I,J,K),            mxLength_Dn(i,j,k) = MIN(GGL90mixingLength(i,j,k),
310       &        mxLength_Dn(I,J,K-1)+drF(k-1))       &        mxLength_Dn(i,j,k-1)+drF(k-1))
311           ENDDO           ENDDO
312          ENDDO          ENDDO
313         ENDDO         ENDDO
314         DO J=jMin,jMax         DO j=jMin,jMax
315          DO I=iMin,iMax          DO i=iMin,iMax
316            GGL90mixingLength(I,J,Nr) = MIN(GGL90mixingLength(I,J,Nr),            GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
317       &       GGL90mixingLengthMin+drF(Nr))       &       GGL90mixingLengthMin+drF(Nr))
318          ENDDO          ENDDO
319         ENDDO         ENDDO
320         DO k=Nr-1,2,-1         DO k=Nr-1,2,-1
321          DO J=jMin,jMax          DO j=jMin,jMax
322           DO I=iMin,iMax           DO i=iMin,iMax
323            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
324       &        GGL90mixingLength(I,J,K+1)+drF(k))       &        GGL90mixingLength(i,j,k+1)+drF(k))
325           ENDDO           ENDDO
326          ENDDO          ENDDO
327         ENDDO         ENDDO
328    
329         DO k=2,Nr         DO k=2,Nr
330          DO J=jMin,jMax          DO j=jMin,jMax
331           DO I=iMin,iMax           DO i=iMin,iMax
332            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
333       &                                  mxLength_Dn(I,J,K))       &                                  mxLength_Dn(i,j,k))
334            tmpmlx = SQRT( GGL90mixingLength(I,J,K)*mxLength_Dn(I,J,K) )            tmpmlx = SQRT( GGL90mixingLength(i,j,k)*mxLength_Dn(i,j,k) )
335            tmpmlx = MAX( tmpmlx, GGL90mixingLengthMin)            tmpmlx = MAX( tmpmlx, GGL90mixingLengthMin)
336            rMixingLength(I,J,K) = 1. _d 0 / tmpmlx            rMixingLength(i,j,k) = 1. _d 0 / tmpmlx
337           ENDDO           ENDDO
338          ENDDO          ENDDO
339         ENDDO         ENDDO
# Line 316  C- Line 344  C-
344    
345  C- Impose minimum mixing length (to avoid division by zero)  C- Impose minimum mixing length (to avoid division by zero)
346  c      DO k=2,Nr  c      DO k=2,Nr
347  c      DO J=jMin,jMax  c      DO j=jMin,jMax
348  c       DO I=iMin,iMax  c       DO i=iMin,iMax
349  c        GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),  c        GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
350  c    &        GGL90mixingLengthMin)  c    &        GGL90mixingLengthMin)
351  c        rMixingLength(I,J,K) = 1. _d 0 /GGL90mixingLength(I,J,K)  c        rMixingLength(i,j,k) = 1. _d 0 /GGL90mixingLength(i,j,k)
352  c       ENDDO  c       ENDDO
353  c      ENDDO  c      ENDDO
354  c     ENDDO  c     ENDDO
355    
356        DO k=2,Nr        DO k=2,Nr
357         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)  
358    
 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)  
359  #ifdef ALLOW_GGL90_HORIZDIFF  #ifdef ALLOW_GGL90_HORIZDIFF
360        IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN        IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
361         DO K=2,Nr  C     horizontal diffusion of TKE (requires an exchange in
362    C      do_fields_blocking_exchanges)
363  C     common factors  C     common factors
364          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
365           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
366            xA(i,j) = _dyG(i,j,bi,bj)            xA(i,j) = _dyG(i,j,bi,bj)*drC(k)*
367       &         *drF(k)*_hFacW(i,j,k,bi,bj)       &                 (min(.5 _d 0,_hFacW(i,j,k-1,bi,bj) ) +
368            yA(i,j) = _dxG(i,j,bi,bj)       &                  min(.5 _d 0,_hFacW(i,j,k  ,bi,bj) ) )
369       &         *drF(k)*_hFacS(i,j,k,bi,bj)            yA(i,j) = _dxG(i,j,bi,bj)*drC(k)*
370         &                 (min(.5 _d 0,_hFacS(i,j,k-1,bi,bj) ) +
371         &                  min(.5 _d 0,_hFacS(i,j,k  ,bi,bj) ) )
372           ENDDO           ENDDO
373          ENDDO          ENDDO
374  C     Compute diffusive fluxes  C     Compute diffusive fluxes
375  C     ... across x-faces  C     ... across x-faces
376          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
377           dfx(1-Olx,j)=0. _d 0           dfx(1-OLx,j)=0. _d 0
378           DO i=1-Olx+1,sNx+Olx           DO i=1-OLx+1,sNx+OLx
379            dfx(i,j) = -GGL90diffTKEh*xA(i,j)            dfx(i,j) = -GGL90diffTKEh*xA(i,j)
380       &      *_recip_dxC(i,j,bi,bj)       &      *_recip_dxC(i,j,bi,bj)
381       &      *(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))
382    #ifdef ISOTROPIC_COS_SCALING
383       &      *CosFacU(j,bi,bj)       &      *CosFacU(j,bi,bj)
384    #endif /* ISOTROPIC_COS_SCALING */
385           ENDDO           ENDDO
386          ENDDO          ENDDO
387  C     ... across y-faces  C     ... across y-faces
388          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
389           dfy(i,1-Oly)=0. _d 0           dfy(i,1-OLy)=0. _d 0
390          ENDDO          ENDDO
391          DO j=1-Oly+1,sNy+Oly          DO j=1-OLy+1,sNy+OLy
392           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
393            dfy(i,j) = -GGL90diffTKEh*yA(i,j)            dfy(i,j) = -GGL90diffTKEh*yA(i,j)
394       &      *_recip_dyC(i,j,bi,bj)       &      *_recip_dyC(i,j,bi,bj)
395       &      *(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 399  C     ... across y-faces
399           ENDDO           ENDDO
400          ENDDO          ENDDO
401  C     Compute divergence of fluxes  C     Compute divergence of fluxes
402          DO j=1-Oly,sNy+Oly-1          DO j=1-OLy,sNy+OLy-1
403           DO i=1-Olx,sNx+Olx-1           DO i=1-OLx,sNx+OLx-1
404            gTKE(i,j,k)=gTKE(i,j,k)  #ifdef ALLOW_GGL90_IDEMIX
405       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj)            gTKE(i,j) = -recip_drC(k)*recip_rA(i,j,bi,bj)
406       &         *( (dfx(i+1,j)-dfx(i,j))       &         *recip_hFacI(i,j,k,bi,bj)
407       &           +(dfy(i,j+1)-dfy(i,j))  #else
408       &           )*deltaTggl90            hFac = MIN(.5 _d 0,_hFacC(i,j,k-1,bi,bj) ) +
409         &           MIN(.5 _d 0,_hFacC(i,j,k  ,bi,bj) )
410              gTKE(i,j) = 0.0
411              IF ( hFac .ne. 0.0 )
412         &      gTKE(i,j) = -recip_drC(k)*recip_rA(i,j,bi,bj)/hFac
413    #endif
414         &         *((dfx(i+1,j)-dfx(i,j))
415         &          +(dfy(i,j+1)-dfy(i,j)) )
416           ENDDO           ENDDO
417          ENDDO          ENDDO
418  C       end of k-loop  C      end if GGL90diffTKEh .eq. 0.
419           ENDIF
420    #endif /* ALLOW_GGL90_HORIZDIFF */
421    
422           DO j=jMin,jMax
423            DO i=iMin,iMax
424    C     vertical shear term (dU/dz)^2+(dV/dz)^2
425             tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)
426         &                 -( uVel(i,j,k  ,bi,bj)+uVel(i+1,j,k  ,bi,bj)) )
427         &        *recip_drC(k)
428    #ifdef ALLOW_GGL90_IDEMIX
429         &        *recip_hFacI(i,j,k,bi,bj)
430    #endif
431             tempV= .5 _d 0*( vVel(i,j,km1,bi,bj)+vVel(i,j+1,km1,bi,bj)
432         &                 -( vVel(i,j,k  ,bi,bj)+vVel(i,j+1,k  ,bi,bj)) )
433         &        *recip_drC(k)
434    #ifdef ALLOW_GGL90_IDEMIX
435         &        *recip_hFacI(i,j,k,bi,bj)
436    #endif
437             verticalShear = tempU*tempU + tempV*tempV
438    
439    C     viscosity and diffusivity
440             KappaM = GGL90ck*GGL90mixingLength(i,j,k)*SQRTTKE(i,j,k)
441             GGL90visctmp(i,j,k) = MAX(KappaM,diffKrNrT(k))
442         &                            * maskC(i,j,k,bi,bj)
443    C        note: storing GGL90visctmp like this, and using it later to compute
444    C              GGL9rdiffKr etc. is robust in case of smoothing (e.g. see OPA)
445             KappaM = MAX(KappaM,viscArNr(k)) * maskC(i,j,k,bi,bj)
446    
447    C     compute Prandtl number (always greater than 0)
448             RiNumber = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)
449    #ifdef ALLOW_GGL90_IDEMIX
450    CML         IDEMIX_RiNumber = 1./GGL90eps
451             IDEMIX_RiNumber = MAX( KappaM*Nsquare(i,j,k), 0. _d 0)/
452         &    (GGL90eps+IDEMIX_tau_d(i,j,k,bi,bj)*IDEMIX_E(i,j,k,bi,bj)**2)
453             prTemp         = MIN(5.*RiNumber, 6.6 _d 0*IDEMIX_RiNumber)
454    #else
455             prTemp = 1. _d 0
456             IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
457    #endif /* ALLOW_GGL90_IDEMIX */
458             TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
459             TKEPrandtlNumber(i,j,k) = MAX( 1. _d 0,TKEPrandtlNumber(i,j,k))
460    
461    C        diffusivity
462             KappaH = KappaM/TKEPrandtlNumber(i,j,k)
463             KappaE(i,j,k) = GGL90alpha * KappaM * maskC(i,j,k,bi,bj)
464    
465    C     dissipation term
466             TKEdissipation = explDissFac*GGL90ceps
467         &        *SQRTTKE(i,j,k)*rMixingLength(i,j,k)
468         &        *GGL90TKE(i,j,k,bi,bj)
469    C     partial update with sum of explicit contributions
470             GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
471         &        + deltaTggl90*(
472         &        + KappaM*verticalShear
473         &        - KappaH*Nsquare(i,j,k)
474         &        - TKEdissipation
475    #ifdef ALLOW_GGL90_IDEMIX
476         &        + IDEMIX_tau_d(i,j,k,bi,bj)*IDEMIX_E(i,j,k,bi,bj)**2
477    #endif
478         &        )
479            ENDDO
480         ENDDO         ENDDO
481  C     end if GGL90diffTKEh .eq. 0.  
482        ENDIF  #ifdef ALLOW_GGL90_HORIZDIFF
483           IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
484    C--    Add horiz. diffusion tendency
485            DO j=jMin,jMax
486             DO i=iMin,iMax
487              GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
488         &                          + gTKE(i,j)*deltaTggl90
489             ENDDO
490            ENDDO
491           ENDIF
492  #endif /* ALLOW_GGL90_HORIZDIFF */  #endif /* ALLOW_GGL90_HORIZDIFF */
493    
494    C--   end of k loop
495          ENDDO
496    
497  C     ============================================  C     ============================================
498  C     Implicit time step to update TKE for k=1,Nr;  C     Implicit time step to update TKE for k=1,Nr;
499  C     TKE(Nr+1)=0 by default  C     TKE(Nr+1)=0 by default
# Line 437  C     set up matrix Line 502  C     set up matrix
502  C--   Lower diagonal  C--   Lower diagonal
503        DO j=jMin,jMax        DO j=jMin,jMax
504         DO i=iMin,iMax         DO i=iMin,iMax
505           a(i,j,1) = 0. _d 0           a3d(i,j,1) = 0. _d 0
506         ENDDO         ENDDO
507        ENDDO        ENDDO
508        DO k=2,Nr        DO k=2,Nr
# Line 447  C--   Lower diagonal Line 512  C--   Lower diagonal
512  C-    We keep recip_hFacC in the diffusive flux calculation,  C-    We keep recip_hFacC in the diffusive flux calculation,
513  C-    but no hFacC in TKE volume control  C-    but no hFacC in TKE volume control
514  C-    No need for maskC(k-1) with recip_hFacC(k-1)  C-    No need for maskC(k-1) with recip_hFacC(k-1)
515           a(i,j,k) = -deltaTggl90           a3d(i,j,k) = -deltaTggl90
516       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)
517       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))
518       &        *recip_drC(k)*maskC(i,j,k,bi,bj)       &        *recip_drC(k)*maskC(i,j,k,bi,bj)
519    #ifdef ALLOW_GGL90_IDEMIX
520         &        *recip_hFacI(i,j,k,bi,bj)
521    #endif
522          ENDDO          ENDDO
523         ENDDO         ENDDO
524        ENDDO        ENDDO
525  C--   Upper diagonal  C--   Upper diagonal
526        DO j=jMin,jMax        DO j=jMin,jMax
527         DO i=iMin,iMax         DO i=iMin,iMax
528           c(i,j,1)  = 0. _d 0           c3d(i,j,1)  = 0. _d 0
529         ENDDO         ENDDO
530        ENDDO        ENDDO
531        DO k=2,Nr        DO k=2,Nr
# Line 467  C--   Upper diagonal Line 535  C--   Upper diagonal
535  C-    We keep recip_hFacC in the diffusive flux calculation,  C-    We keep recip_hFacC in the diffusive flux calculation,
536  C-    but no hFacC in TKE volume control  C-    but no hFacC in TKE volume control
537  C-    No need for maskC(k) with recip_hFacC(k)  C-    No need for maskC(k) with recip_hFacC(k)
538            c(i,j,k) = -deltaTggl90            c3d(i,j,k) = -deltaTggl90
539       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)
540       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))
541       &        *recip_drC(k)*maskC(i,j,k-1,bi,bj)       &        *recip_drC(k)*maskC(i,j,k-1,bi,bj)
542    #ifdef ALLOW_GGL90_IDEMIX
543         &        *recip_hFacI(i,j,k,bi,bj)
544    #endif
545          ENDDO          ENDDO
546         ENDDO         ENDDO
547        ENDDO        ENDDO
548    
549          IF (.NOT.GGL90_dirichlet) THEN
550    C      Neumann bottom boundary condition for TKE: no flux from bottom
551           DO j=jMin,jMax
552            DO i=iMin,iMax
553             kBottom   = MAX(kLowC(i,j,bi,bj),1)
554             c3d(i,j,kBottom) = 0. _d 0
555            ENDDO
556           ENDDO
557          ENDIF
558    
559  C--   Center diagonal  C--   Center diagonal
560        DO k=1,Nr        DO k=1,Nr
561         km1 = MAX(k-1,1)         km1 = MAX(k-1,1)
562         DO j=jMin,jMax         DO j=jMin,jMax
563          DO i=iMin,iMax          DO i=iMin,iMax
564            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)
565       &        + ab15*deltaTggl90*GGL90ceps*SQRTTKE(I,J,K)       &        + implDissFac*deltaTggl90*GGL90ceps*SQRTTKE(i,j,k)
566       &        * rMixingLength(I,J,K)       &        * rMixingLength(i,j,k)
567       &        * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)       &        * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)
568           ENDDO           ENDDO
569         ENDDO         ENDDO
# Line 490  C     end set up matrix Line 572  C     end set up matrix
572    
573  C     Apply boundary condition  C     Apply boundary condition
574        kp1 = MIN(Nr,kSurf+1)        kp1 = MIN(Nr,kSurf+1)
575        DO J=jMin,jMax        DO j=jMin,jMax
576         DO I=iMin,iMax         DO i=iMin,iMax
577  C     estimate friction velocity uStar from surface forcing  C     estimate friction velocity uStar from surface forcing
578          uStarSquare = SQRT(          uStarSquare = SQRT(
579       &    ( .5 _d 0*( surfaceForcingU(I,  J,  bi,bj)       &    ( .5 _d 0*( surfaceForcingU(i,  j,  bi,bj)
580       &              + surfaceForcingU(I+1,J,  bi,bj) ) )**2       &              + surfaceForcingU(i+1,j,  bi,bj) ) )**2
581       &  + ( .5 _d 0*( surfaceForcingV(I,  J,  bi,bj)       &  + ( .5 _d 0*( surfaceForcingV(i,  j,  bi,bj)
582       &              + surfaceForcingV(I,  J+1,bi,bj) ) )**2       &              + surfaceForcingV(i,  j+1,bi,bj) ) )**2
583       &                     )       &                     )
584  C     Dirichlet surface boundary condition for TKE  C     Dirichlet surface boundary condition for TKE
585          gTKE(I,J,kSurf) = MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)          GGL90TKE(i,j,kSurf,bi,bj) = maskC(i,j,kSurf,bi,bj)
586       &                     *maskC(I,J,kSurf,bi,bj)       &           *MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)
587          gTKE(i,j,kp1) = gTKE(i,j,kp1)          GGL90TKE(i,j,kp1,bi,bj) = GGL90TKE(i,j,kp1,bi,bj)
588       &                - a(i,j,kp1)*gTKE(i,j,kSurf)       &               - a3d(i,j,kp1)*GGL90TKE(i,j,kSurf,bi,bj)
589          a(i,j,kp1) = 0. _d 0          a3d(i,j,kp1) = 0. _d 0
 C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom  
         kBottom   = MAX(kLowC(I,J,bi,bj),1)  
         gTKE(I,J,kBottom) = gTKE(I,J,kBottom)  
      &                    - GGL90TKEbottom*c(I,J,kBottom)  
         c(I,J,kBottom) = 0. _d 0  
590         ENDDO         ENDDO
591        ENDDO        ENDDO
592    
593  C     solve tri-diagonal system, and store solution on gTKE (previously rhs)        IF (GGL90_dirichlet) THEN
594    C      Dirichlet bottom boundary condition for TKE = GGL90TKEbottom
595           DO j=jMin,jMax
596            DO i=iMin,iMax
597             kBottom   = MAX(kLowC(i,j,bi,bj),1)
598             GGL90TKE(i,j,kBottom,bi,bj) = GGL90TKE(i,j,kBottom,bi,bj)
599         &                              - GGL90TKEbottom*c3d(i,j,kBottom)
600             c3d(i,j,kBottom) = 0. _d 0
601            ENDDO
602           ENDDO
603          ENDIF
604    
605    C     solve tri-diagonal system
606        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
607       I                        a, b, c,       I                        a3d, b3d, c3d,
608       U                        gTKE,       U                        GGL90TKE(1-OLx,1-OLy,1,bi,bj),
609       O                        errCode,       O                        errCode,
610       I                        bi, bj, myThid )       I                        bi, bj, myThid )
611  c     CALL GGL90_SOLVE( bi, bj, iMin, iMax, jMin, jMax,  
612  c    I     a, b, c,        DO k=1,Nr
613  c    U     gTKE,         DO j=jMin,jMax
614  c    I     myThid )          DO i=iMin,iMax
   
 C     now update TKE  
       DO K=1,Nr  
        DO J=jMin,jMax  
         DO I=iMin,iMax  
615  C     impose minimum TKE to avoid numerical undershoots below zero  C     impose minimum TKE to avoid numerical undershoots below zero
616           GGL90TKE(I,J,K,bi,bj) = MAX( gTKE(I,J,K), GGL90TKEmin )           GGL90TKE(i,j,k,bi,bj) = maskC(i,j,k,bi,bj)
617       &        * maskC(I,J,K,bi,bj)       &                  *MAX( GGL90TKE(i,j,k,bi,bj), GGL90TKEmin )
618          ENDDO          ENDDO
619         ENDDO         ENDDO
620        ENDDO        ENDDO
# Line 538  C     impose minimum TKE to avoid numeri Line 622  C     impose minimum TKE to avoid numeri
622  C     end of time step  C     end of time step
623  C     ===============================  C     ===============================
624    
625          DO k=2,Nr
626           DO j=1,sNy
627            DO i=1,sNx
628  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
629        DO K=1,Nr           tmpVisc=
        DO J=jMin,jMax  
         DO I=iMin,iMax  
          tmpdiffKrS=  
630       &  (       &  (
631       &   p4 *  GGL90viscAr(i  ,j  ,k,bi,bj) * mskCor(i  ,j  ,bi,bj)       &   p4 *  GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
632       &  +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)
633       &       + GGL90viscAr(i  ,j-1,k,bi,bj) * mskCor(i  ,j-1,bi,bj)       &       + GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj)
634       &       + GGL90viscAr(i+1,j  ,k,bi,bj) * mskCor(i+1,j  ,bi,bj)       &       + GGL90visctmp(i+1,j  ,k) * mskCor(i+1,j  ,bi,bj)
635       &       + GGL90viscAr(i  ,j+1,k,bi,bj) * mskCor(i  ,j+1,bi,bj))       &       + GGL90visctmp(i  ,j+1,k) * mskCor(i  ,j+1,bi,bj))
636       &  +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)
637       &       + 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)
638       &       + 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)
639       &       + 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))
640       &  )       &  )
641       & /(p4       & /(p4
642       &  +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 648  C     ===============================
648       &       +       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)
649       &       +       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))
650       &  )*maskC(i,j,k,bi,bj)*mskCor(i,j,bi,bj)       &  )*maskC(i,j,k,bi,bj)*mskCor(i,j,bi,bj)
651       &   /TKEPrandtlNumber(i,j,k)  #else
652           GGL90diffKrS(I,J,K,bi,bj)= MAX( tmpdiffKrS , diffKrNrT(k) )           tmpVisc = GGL90visctmp(i,j,k)
653    #endif
654             tmpVisc = MIN(tmpVisc/TKEPrandtlNumber(i,j,k),GGL90diffMax)
655             GGL90diffKr(i,j,k,bi,bj)= MAX( tmpVisc , diffKrNrT(k) )
656          ENDDO          ENDDO
657         ENDDO         ENDDO
658        ENDDO        ENDDO
659    
660          DO k=2,Nr
661           DO j=1,sNy
662            DO i=1,sNx+1
663    #ifdef ALLOW_GGL90_SMOOTH
664            tmpVisc =
665         & (
666         &   p4 *(GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
667         &       +GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj))
668         &  +p8 *(GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj)
669         &       +GGL90visctmp(i-1,j+1,k) * mskCor(i-1,j+1,bi,bj)
670         &       +GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj)
671         &       +GGL90visctmp(i  ,j+1,k) * mskCor(i  ,j+1,bi,bj))
672         &  )
673         & /(p4 * 2. _d 0
674         &  +p8 *(      maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj)
675         &       +      maskC(i-1,j+1,k,bi,bj) * mskCor(i-1,j+1,bi,bj)
676         &       +      maskC(i  ,j-1,k,bi,bj) * mskCor(i  ,j-1,bi,bj)
677         &       +      maskC(i  ,j+1,k,bi,bj) * mskCor(i  ,j+1,bi,bj))
678         &  )
679         &  *maskC(i  ,j,k,bi,bj)*mskCor(i  ,j,bi,bj)
680         &  *maskC(i-1,j,k,bi,bj)*mskCor(i-1,j,bi,bj)
681    #else
682            tmpVisc = _maskW(i,j,k,bi,bj) *
683         &                   (.5 _d 0*(GGL90visctmp(i,j,k)
684         &                            +GGL90visctmp(i-1,j,k))
685         &                   )
686  #endif  #endif
687            tmpVisc = MIN( tmpVisc , GGL90viscMax )
688            GGL90viscArU(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
689            ENDDO
690           ENDDO
691          ENDDO
692    
693          DO k=2,Nr
694           DO j=1,sNy+1
695            DO i=1,sNx
696    #ifdef ALLOW_GGL90_SMOOTH
697            tmpVisc =
698         & (
699         &   p4 *(GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
700         &       +GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj))
701         &  +p8 *(GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj)
702         &       +GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj)
703         &       +GGL90visctmp(i+1,j  ,k) * mskCor(i+1,j  ,bi,bj)
704         &       +GGL90visctmp(i+1,j-1,k) * mskCor(i+1,j-1,bi,bj))
705         &  )
706         & /(p4 * 2. _d 0
707         &  +p8 *(      maskC(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)
708         &       +      maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj)
709         &       +      maskC(i+1,j  ,k,bi,bj) * mskCor(i+1,j  ,bi,bj)
710         &       +      maskC(i+1,j-1,k,bi,bj) * mskCor(i+1,j-1,bi,bj))
711         &  )
712         &   *maskC(i,j  ,k,bi,bj)*mskCor(i,j  ,bi,bj)
713         &   *maskC(i,j-1,k,bi,bj)*mskCor(i,j-1,bi,bj)
714    #else
715            tmpVisc = _maskS(i,j,k,bi,bj) *
716         &                   (.5 _d 0*(GGL90visctmp(i,j,k)
717         &                            +GGL90visctmp(i,j-1,k))
718         &                   )
719    
720    #endif
721            tmpVisc = MIN( tmpVisc , GGL90viscMax )
722            GGL90viscArV(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
723            ENDDO
724           ENDDO
725          ENDDO
726    
727  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
728        IF ( useDiagnostics ) THEN        IF ( useDiagnostics ) THEN
729           CALL DIAGNOSTICS_FILL( GGL90TKE   ,'GGL90TKE',          CALL DIAGNOSTICS_FILL( GGL90TKE   ,'GGL90TKE',
730       &                          0,Nr, 1, bi, bj, myThid )       &                         0,Nr, 1, bi, bj, myThid )
731           CALL DIAGNOSTICS_FILL( GGL90viscAr,'GGL90Ar ',          CALL DIAGNOSTICS_FILL( GGL90viscArU,'GGL90ArU',
732       &                          0,Nr, 1, bi, bj, myThid )       &                         0,Nr, 1, bi, bj, myThid )
733           CALL DIAGNOSTICS_FILL( GGL90diffKr,'GGL90Kr ',          CALL DIAGNOSTICS_FILL( GGL90viscArV,'GGL90ArV',
734       &                          0,Nr, 1, bi, bj, myThid )       &                         0,Nr, 1, bi, bj, myThid )
735           CALL DIAGNOSTICS_FILL( TKEPrandtlNumber ,'GGL90Prl',          CALL DIAGNOSTICS_FILL( GGL90diffKr,'GGL90Kr ',
736       &                          0,Nr, 2, bi, bj, myThid )       &                         0,Nr, 1, bi, bj, myThid )
737           CALL DIAGNOSTICS_FILL( GGL90mixingLength,'GGL90Lmx',          CALL DIAGNOSTICS_FILL( TKEPrandtlNumber ,'GGL90Prl',
738       &                          0,Nr, 2, bi, bj, myThid )       &                         0,Nr, 2, bi, bj, myThid )
739            CALL DIAGNOSTICS_FILL( GGL90mixingLength,'GGL90Lmx',
740         &                         0,Nr, 2, bi, bj, myThid )
741    
742            kp1 = MIN(Nr,kSurf+1)
743            DO j=jMin,jMax
744             DO i=iMin,iMax
745    C     diagnose surface flux of TKE
746              surf_flx_tke(i,j) =(GGL90TKE(i,j,kSurf,bi,bj)-
747         &                        GGL90TKE(i,j,kp1,bi,bj))
748         &        *recip_drF(kSurf)*recip_hFacC(i,j,kSurf,bi,bj)
749         &        *KappaE(i,j,kp1)
750             ENDDO
751            ENDDO
752            CALL DIAGNOSTICS_FILL( surf_flx_tke,'GGL90flx',
753         &                         0, 1, 2, bi, bj, myThid )
754    
755            k=kSurf
756            DO j=jMin,jMax
757             DO i=iMin,iMax
758    C     diagnose work done by the wind
759              surf_flx_tke(i,j) =
760         &      halfRL*( surfaceForcingU(i,  j,bi,bj)*uVel(i  ,j,k,bi,bj)
761         &              +surfaceForcingU(i+1,j,bi,bj)*uVel(i+1,j,k,bi,bj))
762         &    + halfRL*( surfaceForcingV(i,j,  bi,bj)*vVel(i,j  ,k,bi,bj)
763         &              +surfaceForcingV(i,j+1,bi,bj)*vVel(i,j+1,k,bi,bj))
764             ENDDO
765            ENDDO
766            CALL DIAGNOSTICS_FILL( surf_flx_tke,'GGL90tau',
767         &                         0, 1, 2, bi, bj, myThid )
768    
769        ENDIF        ENDIF
770  #endif  #endif /* ALLOW_DIAGNOSTICS */
771    
772  #endif /* ALLOW_GGL90 */  #endif /* ALLOW_GGL90 */
773    
Legend:
Removed from v.1.15  
changed lines
  Added in v.1.30
  ViewVC Help
Powered by ViewVC 1.1.22