/[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.18 by gforget, Wed Aug 11 03:32:29 2010 UTC revision 1.34 by jmc, Thu Jan 14 17:50:25 2016 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 35  C !USES: =============================== Line 35  C !USES: ===============================
35  #include "EEPARAMS.h"  #include "EEPARAMS.h"
36  #include "PARAMS.h"  #include "PARAMS.h"
37  #include "DYNVARS.h"  #include "DYNVARS.h"
 #include "GGL90.h"  
38  #include "FFIELDS.h"  #include "FFIELDS.h"
39  #include "GRID.h"  #include "GRID.h"
40    #include "GGL90.h"
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     hFac/hFacI       :: fractional thickness of W-cell
63  C     verticalShear    - (squared) vertical shear of horizontal velocity  C     explDissFac      :: explicit Dissipation Factor (in [0-1])
64  C     Nsquare          - squared buoyancy freqency  C     implDissFac      :: implicit Dissipation Factor (in [0-1])
65  C     RiNumber         - local Richardson number  C     uStarSquare      :: square of friction velocity
66  C     KappaM           - (local) viscosity parameter (eq.10)  C     verticalShear    :: (squared) vertical shear of horizontal velocity
67  C     KappaH           - (local) diffusivity parameter for temperature (eq.11)  C     Nsquare          :: squared buoyancy freqency
68  C     KappaE           - (local) diffusivity parameter for TKE (eq.15)  C     RiNumber         :: local Richardson number
69  C     TKEdissipation   - dissipation of TKE  C     KappaM           :: (local) viscosity parameter (eq.10)
70  C     GGL90mixingLength- mixing length of scheme following Banke+Delecuse  C     KappaH           :: (local) diffusivity parameter for temperature (eq.11)
71  C         rMixingLength- inverse of mixing length  C     KappaE           :: (local) diffusivity parameter for TKE (eq.15)
72  C     totalDepth       - thickness of water column (inverse of recip_Rcol)  C     TKEdissipation   :: dissipation of TKE
73  C     TKEPrandtlNumber - here, an empirical function of the Richardson number  C     GGL90mixingLength:: mixing length of scheme following Banke+Delecuse
74  C     rhoK, rhoKm1     - density at layer K and Km1 (relative to K)  C         rMixingLength:: inverse of mixing length
75  C     gTKE             - right hand side of implicit equation  C     totalDepth       :: thickness of water column (inverse of recip_Rcol)
76    C     TKEPrandtlNumber :: here, an empirical function of the Richardson number
77        INTEGER iMin ,iMax ,jMin ,jMax        INTEGER iMin ,iMax ,jMin ,jMax
78        INTEGER I, J, K, Kp1, Km1, kSurf, kBottom        INTEGER i, j, k, kp1, km1, kSurf, kBottom
79        _RL     ab15, ab05        _RL     explDissFac, implDissFac
80        _RL     uStarSquare        _RL     uStarSquare
81        _RL     verticalShear        _RL     verticalShear(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
82        _RL     KappaM, KappaH        _RL     KappaM(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
83          _RL     KappaH
84  c     _RL     Nsquare  c     _RL     Nsquare
85        _RL     Nsquare(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     Nsquare(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
86        _RL     deltaTggl90        _RL     deltaTggl90
87  c     _RL     SQRTTKE  c     _RL     SQRTTKE
88        _RL     SQRTTKE(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     SQRTTKE(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
89        _RL     RiNumber        _RL     RiNumber
90    #ifdef ALLOW_GGL90_IDEMIX
91          _RL     IDEMIX_RiNumber
92    #endif
93        _RL     TKEdissipation        _RL     TKEdissipation
94        _RL     tempU, tempV, prTemp        _RL     tempU, tempUp, tempV, tempVp, prTemp
95        _RL     MaxLength, tmpmlx, tmpVisc        _RL     MaxLength, tmpmlx, tmpVisc
96        _RL     TKEPrandtlNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     TKEPrandtlNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
97        _RL     GGL90mixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     GGL90mixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
98        _RL     rMixingLength    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     rMixingLength    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
99        _RL     mxLength_Dn      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     mxLength_Dn      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
100        _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)  
101        _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL     gTKE             (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
102        _RL     GGL90visctmp     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     GGL90visctmp     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
103    #ifdef ALLOW_DIAGNOSTICS
104          _RL     surf_flx_tke     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105    #endif /* ALLOW_DIAGNOSTICS */
106    C     hFac(I)  :: fractional thickness of W-cell
107          _RL       hFac
108    #ifdef ALLOW_GGL90_IDEMIX
109          _RL       hFacI(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
110    #endif /* ALLOW_GGL90_IDEMIX */
111          _RL recip_hFacI(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
112  C-    tri-diagonal matrix  C-    tri-diagonal matrix
113        _RL     a(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     a3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
114        _RL     b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     b3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
115        _RL     c(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     c3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
116        INTEGER errCode        INTEGER errCode
117  #ifdef ALLOW_GGL90_HORIZDIFF  #ifdef ALLOW_GGL90_HORIZDIFF
118  C-    xA, yA   - area of lateral faces  C     xA, yA   :: area of lateral faces
119  C-    dfx, dfy - diffusive flux across lateral faces  C     dfx, dfy :: diffusive flux across lateral faces
120    C     gTKE     :: right hand side of diffusion equation
121        _RL     xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
122        _RL     yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
123        _RL     dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
124        _RL     dfy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     dfy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
125          _RL    gTKE(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
126  #endif /* ALLOW_GGL90_HORIZDIFF */  #endif /* ALLOW_GGL90_HORIZDIFF */
127  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
128        _RL p4, p8, p16        _RL p4, p8, p16
       p4=0.25 _d 0  
       p8=0.125 _d 0  
       p16=0.0625 _d 0  
129  #endif  #endif
130        iMin = 2-OLx  CEOP
131        iMax = sNx+OLx-1  
132        jMin = 2-OLy        PARAMETER( iMin = 2-OLx, iMax = sNx+OLx-1 )
133        jMax = sNy+OLy-1        PARAMETER( jMin = 2-OLy, jMax = sNy+OLy-1 )
134    #ifdef ALLOW_GGL90_SMOOTH
135          p4  = 0.25   _d 0
136          p8  = 0.125  _d 0
137          p16 = 0.0625 _d 0
138    #endif
139    
140  C     set separate time step (should be deltaTtracer)  C     set separate time step (should be deltaTtracer)
141        deltaTggl90 = dTtracerLev(1)        deltaTggl90 = dTtracerLev(1)
142    
143        kSurf = 1        kSurf = 1
144  C     implicit timestepping weights for dissipation  C     explicit/implicit timestepping weights for dissipation
145        ab15 =  1.5 _d 0        explDissFac = 0. _d 0
146        ab05 = -0.5 _d 0        implDissFac = 1. _d 0 - explDissFac
147        ab15 =  1. _d 0  
148        ab05 =  0. _d 0  C     For nonlinear free surface and especially with r*-coordinates, the
149    C     hFacs change every timestep, so we need to update them here in the
150    C     case of using IDEMIX.
151           DO K=1,Nr
152            Km1 = MAX(K-1,1)
153            DO j=1-OLy,sNy+OLy
154             DO i=1-OLx,sNx+OLx
155              hFac =
156         &         MIN(.5 _d 0,_hFacC(i,j,km1,bi,bj) ) +
157         &         MIN(.5 _d 0,_hFacC(i,j,k  ,bi,bj) )
158              recip_hFacI(I,J,K)=0. _d 0
159              IF ( hFac .NE. 0. _d 0 )
160         &         recip_hFacI(I,J,K)=1. _d 0/hFac
161    #ifdef ALLOW_GGL90_IDEMIX
162              hFacI(i,j,k) = hFac
163    #endif /* ALLOW_GGL90_IDEMIX */
164             ENDDO
165            ENDDO
166           ENDDO
167    
168  C     Initialize local fields  C     Initialize local fields
169        DO K = 1, Nr        DO k = 1, Nr
170         DO J=1-Oly,sNy+Oly         DO j=1-OLy,sNy+OLy
171          DO I=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
172           gTKE(I,J,K)              = 0. _d 0           rMixingLength(i,j,k)     = 0. _d 0
173           KappaE(I,J,K)            = 0. _d 0           mxLength_Dn(i,j,k)       = 0. _d 0
174           TKEPrandtlNumber(I,J,K)  = 1. _d 0           GGL90visctmp(i,j,k)      = 0. _d 0
175           GGL90mixingLength(I,J,K) = GGL90mixingLengthMin           KappaE(i,j,k)            = 0. _d 0
176           GGL90visctmp(I,J,K)      = 0. _d 0           TKEPrandtlNumber(i,j,k)  = 1. _d 0
177             GGL90mixingLength(i,j,k) = GGL90mixingLengthMin
178             GGL90visctmp(i,j,k)      = 0. _d 0
179    #ifndef SOLVE_DIAGONAL_LOWMEMORY
180             a3d(i,j,k) = 0. _d 0
181             b3d(i,j,k) = 1. _d 0
182             c3d(i,j,k) = 0. _d 0
183    #endif
184             Nsquare(i,j,k) = 0. _d 0
185             SQRTTKE(i,j,k) = 0. _d 0
186          ENDDO          ENDDO
187         ENDDO         ENDDO
188        ENDDO        ENDDO
189        DO J=1-Oly,sNy+Oly        DO j=1-OLy,sNy+OLy
190         DO I=1-Olx,sNx+Olx         DO i=1-OLx,sNx+OLx
191          rhoK(I,J)          = 0. _d 0          KappaM(i,j)        = 0. _d 0
192          rhoKm1(I,J)        = 0. _d 0          verticalShear(i,j) = 0. _d 0
193          totalDepth(I,J)    = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)          totalDepth(i,j)    = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
194          rMixingLength(i,j,1) = 0. _d 0          rMixingLength(i,j,1) = 0. _d 0
195          mxLength_Dn(I,J,1) = GGL90mixingLengthMin          mxLength_Dn(i,j,1) = GGL90mixingLengthMin
196          SQRTTKE(i,j,1) = SQRT( GGL90TKE(i,j,1,bi,bj) )          SQRTTKE(i,j,1) = SQRT( GGL90TKE(i,j,1,bi,bj) )
197    #ifdef ALLOW_GGL90_HORIZDIFF
198            xA(i,j)  = 0. _d 0
199            yA(i,j)  = 0. _d 0
200            dfx(i,j) = 0. _d 0
201            dfy(i,j) = 0. _d 0
202            gTKE(i,j) = 0. _d 0
203    #endif /* ALLOW_GGL90_HORIZDIFF */
204         ENDDO         ENDDO
205        ENDDO        ENDDO
206    
207  C     start k-loop  #ifdef ALLOW_GGL90_IDEMIX
208        DO K = 2, Nr        IF ( useIDEMIX) CALL GGL90_IDEMIX(
209         Km1 = K-1       &   bi, bj, hFacI, recip_hFacI, sigmaR, myTime, myIter, myThid )
210  c      Kp1 = MIN(Nr,K+1)  #endif /* ALLOW_GGL90_IDEMIX */
211         CALL FIND_RHO_2D(  
212       I      iMin, iMax, jMin, jMax, K,        DO k = 2, Nr
213       I      theta(1-OLx,1-OLy,Km1,bi,bj), salt(1-OLx,1-OLy,Km1,bi,bj),         DO j=jMin,jMax
214       O      rhoKm1,          DO i=iMin,iMax
215       I      Km1, bi, bj, myThid )           SQRTTKE(i,j,k)=SQRT( GGL90TKE(i,j,k,bi,bj) )
   
        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) )  
216    
217  C     buoyancy frequency  C     buoyancy frequency
218           Nsquare(i,j,k) = - gravity*recip_rhoConst*recip_drC(K)           Nsquare(i,j,k) = gravity*gravitySign*recip_rhoConst
219       &        * ( rhoKm1(I,J) - rhoK(I,J) )*maskC(I,J,K,bi,bj)       &                  * sigmaR(i,j,k)
220  cC     vertical shear term (dU/dz)^2+(dV/dz)^2  C     vertical shear term (dU/dz)^2+(dV/dz)^2 is computed later
221  c         tempU= .5 _d 0*( uVel(I,J,Km1,bi,bj)+uVel(I+1,J,Km1,bi,bj)  C     to save some memory
 c     &                 -( uVel(I,J,K  ,bi,bj)+uVel(I+1,J,K  ,bi,bj)) )  
 c     &        *recip_drC(K)  
 c         tempV= .5 _d 0*( vVel(I,J,Km1,bi,bj)+vVel(I,J+1,Km1,bi,bj)  
 c     &                 -( vVel(I,J,K  ,bi,bj)+vVel(I,J+1,K  ,bi,bj)) )  
 c     &        *recip_drC(K)  
 c         verticalShear = tempU*tempU + tempV*tempV  
 c         RiNumber   = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)  
 cC     compute Prandtl number (always greater than 0)  
 c         prTemp = 1. _d 0  
 c         IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber  
 c         TKEPrandtlNumber(I,J,K) = MIN(10. _d 0,prTemp)  
222  C     mixing length  C     mixing length
223           GGL90mixingLength(I,J,K) = SQRTTWO *           GGL90mixingLength(i,j,k) = SQRTTWO *
224       &        SQRTTKE(i,j,k)/SQRT( MAX(Nsquare(i,j,k),GGL90eps) )       &        SQRTTKE(i,j,k)/SQRT( MAX(Nsquare(i,j,k),GGL90eps) )
225          ENDDO          ENDDO
226         ENDDO         ENDDO
227        ENDDO        ENDDO
228    
229  C- Impose upper and lower bound for mixing length  C- ensure mixing between first and second level
230          IF (mxlSurfFlag) THEN
231           DO j=jMin,jMax
232            DO i=iMin,iMax
233             GGL90mixingLength(i,j,2)=drF(1)
234            ENDDO
235           ENDDO
236          ENDIF
237    
238    C--   Impose upper and lower bound for mixing length
239    C--   Impose minimum mixing length to avoid division by zero
240        IF ( mxlMaxFlag .EQ. 0 ) THEN        IF ( mxlMaxFlag .EQ. 0 ) THEN
241  C-  
242         DO k=2,Nr         DO k=2,Nr
243          DO J=jMin,jMax          DO j=jMin,jMax
244           DO I=iMin,iMax           DO i=iMin,iMax
245            MaxLength=totalDepth(I,J)            MaxLength=totalDepth(i,j)
246            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
247       &                                   MaxLength)       &                                   MaxLength)
248           ENDDO           ENDDO
249          ENDDO          ENDDO
250         ENDDO         ENDDO
251    
252         DO k=2,Nr         DO k=2,Nr
253          DO J=jMin,jMax          DO j=jMin,jMax
254           DO I=iMin,iMax           DO i=iMin,iMax
255            GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
256       &                                   GGL90mixingLengthMin)       &                                   GGL90mixingLengthMin)
257            rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K)            rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
258           ENDDO           ENDDO
259          ENDDO          ENDDO
260         ENDDO         ENDDO
261    
262        ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN
263  C-  
264         DO k=2,Nr         DO k=2,Nr
265          DO J=jMin,jMax          DO j=jMin,jMax
266           DO I=iMin,iMax           DO i=iMin,iMax
267            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))
268  c         MaxLength=MAX(MaxLength,20. _d 0)  c         MaxLength=MAX(MaxLength,20. _d 0)
269            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
270       &                                   MaxLength)       &                                   MaxLength)
271           ENDDO           ENDDO
272          ENDDO          ENDDO
273         ENDDO         ENDDO
274    
275         DO k=2,Nr         DO k=2,Nr
276          DO J=jMin,jMax          DO j=jMin,jMax
277           DO I=iMin,iMax           DO i=iMin,iMax
278            GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
279       &                                   GGL90mixingLengthMin)       &                                   GGL90mixingLengthMin)
280            rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K)            rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
281           ENDDO           ENDDO
282          ENDDO          ENDDO
283         ENDDO         ENDDO
284    
285        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN
286  C-  
 cgf ensure mixing between first and second level  
 c      DO J=jMin,jMax  
 c        DO I=iMin,iMax  
 c         GGL90mixingLength(I,J,2)=drF(1)  
 c        ENDDO  
 c      ENDDO  
 cgf  
