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