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