287         DO k=2,Nr         DO k=2,Nr
288          DO J=jMin,jMax          DO j=jMin,jMax
289           DO I=iMin,iMax           DO i=iMin,iMax
290            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
291       &        GGL90mixingLength(I,J,K-1)+drF(k-1))       &        GGL90mixingLength(i,j,k-1)+drF(k-1))
292           ENDDO           ENDDO
293          ENDDO          ENDDO
294         ENDDO         ENDDO
295         DO J=jMin,jMax         DO j=jMin,jMax
296          DO I=iMin,iMax          DO i=iMin,iMax
297            GGL90mixingLength(I,J,Nr) = MIN(GGL90mixingLength(I,J,Nr),            GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
298       &       GGL90mixingLengthMin+drF(Nr))       &       GGL90mixingLengthMin+drF(Nr))
299          ENDDO          ENDDO
300         ENDDO         ENDDO
301         DO k=Nr-1,2,-1         DO k=Nr-1,2,-1
302          DO J=jMin,jMax          DO j=jMin,jMax
303           DO I=iMin,iMax           DO i=iMin,iMax
304            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
305       &        GGL90mixingLength(I,J,K+1)+drF(k))       &        GGL90mixingLength(i,j,k+1)+drF(k))
306           ENDDO           ENDDO
307          ENDDO          ENDDO
308         ENDDO         ENDDO
309    
310         DO k=2,Nr         DO k=2,Nr
311          DO J=jMin,jMax          DO j=jMin,jMax
312           DO I=iMin,iMax           DO i=iMin,iMax
313            GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
314       &                                   GGL90mixingLengthMin)       &                                   GGL90mixingLengthMin)
315            rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K)            rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
316           ENDDO           ENDDO
317          ENDDO          ENDDO
318         ENDDO         ENDDO
319    
320        ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN
321  C-  
322         DO k=2,Nr         DO k=2,Nr
323          DO J=jMin,jMax          DO j=jMin,jMax
324           DO I=iMin,iMax           DO i=iMin,iMax
325            mxLength_Dn(I,J,K) = MIN(GGL90mixingLength(I,J,K),            mxLength_Dn(i,j,k) = MIN(GGL90mixingLength(i,j,k),
326       &        mxLength_Dn(I,J,K-1)+drF(k-1))       &        mxLength_Dn(i,j,k-1)+drF(k-1))
327           ENDDO           ENDDO
328          ENDDO          ENDDO
329         ENDDO         ENDDO
330         DO J=jMin,jMax         DO j=jMin,jMax
331          DO I=iMin,iMax          DO i=iMin,iMax
332            GGL90mixingLength(I,J,Nr) = MIN(GGL90mixingLength(I,J,Nr),            GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
333       &       GGL90mixingLengthMin+drF(Nr))       &       GGL90mixingLengthMin+drF(Nr))
334          ENDDO          ENDDO
335         ENDDO         ENDDO
336         DO k=Nr-1,2,-1         DO k=Nr-1,2,-1
337          DO J=jMin,jMax          DO j=jMin,jMax
338           DO I=iMin,iMax           DO i=iMin,iMax
339            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
340       &        GGL90mixingLength(I,J,K+1)+drF(k))       &        GGL90mixingLength(i,j,k+1)+drF(k))
341           ENDDO           ENDDO
342          ENDDO          ENDDO
343         ENDDO         ENDDO
344    
345         DO k=2,Nr         DO k=2,Nr
346          DO J=jMin,jMax          DO j=jMin,jMax
347           DO I=iMin,iMax           DO i=iMin,iMax
348            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
349       &                                  mxLength_Dn(I,J,K))       &                                  mxLength_Dn(i,j,k))
350            tmpmlx = SQRT( GGL90mixingLength(I,J,K)*mxLength_Dn(I,J,K) )            tmpmlx = SQRT( GGL90mixingLength(i,j,k)*mxLength_Dn(i,j,k) )
351            tmpmlx = MAX( tmpmlx, GGL90mixingLengthMin)            tmpmlx = MAX( tmpmlx, GGL90mixingLengthMin)
352            rMixingLength(I,J,K) = 1. _d 0 / tmpmlx            rMixingLength(i,j,k) = 1. _d 0 / tmpmlx
353           ENDDO           ENDDO
354          ENDDO          ENDDO
355         ENDDO         ENDDO
# Line 323  C- Line 358  C-
358         STOP 'GGL90_CALC: Wrong mxlMaxFlag (mixing length limit)'         STOP 'GGL90_CALC: Wrong mxlMaxFlag (mixing length limit)'
359        ENDIF        ENDIF
360    
361  C- Impose minimum mixing length (to avoid division by zero)  C     start "proper" k-loop (the code above was moved out and up to
362  c      DO k=2,Nr  C     implemement various mixing parameters efficiently)
 c      DO J=jMin,jMax  
 c       DO I=iMin,iMax  
 c        GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),  
 c    &        GGL90mixingLengthMin)  
 c        rMixingLength(I,J,K) = 1. _d 0 /GGL90mixingLength(I,J,K)  
 c       ENDDO  
 c      ENDDO  
 c     ENDDO  
   
   
363        DO k=2,Nr        DO k=2,Nr
364         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)  
          GGL90visctmp(I,J,K) = MAX(KappaM,diffKrNrT(k))  
      &                            * maskC(I,J,K,bi,bj)  
 c        note: storing GGL90visctmp like this, and using it later to compute  
 c              GGL9rdiffKr etc. is robust in case of smoothing (e.g. see OPA)  
          KappaM = MAX(KappaM,viscArNr(k)) * maskC(I,J,K,bi,bj)  
          KappaH = KappaM/TKEPrandtlNumber(I,J,K)  
          KappaE(I,J,K) = GGL90alpha * KappaM * maskC(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  
365    
 C     horizontal diffusion of TKE (requires an exchange in  
 C      do_fields_blocking_exchanges)  
366  #ifdef ALLOW_GGL90_HORIZDIFF  #ifdef ALLOW_GGL90_HORIZDIFF
367        IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN        IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
368         DO K=2,Nr  C     horizontal diffusion of TKE (requires an exchange in
369    C      do_fields_blocking_exchanges)
370  C     common factors  C     common factors
371          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
372           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
373            xA(i,j) = _dyG(i,j,bi,bj)            xA(i,j) = _dyG(i,j,bi,bj)*drC(k)*
374       &         *drF(k)*_hFacW(i,j,k,bi,bj)       &                 (min(.5 _d 0,_hFacW(i,j,k-1,bi,bj) ) +
375            yA(i,j) = _dxG(i,j,bi,bj)       &                  min(.5 _d 0,_hFacW(i,j,k  ,bi,bj) ) )
376       &         *drF(k)*_hFacS(i,j,k,bi,bj)            yA(i,j) = _dxG(i,j,bi,bj)*drC(k)*
377         &                 (min(.5 _d 0,_hFacS(i,j,k-1,bi,bj) ) +
378         &                  min(.5 _d 0,_hFacS(i,j,k  ,bi,bj) ) )
379           ENDDO           ENDDO
380          ENDDO          ENDDO
381  C     Compute diffusive fluxes  C     Compute diffusive fluxes
382  C     ... across x-faces  C     ... across x-faces
383          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
384           dfx(1-Olx,j)=0. _d 0           dfx(1-OLx,j)=0. _d 0
385           DO i=1-Olx+1,sNx+Olx           DO i=1-OLx+1,sNx+OLx
386            dfx(i,j) = -GGL90diffTKEh*xA(i,j)            dfx(i,j) = -GGL90diffTKEh*xA(i,j)
387       &      *_recip_dxC(i,j,bi,bj)       &      *_recip_dxC(i,j,bi,bj)
388       &      *(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))
389    #ifdef ISOTROPIC_COS_SCALING
390       &      *CosFacU(j,bi,bj)       &      *CosFacU(j,bi,bj)
391    #endif /* ISOTROPIC_COS_SCALING */
392           ENDDO           ENDDO
393          ENDDO          ENDDO
394  C     ... across y-faces  C     ... across y-faces
395          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
396           dfy(i,1-Oly)=0. _d 0           dfy(i,1-OLy)=0. _d 0
397          ENDDO          ENDDO
398          DO j=1-Oly+1,sNy+Oly          DO j=1-OLy+1,sNy+OLy
399           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
400            dfy(i,j) = -GGL90diffTKEh*yA(i,j)            dfy(i,j) = -GGL90diffTKEh*yA(i,j)
401       &      *_recip_dyC(i,j,bi,bj)       &      *_recip_dyC(i,j,bi,bj)
402       &      *(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 419  C     ... across y-faces Line 406  C     ... across y-faces
406           ENDDO           ENDDO
407          ENDDO          ENDDO
408  C     Compute divergence of fluxes  C     Compute divergence of fluxes
409          DO j=1-Oly,sNy+Oly-1          DO j=1-OLy,sNy+OLy-1
410           DO i=1-Olx,sNx+Olx-1           DO i=1-OLx,sNx+OLx-1
411            gTKE(i,j,k)=gTKE(i,j,k)            gTKE(i,j) = -recip_drC(k)*recip_rA(i,j,bi,bj)
412       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj)       &         *recip_hFacI(i,j,k)
413       &         *( (dfx(i+1,j)-dfx(i,j))       &         *((dfx(i+1,j)-dfx(i,j))
414       &           +(dfy(i,j+1)-dfy(i,j))       &         + (dfy(i,j+1)-dfy(i,j)) )
      &           )*deltaTggl90  
415           ENDDO           ENDDO
416          ENDDO          ENDDO
 C       end of k-loop  
        ENDDO  
417  C     end if GGL90diffTKEh .eq. 0.  C     end if GGL90diffTKEh .eq. 0.
418        ENDIF         ENDIF
419  #endif /* ALLOW_GGL90_HORIZDIFF */  #endif /* ALLOW_GGL90_HORIZDIFF */
420    
421    C     viscosity and diffusivity
422           DO j=jMin,jMax
423            DO i=iMin,iMax
424             KappaM(i,j) = GGL90ck*GGL90mixingLength(i,j,k)*SQRTTKE(i,j,k)
425             GGL90visctmp(i,j,k) = MAX(KappaM(i,j),diffKrNrS(k))
426         &                            * maskC(i,j,k,bi,bj)
427    C        note: storing GGL90visctmp like this, and using it later to compute
428    C              GGL9rdiffKr etc. is robust in case of smoothing (e.g. see OPA)
429             KappaM(i,j) = MAX(KappaM(i,j),viscArNr(k)) * maskC(i,j,k,bi,bj)
430            ENDDO
431           ENDDO
432    
433    C     compute vertical shear (dU/dz)^2+(dV/dz)^2
434           IF ( calcMeanVertShear ) THEN
435    C     by averaging (@ grid-cell center) the 4 vertical shear compon @ U,V pos.
436            DO j=jMin,jMax
437             DO i=iMin,iMax
438              tempU  = ( uVel( i ,j,km1,bi,bj) - uVel( i ,j,k,bi,bj) )
439              tempUp = ( uVel(i+1,j,km1,bi,bj) - uVel(i+1,j,k,bi,bj) )
440              tempV  = ( vVel(i, j ,km1,bi,bj) - vVel(i, j ,k,bi,bj) )
441              tempVp = ( vVel(i,j+1,km1,bi,bj) - vVel(i,j+1,k,bi,bj) )
442              verticalShear(i,j) = (
443         &                 ( tempU*tempU + tempUp*tempUp )*halfRL
444         &               + ( tempV*tempV + tempVp*tempVp )*halfRL
445         &                         )*recip_drC(k)*recip_drC(k)
446             ENDDO
447            ENDDO
448           ELSE
449    C     from the averaged flow at grid-cell center (2 compon x 2 pos.)
450            DO j=jMin,jMax
451             DO i=iMin,iMax
452              tempU = ( ( uVel(i,j,km1,bi,bj) + uVel(i+1,j,km1,bi,bj) )
453         &             -( uVel(i,j,k  ,bi,bj) + uVel(i+1,j,k  ,bi,bj) )
454         &            )*halfRL*recip_drC(k)
455              tempV = ( ( vVel(i,j,km1,bi,bj) + vVel(i,j+1,km1,bi,bj) )
456         &             -( vVel(i,j,k  ,bi,bj) + vVel(i,j+1,k  ,bi,bj) )
457         &            )*halfRL*recip_drC(k)
458              verticalShear(i,j) = tempU*tempU + tempV*tempV
459             ENDDO
460            ENDDO
461           ENDIF
462    
463    C     compute Prandtl number (always greater than 0)
464    #ifdef ALLOW_GGL90_IDEMIX
465           IF ( useIDEMIX ) THEN
466            DO j=jMin,jMax
467             DO i=iMin,iMax
468    C     account for partical cell factor in vertical shear:
469              verticalShear(i,j) = verticalShear(i,j)
470         &                       * recip_hFacI(i,j,k)*recip_hFacI(i,j,k)
471              RiNumber = MAX(Nsquare(i,j,k),0. _d 0)
472         &         /(verticalShear(i,j)+GGL90eps)
473    CML         IDEMIX_RiNumber = 1./GGL90eps
474              IDEMIX_RiNumber = MAX( KappaM(i,j)*Nsquare(i,j,k), 0. _d 0)/
475         &     (GGL90eps+IDEMIX_tau_d(i,j,k,bi,bj)*IDEMIX_E(i,j,k,bi,bj)**2)
476              prTemp         = MIN(5.*RiNumber, 6.6 _d 0*IDEMIX_RiNumber)
477              TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
478              TKEPrandtlNumber(i,j,k) = MAX( 1. _d 0,TKEPrandtlNumber(i,j,k))
479             ENDDO
480            ENDDO
481           ELSE
482    #else /* ndef ALLOW_GGL90_IDEMIX */
483           IF (.TRUE.) THEN
484    #endif /* ALLOW_GGL90_IDEMIX */
485            DO j=jMin,jMax
486             DO i=iMin,iMax
487              RiNumber = MAX(Nsquare(i,j,k),0. _d 0)
488         &         /(verticalShear(i,j)+GGL90eps)
489              prTemp = 1. _d 0
490              IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
491              TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
492             ENDDO
493            ENDDO
494           ENDIF
495    
496           DO j=jMin,jMax
497            DO i=iMin,iMax
498    C        diffusivity
499             KappaH = KappaM(i,j)/TKEPrandtlNumber(i,j,k)
500             KappaE(i,j,k) = GGL90alpha * KappaM(i,j) * maskC(i,j,k,bi,bj)
501    
502    C     dissipation term
503             TKEdissipation = explDissFac*GGL90ceps
504         &        *SQRTTKE(i,j,k)*rMixingLength(i,j,k)
505         &        *GGL90TKE(i,j,k,bi,bj)
506    C     partial update with sum of explicit contributions
507             GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
508         &        + deltaTggl90*(
509         &        + KappaM(i,j)*verticalShear(i,j)
510         &        - KappaH*Nsquare(i,j,k)
511         &        - TKEdissipation
512         &        )
513            ENDDO
514           ENDDO
515    
516    #ifdef ALLOW_GGL90_IDEMIX
517           IF ( useIDEMIX ) THEN
518    C     add IDEMIX contribution to the turbulent kinetic energy
519            DO j=jMin,jMax
520             DO i=iMin,iMax
521              GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
522         &         + deltaTggl90*(
523         &         + IDEMIX_tau_d(i,j,k,bi,bj)*IDEMIX_E(i,j,k,bi,bj)**2
524         &         )
525             ENDDO
526            ENDDO
527           ENDIF
528    #endif /* ALLOW_GGL90_IDEMIX */
529    
530    #ifdef ALLOW_GGL90_HORIZDIFF
531           IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
532    C--    Add horiz. diffusion tendency
533            DO j=jMin,jMax
534             DO i=iMin,iMax
535              GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
536         &                          + gTKE(i,j)*deltaTggl90
537             ENDDO
538            ENDDO
539           ENDIF
540    #endif /* ALLOW_GGL90_HORIZDIFF */
541    
542    C--   end of k loop
543          ENDDO
544    
545  C     ============================================  C     ============================================
546  C     Implicit time step to update TKE for k=1,Nr;  C     Implicit time step to update TKE for k=1,Nr;
547  C     TKE(Nr+1)=0 by default  C     TKE(Nr+1)=0 by default
# Line 442  C     set up matrix Line 550  C     set up matrix
550  C--   Lower diagonal  C--   Lower diagonal
551        DO j=jMin,jMax        DO j=jMin,jMax
552         DO i=iMin,iMax         DO i=iMin,iMax
553           a(i,j,1) = 0. _d 0           a3d(i,j,1) = 0. _d 0
554         ENDDO         ENDDO
555        ENDDO        ENDDO
556        DO k=2,Nr        DO k=2,Nr
# Line 452  C--   Lower diagonal Line 560  C--   Lower diagonal
560  C-    We keep recip_hFacC in the diffusive flux calculation,  C-    We keep recip_hFacC in the diffusive flux calculation,
561  C-    but no hFacC in TKE volume control  C-    but no hFacC in TKE volume control
562  C-    No need for maskC(k-1) with recip_hFacC(k-1)  C-    No need for maskC(k-1) with recip_hFacC(k-1)
563           a(i,j,k) = -deltaTggl90           a3d(i,j,k) = -deltaTggl90
564       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)
565       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))
566       &        *recip_drC(k)*maskC(i,j,k,bi,bj)       &        *recip_drC(k)*maskC(i,j,k,bi,bj)
# Line 462  C-    No need for maskC(k-1) with recip_ Line 570  C-    No need for maskC(k-1) with recip_
570  C--   Upper diagonal  C--   Upper diagonal
571        DO j=jMin,jMax        DO j=jMin,jMax
572         DO i=iMin,iMax         DO i=iMin,iMax
573           c(i,j,1)  = 0. _d 0           c3d(i,j,1)  = 0. _d 0
574         ENDDO         ENDDO
575        ENDDO        ENDDO
576        DO k=2,Nr        DO k=2,Nr
577         DO j=jMin,jMax         DO j=jMin,jMax
578          DO i=iMin,iMax          DO i=iMin,iMax
579            kp1=MAX(1,MIN(klowC(i,j,bi,bj),k+1))           kp1=MAX(1,MIN(klowC(i,j,bi,bj),k+1))
580  C-    We keep recip_hFacC in the diffusive flux calculation,  C-    We keep recip_hFacC in the diffusive flux calculation,
581  C-    but no hFacC in TKE volume control  C-    but no hFacC in TKE volume control
582  C-    No need for maskC(k) with recip_hFacC(k)  C-    No need for maskC(k) with recip_hFacC(k)
583            c(i,j,k) = -deltaTggl90           c3d(i,j,k) = -deltaTggl90
584       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)
585       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))
586       &        *recip_drC(k)*maskC(i,j,k-1,bi,bj)       &        *recip_drC(k)*maskC(i,j,k-1,bi,bj)
587          ENDDO          ENDDO
588         ENDDO         ENDDO
589        ENDDO        ENDDO
590    
591    #ifdef ALLOW_GGL90_IDEMIX
592          IF ( useIDEMIX ) THEN
593           DO k=2,Nr
594            DO j=jMin,jMax
595             DO i=iMin,iMax
596              a3d(i,j,k) = a3d(i,j,k)*recip_hFacI(i,j,k)
597              c3d(i,j,k) = c3d(i,j,k)*recip_hFacI(i,j,k)
598             ENDDO
599            ENDDO
600           ENDDO
601          ENDIF
602    #endif /* ALLOW_GGL90_IDEMIX */
603    
604          IF (.NOT.GGL90_dirichlet) THEN
605    C      Neumann bottom boundary condition for TKE: no flux from bottom
606           DO j=jMin,jMax
607            DO i=iMin,iMax
608             kBottom   = MAX(kLowC(i,j,bi,bj),1)
609             c3d(i,j,kBottom) = 0. _d 0
610            ENDDO
611           ENDDO
612          ENDIF
613    
614  C--   Center diagonal  C--   Center diagonal
615        DO k=1,Nr        DO k=1,Nr
616         km1 = MAX(k-1,1)         km1 = MAX(k-1,1)
617         DO j=jMin,jMax         DO j=jMin,jMax
618          DO i=iMin,iMax          DO i=iMin,iMax
619            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)
620       &        + ab15*deltaTggl90*GGL90ceps*SQRTTKE(I,J,K)       &        + implDissFac*deltaTggl90*GGL90ceps*SQRTTKE(i,j,k)
621       &        * rMixingLength(I,J,K)       &        * rMixingLength(i,j,k)
622       &        * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)       &        * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)
623           ENDDO          ENDDO
624         ENDDO         ENDDO
625        ENDDO        ENDDO
626  C     end set up matrix  C     end set up matrix
627    
628  C     Apply boundary condition  C     Apply boundary condition
629        kp1 = MIN(Nr,kSurf+1)        kp1 = MIN(Nr,kSurf+1)
630        DO J=jMin,jMax        DO j=jMin,jMax
631         DO I=iMin,iMax         DO i=iMin,iMax
632  C     estimate friction velocity uStar from surface forcing  C     estimate friction velocity uStar from surface forcing
633          uStarSquare = SQRT(          uStarSquare = SQRT(
634       &    ( .5 _d 0*( surfaceForcingU(I,  J,  bi,bj)       &    ( .5 _d 0*( surfaceForcingU(i,  j,  bi,bj)
635       &              + surfaceForcingU(I+1,J,  bi,bj) ) )**2       &              + surfaceForcingU(i+1,j,  bi,bj) ) )**2
636       &  + ( .5 _d 0*( surfaceForcingV(I,  J,  bi,bj)       &  + ( .5 _d 0*( surfaceForcingV(i,  j,  bi,bj)
637       &              + surfaceForcingV(I,  J+1,bi,bj) ) )**2       &              + surfaceForcingV(i,  j+1,bi,bj) ) )**2
638       &                     )       &                     )
639  C     Dirichlet surface boundary condition for TKE  C     Dirichlet surface boundary condition for TKE
640          gTKE(I,J,kSurf) = MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)          GGL90TKE(i,j,kSurf,bi,bj) = maskC(i,j,kSurf,bi,bj)
641       &                     *maskC(I,J,kSurf,bi,bj)       &           *MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)
642          gTKE(i,j,kp1) = gTKE(i,j,kp1)          GGL90TKE(i,j,kp1,bi,bj) = GGL90TKE(i,j,kp1,bi,bj)
643       &                - a(i,j,kp1)*gTKE(i,j,kSurf)       &               - a3d(i,j,kp1)*GGL90TKE(i,j,kSurf,bi,bj)
644          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  
645         ENDDO         ENDDO
646        ENDDO        ENDDO
647    
648  C     solve tri-diagonal system, and store solution on gTKE (previously rhs)        IF (GGL90_dirichlet) THEN
649    C      Dirichlet bottom boundary condition for TKE = GGL90TKEbottom
650           DO j=jMin,jMax
651            DO i=iMin,iMax
652             kBottom   = MAX(kLowC(i,j,bi,bj),1)
653             GGL90TKE(i,j,kBottom,bi,bj) = GGL90TKE(i,j,kBottom,bi,bj)
654         &                              - GGL90TKEbottom*c3d(i,j,kBottom)
655             c3d(i,j,kBottom) = 0. _d 0
656            ENDDO
657           ENDDO
658          ENDIF
659    
660    C     solve tri-diagonal system
661        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
662       I                        a, b, c,       I                        a3d, b3d, c3d,
663       U                        gTKE,       U                        GGL90TKE(1-OLx,1-OLy,1,bi,bj),
664       O                        errCode,       O                        errCode,
665       I                        1, 1, myThid )       I                        bi, bj, myThid )
666    
667  C     now update TKE        DO k=1,Nr
668        DO K=1,Nr         DO j=jMin,jMax
669         DO J=jMin,jMax          DO i=iMin,iMax
         DO I=iMin,iMax  
670  C     impose minimum TKE to avoid numerical undershoots below zero  C     impose minimum TKE to avoid numerical undershoots below zero
671           GGL90TKE(I,J,K,bi,bj) = MAX( gTKE(I,J,K), GGL90TKEmin )           GGL90TKE(i,j,k,bi,bj) = maskC(i,j,k,bi,bj)
672       &        * maskC(I,J,K,bi,bj)       &                  *MAX( GGL90TKE(i,j,k,bi,bj), GGL90TKEmin )
673          ENDDO          ENDDO
674         ENDDO         ENDDO
675        ENDDO        ENDDO
# Line 539  C     impose minimum TKE to avoid numeri Line 677  C     impose minimum TKE to avoid numeri
677  C     end of time step  C     end of time step
678  C     ===============================  C     ===============================
679    
680        DO K=2,Nr        DO k=2,Nr
681         DO J=1,sNy         DO j=1,sNy
682          DO I=1,sNx          DO i=1,sNx
683  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
684           tmpVisc=           tmpVisc = (
685       &  (       &     p4 *    GGL90visctmp(i  ,j  ,k)*mskCor(i  ,j  ,bi,bj)
686       &   p4 *  GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)       &    +p8 *( ( GGL90visctmp(i-1,j  ,k)*mskCor(i-1,j  ,bi,bj)
687       &  +p8 *( GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj)       &           + GGL90visctmp(i+1,j  ,k)*mskCor(i+1,j  ,bi,bj) )
688       &       + GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj)       &         + ( GGL90visctmp(i  ,j-1,k)*mskCor(i  ,j-1,bi,bj)
689       &       + GGL90visctmp(i+1,j  ,k) * mskCor(i+1,j  ,bi,bj)       &           + GGL90visctmp(i  ,j+1,k)*mskCor(i  ,j+1,bi,bj) ) )
690       &       + GGL90visctmp(i  ,j+1,k) * mskCor(i  ,j+1,bi,bj))       &    +p16*( ( GGL90visctmp(i+1,j+1,k)*mskCor(i+1,j+1,bi,bj)
691       &  +p16*( GGL90visctmp(i+1,j+1,k) * mskCor(i+1,j+1,bi,bj)       &           + GGL90visctmp(i-1,j-1,k)*mskCor(i-1,j-1,bi,bj) )
692       &       + GGL90visctmp(i+1,j-1,k) * mskCor(i+1,j-1,bi,bj)       &         + ( GGL90visctmp(i+1,j-1,k)*mskCor(i+1,j-1,bi,bj)
693       &       + GGL90visctmp(i-1,j+1,k) * mskCor(i-1,j+1,bi,bj)       &           + GGL90visctmp(i-1,j+1,k)*mskCor(i-1,j+1,bi,bj) ) )
694       &       + GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj))       &             )/(
695       &  )       &     p4
696       & /(p4       &    +p8 *( (  maskC(i-1,j  ,k,bi,bj)*mskCor(i-1,j  ,bi,bj)
697       &  +p8 *(       maskC(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)       &           +  maskC(i+1,j  ,k,bi,bj)*mskCor(i+1,j  ,bi,bj) )
698       &       +       maskC(i  ,j-1,k,bi,bj) * mskCor(i  ,j-1,bi,bj)       &         + (  maskC(i  ,j-1,k,bi,bj)*mskCor(i  ,j-1,bi,bj)
699       &       +       maskC(i+1,j  ,k,bi,bj) * mskCor(i+1,j  ,bi,bj)       &           +  maskC(i  ,j+1,k,bi,bj)*mskCor(i  ,j+1,bi,bj) ) )
700       &       +       maskC(i  ,j+1,k,bi,bj) * mskCor(i  ,j+1,bi,bj))       &    +p16*( (  maskC(i+1,j+1,k,bi,bj)* mskCor(i+1,j+1,bi,bj)
701       &  +p16*(       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) )
702       &       +       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)
703       &       +       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) ) )
704       &       +       maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj))       &               )*maskC(i,j,k,bi,bj)*mskCor(i,j,bi,bj)
      &  )*maskC(i,j,k,bi,bj)*mskCor(i,j,bi,bj)  
705  #else  #else
706           tmpVisc = GGL90visctmp(I,J,K)           tmpVisc = GGL90visctmp(i,j,k)
707  #endif  #endif
708           tmpVisc = MIN(tmpVisc/TKEPrandtlNumber(i,j,k),GGL90diffMax)           tmpVisc = MIN(tmpVisc/TKEPrandtlNumber(i,j,k),GGL90diffMax)
709           GGL90diffKr(I,J,K,bi,bj)= MAX( tmpVisc , diffKrNrT(k) )           GGL90diffKr(i,j,k,bi,bj)= MAX( tmpVisc , diffKrNrS(k) )
710          ENDDO          ENDDO
711         ENDDO         ENDDO
712        ENDDO        ENDDO
713    
714          DO k=2,Nr
715           DO j=1,sNy
716        DO K=2,Nr          DO i=1,sNx+1
        DO J=1,sNy  
         DO I=1,sNx+1  
717  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
718          tmpVisc =           tmpVisc = (
719       & (       &     p4 *(   GGL90visctmp(i-1,j  ,k)*mskCor(i-1,j  ,bi,bj)
720       &   p4 *(GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)       &           + GGL90visctmp(i  ,j  ,k)*mskCor(i  ,j  ,bi,bj) )
721       &       +GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj))       &    +p8 *( ( GGL90visctmp(i-1,j-1,k)*mskCor(i-1,j-1,bi,bj)
722       &  +p8 *(GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj)       &           + GGL90visctmp(i  ,j-1,k)*mskCor(i  ,j-1,bi,bj) )
723       &       +GGL90visctmp(i-1,j+1,k) * mskCor(i-1,j+1,bi,bj)       &         + ( GGL90visctmp(i-1,j+1,k)*mskCor(i-1,j+1,bi,bj)
724       &       +GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj)       &           + GGL90visctmp(i  ,j+1,k)*mskCor(i  ,j+1,bi,bj) ) )
725       &       +GGL90visctmp(i  ,j+1,k) * mskCor(i  ,j+1,bi,bj))       &             )/(
726       &  )       &     p4 * 2. _d 0
727       & /(p4 * 2. _d 0       &    +p8 *( (  maskC(i-1,j-1,k,bi,bj)*mskCor(i-1,j-1,bi,bj)
728       &  +p8 *(      maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj)       &           +  maskC(i  ,j-1,k,bi,bj)*mskCor(i  ,j-1,bi,bj) )
729       &       +      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)
730       &       +      maskC(i  ,j-1,k,bi,bj) * mskCor(i  ,j-1,bi,bj)       &           +  maskC(i  ,j+1,k,bi,bj)*mskCor(i  ,j+1,bi,bj) ) )
731       &       +      maskC(i  ,j+1,k,bi,bj) * mskCor(i  ,j+1,bi,bj))       &               )*maskC(i-1,j,k,bi,bj)*mskCor(i-1,j,bi,bj)
732       &  )       &                *maskC(i  ,j,k,bi,bj)*mskCor(i  ,j,bi,bj)
      &  *maskC(i  ,j,k,bi,bj)*mskCor(i  ,j,bi,bj)  
      &  *maskC(i-1,j,k,bi,bj)*mskCor(i-1,j,bi,bj)  
733  #else  #else
734          tmpVisc = _maskW(i,j,k,bi,bj) *           tmpVisc = _maskW(i,j,k,bi,bj) * halfRL
735       &                   (.5 _d 0*(GGL90visctmp(i,j,k)       &          *( GGL90visctmp(i-1,j,k)
736       &                            +GGL90visctmp(i-1,j,k))       &           + GGL90visctmp(i,j,k) )
      &                   )  
737  #endif  #endif
738          tmpVisc = MIN( tmpVisc , GGL90viscMax )           tmpVisc = MIN( tmpVisc , GGL90viscMax )
739          GGL90viscArU(i,j,k,bi,bj) = MAX( tmpVisc , viscArNr(k)  )           GGL90viscArU(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
740          ENDDO          ENDDO
741         ENDDO         ENDDO
742        ENDDO        ENDDO
743    
744          DO k=2,Nr
745        DO K=2,Nr         DO j=1,sNy+1
746         DO J=1,sNy+1          DO i=1,sNx
         DO I=1,sNx  
747  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
748          tmpVisc =           tmpVisc = (
749       & (       &     p4 *(   GGL90visctmp(i  ,j-1,k)*mskCor(i  ,j-1,bi,bj)
750       &   p4 *(GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)       &           + GGL90visctmp(i  ,j  ,k)*mskCor(i  ,j  ,bi,bj) )
751       &       +GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj))       &    +p8 *( ( GGL90visctmp(i-1,j-1,k)*mskCor(i-1,j-1,bi,bj)
752       &  +p8 *(GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj)       &           + GGL90visctmp(i-1,j  ,k)*mskCor(i-1,j  ,bi,bj) )
753       &       +GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj)       &         + ( GGL90visctmp(i+1,j-1,k)*mskCor(i+1,j-1,bi,bj)
754       &       +GGL90visctmp(i+1,j  ,k) * mskCor(i+1,j  ,bi,bj)       &           + GGL90visctmp(i+1,j  ,k)*mskCor(i+1,j  ,bi,bj) ) )
755       &       +GGL90visctmp(i+1,j-1,k) * mskCor(i+1,j-1,bi,bj))       &             )/(
756       &  )       &     p4 * 2. _d 0
757       & /(p4 * 2. _d 0       &    +p8 *( (  maskC(i-1,j-1,k,bi,bj)*mskCor(i-1,j-1,bi,bj)
758       &  +p8 *(      maskC(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)       &           +  maskC(i-1,j  ,k,bi,bj)*mskCor(i-1,j  ,bi,bj) )
759       &       +      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)
760       &       +      maskC(i+1,j  ,k,bi,bj) * mskCor(i+1,j  ,bi,bj)       &           +  maskC(i+1,j  ,k,bi,bj)*mskCor(i+1,j  ,bi,bj) ) )
761       &       +      maskC(i+1,j-1,k,bi,bj) * mskCor(i+1,j-1,bi,bj))       &               )*maskC(i,j-1,k,bi,bj)*mskCor(i,j-1,bi,bj)
762       &  )       &                *maskC(i,j  ,k,bi,bj)*mskCor(i,j  ,bi,bj)
      &   *maskC(i,j  ,k,bi,bj)*mskCor(i,j  ,bi,bj)  
      &   *maskC(i,j-1,k,bi,bj)*mskCor(i,j-1,bi,bj)  
763  #else  #else
764          tmpVisc = _maskS(i,j,k,bi,bj) *           tmpVisc = _maskS(i,j,k,bi,bj) * halfRL
765       &                   (.5 _d 0*(GGL90visctmp(i,j,k)       &          *( GGL90visctmp(i,j-1,k)
766       &                            +GGL90visctmp(i,j-1,k))       &           + GGL90visctmp(i,j,k) )
      &                   )  
   
767  #endif  #endif
768          tmpVisc = MIN( tmpVisc , GGL90viscMax )           tmpVisc = MIN( tmpVisc , GGL90viscMax )
769          GGL90viscArV(i,j,k,bi,bj) = MAX( tmpVisc , viscArNr(k)  )           GGL90viscArV(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
770          ENDDO          ENDDO
771         ENDDO         ENDDO
772        ENDDO        ENDDO
773    
774  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
775        IF ( useDiagnostics ) THEN        IF ( useDiagnostics ) THEN
776           CALL DIAGNOSTICS_FILL( GGL90TKE   ,'GGL90TKE',          CALL DIAGNOSTICS_FILL( GGL90TKE   ,'GGL90TKE',
777       &                          0,Nr, 1, bi, bj, myThid )       &                         0,Nr, 1, bi, bj, myThid )
778           CALL DIAGNOSTICS_FILL( GGL90viscArU,'GGL90ArU',          CALL DIAGNOSTICS_FILL( GGL90viscArU,'GGL90ArU',
779       &                          0,Nr, 1, bi, bj, myThid )       &                         0,Nr, 1, bi, bj, myThid )
780           CALL DIAGNOSTICS_FILL( GGL90viscArV,'GGL90ArV',          CALL DIAGNOSTICS_FILL( GGL90viscArV,'GGL90ArV',
781       &                          0,Nr, 1, bi, bj, myThid )       &                         0,Nr, 1, bi, bj, myThid )
782           CALL DIAGNOSTICS_FILL( GGL90diffKr,'GGL90Kr ',          CALL DIAGNOSTICS_FILL( GGL90diffKr,'GGL90Kr ',
783       &                          0,Nr, 1, bi, bj, myThid )       &                         0,Nr, 1, bi, bj, myThid )
784           CALL DIAGNOSTICS_FILL( TKEPrandtlNumber ,'GGL90Prl',          CALL DIAGNOSTICS_FILL( TKEPrandtlNumber ,'GGL90Prl',
785       &                          0,Nr, 2, bi, bj, myThid )       &                         0,Nr, 2, bi, bj, myThid )
786           CALL DIAGNOSTICS_FILL( GGL90mixingLength,'GGL90Lmx',          CALL DIAGNOSTICS_FILL( GGL90mixingLength,'GGL90Lmx',
787       &                          0,Nr, 2, bi, bj, myThid )       &                         0,Nr, 2, bi, bj, myThid )
788    
789            kp1 = MIN(Nr,kSurf+1)
790            DO j=jMin,jMax
791             DO i=iMin,iMax
792    C     diagnose surface flux of TKE
793              surf_flx_tke(i,j) =(GGL90TKE(i,j,kSurf,bi,bj)-
794         &                        GGL90TKE(i,j,kp1,bi,bj))
795         &        *recip_drF(kSurf)*recip_hFacC(i,j,kSurf,bi,bj)
796         &        *KappaE(i,j,kp1)
797             ENDDO
798            ENDDO
799            CALL DIAGNOSTICS_FILL( surf_flx_tke,'GGL90flx',
800         &                         0, 1, 2, bi, bj, myThid )
801    
802            k=kSurf
803            DO j=jMin,jMax
804             DO i=iMin,iMax
805    C     diagnose work done by the wind
806              surf_flx_tke(i,j) =
807         &      halfRL*( surfaceForcingU(i,  j,bi,bj)*uVel(i  ,j,k,bi,bj)
808         &              +surfaceForcingU(i+1,j,bi,bj)*uVel(i+1,j,k,bi,bj))
809         &    + halfRL*( surfaceForcingV(i,j,  bi,bj)*vVel(i,j  ,k,bi,bj)
810         &              +surfaceForcingV(i,j+1,bi,bj)*vVel(i,j+1,k,bi,bj))
811             ENDDO
812            ENDDO
813            CALL DIAGNOSTICS_FILL( surf_flx_tke,'GGL90tau',
814         &                         0, 1, 2, bi, bj, myThid )
815    
816        ENDIF        ENDIF
817  #endif  #endif /* ALLOW_DIAGNOSTICS */
818    
819  #endif /* ALLOW_GGL90 */  #endif /* ALLOW_GGL90 */
820    
Legend:
Removed from v.1.18  
changed lines
  Added in v.1.34
  ViewVC Help
Powered by ViewVC 1.1.22