/[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.8 by jmc, Mon Aug 11 22:28:06 2008 UTC revision 1.26 by jmc, Thu Aug 14 16:42:35 2014 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    CEOP
55    
56  #ifdef ALLOW_GGL90  #ifdef ALLOW_GGL90
57    
58  C !LOCAL VARIABLES: ====================================================  C !LOCAL VARIABLES: ====================================================
59  c Local constants  C Local constants
60  C     iMin, iMax, jMin, jMax, I, J - array computation indices  C     iMin,iMax,jMin,jMax :: index boundaries of computation domain
61  C     K, Kp1, km1, kSurf, kBottom  - vertical loop indices  C     i, j, k, kp1,km1 :: array computation indices
62  C     ab15, ab05       - weights for implicit timestepping  C     kSurf, kBottom   :: vertical indices of domain boundaries
63  C     uStarSquare      - square of friction velocity  C     explDissFac      :: explicit Dissipation Factor (in [0-1])
64  C     verticalShear    - (squared) vertical shear of horizontal velocity  C     implDissFac      :: implicit Dissipation Factor (in [0-1])
65  C     Nsquare          - squared buoyancy freqency  C     uStarSquare      :: square of friction velocity
66  C     RiNumber         - local Richardson number  C     verticalShear    :: (squared) vertical shear of horizontal velocity
67  C     KappaM           - (local) viscosity parameter (eq.10)  C     Nsquare          :: squared buoyancy freqency
68  C     KappaH           - (local) diffusivity parameter for temperature (eq.11)  C     RiNumber         :: local Richardson number
69  C     KappaE           - (local) diffusivity parameter for TKE (eq.15)  C     KappaM           :: (local) viscosity parameter (eq.10)
70  C     buoyFreq         - buoyancy freqency  C     KappaH           :: (local) diffusivity parameter for temperature (eq.11)
71  C     TKEdissipation   - dissipation of TKE  C     KappaE           :: (local) diffusivity parameter for TKE (eq.15)
72  C     GGL90mixingLength- mixing length of scheme following Banke+Delecuse  C     TKEdissipation   :: dissipation of TKE
73  C         rMixingLength- inverse of mixing length  C     GGL90mixingLength:: mixing length of scheme following Banke+Delecuse
74  C     totalDepth       - thickness of water column (inverse of recip_Rcol)  C         rMixingLength:: inverse of mixing length
75  C     TKEPrandtlNumber - here, an empirical function of the Richardson number  C     totalDepth       :: thickness of water column (inverse of recip_Rcol)
76  C     rhoK, rhoKm1     - density at layer K and Km1 (relative to K)  C     TKEPrandtlNumber :: here, an empirical function of the Richardson number
 C     gTKE             - right hand side of implicit equation  
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
82        _RL     KappaM, KappaH        _RL     KappaM, KappaH
83        _RL     Nsquare  c     _RL     Nsquare
84          _RL     Nsquare(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
85        _RL     deltaTggl90        _RL     deltaTggl90
86        _RL     SQRTTKE  c     _RL     SQRTTKE
87          _RL     SQRTTKE(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
88        _RL     RiNumber        _RL     RiNumber
89        _RL     TKEdissipation        _RL     TKEdissipation
90        _RL     tempU, tempV, prTemp        _RL     tempU, tempV, prTemp
91          _RL     MaxLength, tmpmlx, tmpVisc
92        _RL     TKEPrandtlNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     TKEPrandtlNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
93        _RL     GGL90mixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     GGL90mixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
94        _RL         rMixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     rMixingLength    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
95          _RL     mxLength_Dn      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
96        _RL     KappaE           (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     KappaE           (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
       _RL     rhoK             (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL     rhoKm1           (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
97        _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98        _RL     gTKE             (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     GGL90visctmp     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
99  C     tri-diagonal matrix  C-    tri-diagonal matrix
100        _RL     a(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     a3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
101        _RL     b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     b3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
102        _RL     c(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     c3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
103          INTEGER errCode
104  #ifdef ALLOW_GGL90_HORIZDIFF  #ifdef ALLOW_GGL90_HORIZDIFF
105  C     xA, yA   - area of lateral faces  C     hFac     :: fractional thickness of W-cell
106  C     dfx, dfy - diffusive flux across lateral faces  C     xA, yA   :: area of lateral faces
107    C     dfx, dfy :: diffusive flux across lateral faces
108    C     gTKE     :: right hand side of diffusion equation
109          _RL     hFac
110        _RL     xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
111        _RL     yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
112        _RL     dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
113        _RL     dfy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     dfy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
114          _RL    gTKE(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
115  #endif /* ALLOW_GGL90_HORIZDIFF */  #endif /* ALLOW_GGL90_HORIZDIFF */
116  CEOP  #ifdef ALLOW_GGL90_SMOOTH
117          _RL p4, p8, p16
118          p4=0.25 _d 0
119          p8=0.125 _d 0
120          p16=0.0625 _d 0
121    #endif
122        iMin = 2-OLx        iMin = 2-OLx
123        iMax = sNx+OLx-1        iMax = sNx+OLx-1
124        jMin = 2-OLy        jMin = 2-OLy
# Line 113  CEOP Line 126  CEOP
126    
127  C     set separate time step (should be deltaTtracer)  C     set separate time step (should be deltaTtracer)
128        deltaTggl90 = dTtracerLev(1)        deltaTggl90 = dTtracerLev(1)
129  C  
130        kSurf = 1        kSurf = 1
131  C     implicit timestepping weights for dissipation  C     explicit/implicit timestepping weights for dissipation
132        ab15 =  1.5 _d 0        explDissFac = 0. _d 0
133        ab05 = -0.5 _d 0        implDissFac = 1. _d 0 - explDissFac
       ab15 =  1. _d 0  
       ab05 =  0. _d 0  
134    
135  C     Initialize local fields  C     Initialize local fields
136        DO K = 1, Nr        DO k = 1, Nr
137         DO J=1-Oly,sNy+Oly         DO j=1-OLy,sNy+OLy
138          DO I=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
139           gTKE(I,J,K)              = 0. _d 0           KappaE(i,j,k)            = 0. _d 0
140           KappaE(I,J,K)            = 0. _d 0           TKEPrandtlNumber(i,j,k)  = 1. _d 0
141           TKEPrandtlNumber(I,J,K)  = 0. _d 0           GGL90mixingLength(i,j,k) = GGL90mixingLengthMin
142           GGL90mixingLength(I,J,K) = GGL90mixingLengthMin           GGL90visctmp(i,j,k)      = 0. _d 0
143               rMixingLength(I,J,K) = 0. _d 0  #ifndef SOLVE_DIAGONAL_LOWMEMORY
144             a3d(i,j,k) = 0. _d 0
145             b3d(i,j,k) = 1. _d 0
146             c3d(i,j,k) = 0. _d 0
147    #endif
148          ENDDO          ENDDO
149         ENDDO         ENDDO
150        ENDDO        ENDDO
151        DO J=1-Oly,sNy+Oly        DO j=1-OLy,sNy+OLy
152         DO I=1-Olx,sNx+Olx         DO i=1-OLx,sNx+OLx
153          rhoK    (I,J) = 0. _d 0          totalDepth(i,j)    = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
154          rhoKm1  (I,J) = 0. _d 0          rMixingLength(i,j,1) = 0. _d 0
155          totalDepth(I,J)   = 0. _d 0          mxLength_Dn(i,j,1) = GGL90mixingLengthMin
156          IF ( recip_Rcol(I,J,bi,bj) .NE. 0. )          SQRTTKE(i,j,1) = SQRT( GGL90TKE(i,j,1,bi,bj) )
      &       totalDepth(I,J) = 1./recip_Rcol(I,J,bi,bj)  
157         ENDDO         ENDDO
158        ENDDO        ENDDO
159    
160  C     start k-loop  C     start k-loop
161        DO K = 2, Nr        DO k = 2, Nr
162         Km1 = K-1  c      km1 = k-1
163         Kp1 = MIN(Nr,K+1)  c      kp1 = MIN(Nr,k+1)
164         CALL FIND_RHO_2D(         DO j=jMin,jMax
165       I      iMin, iMax, jMin, jMax, K,          DO i=iMin,iMax
166       I      theta(1-OLx,1-OLy,Km1,bi,bj), salt(1-OLx,1-OLy,Km1,bi,bj),           SQRTTKE(i,j,k)=SQRT( GGL90TKE(i,j,k,bi,bj) )
167       O      rhoKm1,  
      I      Km1, bi, bj, myThid )  
   
        CALL FIND_RHO_2D(  
      I      iMin, iMax, jMin, jMax, K,  
      I      theta(1-OLx,1-OLy,K,bi,bj), salt(1-OLx,1-OLy,K,bi,bj),  
      O      rhoK,  
      I      K, bi, bj, myThid )  
        DO J=jMin,jMax  
         DO I=iMin,iMax  
          SQRTTKE=SQRT( GGL90TKE(I,J,K,bi,bj) )  
 C  
168  C     buoyancy frequency  C     buoyancy frequency
169  C           Nsquare(i,j,k) = gravity*gravitySign*recip_rhoConst
170           Nsquare = - gravity*recip_rhoConst*recip_drC(K)       &                  * sigmaR(i,j,k)
171       &        * ( rhoKm1(I,J) - rhoK(I,J) )*maskC(I,J,K,bi,bj)  cC     vertical shear term (dU/dz)^2+(dV/dz)^2
172  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)
173           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)) )
174       &             - (uVel(I,J,K  ,bi,bj)+uVel(I+1,J,K  ,bi,bj)) )  c     &        *recip_drC(k)
175       &        *recip_drC(K)  c         tempV= .5 _d 0*( vVel(i,j,km1,bi,bj)+vVel(i,j+1,km1,bi,bj)
176           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)) )
177       &             - (vVel(I,J,K  ,bi,bj)+vVel(I,J+1,K  ,bi,bj)) )  c     &        *recip_drC(k)
178       &        *recip_drC(K)  c         verticalShear = tempU*tempU + tempV*tempV
179           verticalShear = tempU*tempU + tempV*tempV  c         RiNumber   = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)
180           RiNumber   = MAX(Nsquare,0. _d 0)/(verticalShear+GGL90eps)  cC     compute Prandtl number (always greater than 0)
181  C     compute Prandtl number (always greater than 0)  c         prTemp = 1. _d 0
182           prTemp = 1. _d 0  c         IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
183           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.0 _d 0,prTemp)  
184  C     mixing length  C     mixing length
185           GGL90mixingLength(I,J,K) =           GGL90mixingLength(i,j,k) = SQRTTWO *
186       &        SQRTTKE/SQRT( MAX(Nsquare,GGL90eps) )       &        SQRTTKE(i,j,k)/SQRT( MAX(Nsquare(i,j,k),GGL90eps) )
 C     impose upper bound for mixing length (total depth)  
          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)  
          rMixingLength(I,J,K) = 1. _d 0 /GGL90mixingLength(I,J,K)  
 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*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  
      &        - KappaM*Nsquare/TKEPrandtlNumber(I,J,K)  
      &        - TKEdissipation  
      &        )  
187          ENDDO          ENDDO
188         ENDDO         ENDDO
189        ENDDO        ENDDO
190  C     horizontal diffusion of TKE (requires an exchange in  
191  C      do_fields_blocking_exchanges)  C- ensure mixing between first and second level
192          IF (mxlSurfFlag) THEN
193           DO j=jMin,jMax
194            DO i=iMin,iMax
195             GGL90mixingLength(i,j,2)=drF(1)
196            ENDDO
197           ENDDO
198          ENDIF
199    
200    C- Impose upper and lower bound for mixing length
201          IF ( mxlMaxFlag .EQ. 0 ) THEN
202    
203           DO k=2,Nr
204            DO j=jMin,jMax
205             DO i=iMin,iMax
206              MaxLength=totalDepth(i,j)
207              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
208         &                                   MaxLength)
209             ENDDO
210            ENDDO
211           ENDDO
212    
213           DO k=2,Nr
214            DO j=jMin,jMax
215             DO i=iMin,iMax
216              GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
217         &                                   GGL90mixingLengthMin)
218              rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
219             ENDDO
220            ENDDO
221           ENDDO
222    
223          ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN
224    
225           DO k=2,Nr
226            DO j=jMin,jMax
227             DO i=iMin,iMax
228              MaxLength=MIN(Ro_surf(i,j,bi,bj)-rF(k),rF(k)-R_low(i,j,bi,bj))
229    c         MaxLength=MAX(MaxLength,20. _d 0)
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. 2 ) THEN
247    
248           DO k=2,Nr
249            DO j=jMin,jMax
250             DO i=iMin,iMax
251              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
252         &        GGL90mixingLength(i,j,k-1)+drF(k-1))
253             ENDDO
254            ENDDO
255           ENDDO
256           DO j=jMin,jMax
257            DO i=iMin,iMax
258              GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
259         &       GGL90mixingLengthMin+drF(Nr))
260            ENDDO
261           ENDDO
262           DO k=Nr-1,2,-1
263            DO j=jMin,jMax
264             DO i=iMin,iMax
265              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
266         &        GGL90mixingLength(i,j,k+1)+drF(k))
267             ENDDO
268            ENDDO
269           ENDDO
270    
271           DO k=2,Nr
272            DO j=jMin,jMax
273             DO i=iMin,iMax
274              GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
275         &                                   GGL90mixingLengthMin)
276              rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
277             ENDDO
278            ENDDO
279           ENDDO
280    
281          ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN
282    
283           DO k=2,Nr
284            DO j=jMin,jMax
285             DO i=iMin,iMax
286              mxLength_Dn(i,j,k) = MIN(GGL90mixingLength(i,j,k),
287         &        mxLength_Dn(i,j,k-1)+drF(k-1))
288             ENDDO
289            ENDDO
290           ENDDO
291           DO j=jMin,jMax
292            DO i=iMin,iMax
293              GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
294         &       GGL90mixingLengthMin+drF(Nr))
295            ENDDO
296           ENDDO
297           DO k=Nr-1,2,-1
298            DO j=jMin,jMax
299             DO i=iMin,iMax
300              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
301         &        GGL90mixingLength(i,j,k+1)+drF(k))
302             ENDDO
303            ENDDO
304           ENDDO
305    
306           DO k=2,Nr
307            DO j=jMin,jMax
308             DO i=iMin,iMax
309              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
310         &                                  mxLength_Dn(i,j,k))
311              tmpmlx = SQRT( GGL90mixingLength(i,j,k)*mxLength_Dn(i,j,k) )
312              tmpmlx = MAX( tmpmlx, GGL90mixingLengthMin)
313              rMixingLength(i,j,k) = 1. _d 0 / tmpmlx
314             ENDDO
315            ENDDO
316           ENDDO
317    
318          ELSE
319           STOP 'GGL90_CALC: Wrong mxlMaxFlag (mixing length limit)'
320          ENDIF
321    
322    C- Impose minimum mixing length (to avoid division by zero)
323    c      DO k=2,Nr
324    c      DO j=jMin,jMax
325    c       DO i=iMin,iMax
326    c        GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
327    c    &        GGL90mixingLengthMin)
328    c        rMixingLength(i,j,k) = 1. _d 0 /GGL90mixingLength(i,j,k)
329    c       ENDDO
330    c      ENDDO
331    c     ENDDO
332    
333          DO k=2,Nr
334           km1 = k-1
335    
336  #ifdef ALLOW_GGL90_HORIZDIFF  #ifdef ALLOW_GGL90_HORIZDIFF
337        IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN        IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
338         DO K=2,Nr  C     horizontal diffusion of TKE (requires an exchange in
339    C      do_fields_blocking_exchanges)
340  C     common factors  C     common factors
341          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
342           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
343            xA(i,j) = _dyG(i,j,bi,bj)            xA(i,j) = _dyG(i,j,bi,bj)*drC(k)*
344       &         *drF(k)*_hFacW(i,j,k,bi,bj)       &                 (min(.5 _d 0,_hFacW(i,j,k-1,bi,bj) ) +
345            yA(i,j) = _dxG(i,j,bi,bj)       &                  min(.5 _d 0,_hFacW(i,j,k  ,bi,bj) ) )
346       &         *drF(k)*_hFacS(i,j,k,bi,bj)            yA(i,j) = _dxG(i,j,bi,bj)*drC(k)*
347         &                 (min(.5 _d 0,_hFacS(i,j,k-1,bi,bj) ) +
348         &                  min(.5 _d 0,_hFacS(i,j,k  ,bi,bj) ) )
349           ENDDO           ENDDO
350          ENDDO          ENDDO
351  C     Compute diffusive fluxes  C     Compute diffusive fluxes
352  C     ... across x-faces  C     ... across x-faces
353          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
354           dfx(1-Olx,j)=0.           dfx(1-OLx,j)=0. _d 0
355           DO i=1-Olx+1,sNx+Olx           DO i=1-OLx+1,sNx+OLx
356            dfx(i,j) = -GGL90diffTKEh*xA(i,j)            dfx(i,j) = -GGL90diffTKEh*xA(i,j)
357       &      *_recip_dxC(i,j,bi,bj)       &      *_recip_dxC(i,j,bi,bj)
358       &      *(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))
359    #ifdef ISOTROPIC_COS_SCALING
360       &      *CosFacU(j,bi,bj)       &      *CosFacU(j,bi,bj)
361    #endif /* ISOTROPIC_COS_SCALING */
362           ENDDO           ENDDO
363          ENDDO          ENDDO
364  C     ... across y-faces  C     ... across y-faces
365          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
366           dfy(i,1-Oly)=0.           dfy(i,1-OLy)=0. _d 0
367          ENDDO          ENDDO
368          DO j=1-Oly+1,sNy+Oly          DO j=1-OLy+1,sNy+OLy
369           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
370            dfy(i,j) = -GGL90diffTKEh*yA(i,j)            dfy(i,j) = -GGL90diffTKEh*yA(i,j)
371       &      *_recip_dyC(i,j,bi,bj)       &      *_recip_dyC(i,j,bi,bj)
372       &      *(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 246  C     ... across y-faces Line 376  C     ... across y-faces
376           ENDDO           ENDDO
377          ENDDO          ENDDO
378  C     Compute divergence of fluxes  C     Compute divergence of fluxes
379          DO j=1-Oly,sNy+Oly-1          DO j=1-OLy,sNy+OLy-1
380           DO i=1-Olx,sNx+Olx-1           DO i=1-OLx,sNx+OLx-1
381            gTKE(i,j,k)=gTKE(i,j,k)            hFac = min(.5 _d 0,_hFacC(i,j,k-1,bi,bj) ) +
382       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj)       &          min(.5 _d 0,_hFacC(i,j,k  ,bi,bj) )
383       &         *( (dfx(i+1,j)-dfx(i,j))            gTKE(i,j) = 0.0
384       &           +(dfy(i,j+1)-dfy(i,j))            if ( hFac .ne. 0.0 )
385       &           )       &      gTKE(i,j) = -recip_drC(k)*recip_rA(i,j,bi,bj)/hFac
386         &         *((dfx(i+1,j)-dfx(i,j))
387         &          +(dfy(i,j+1)-dfy(i,j)) )
388           ENDDO           ENDDO
389          ENDDO          ENDDO
390  C       end of k-loop  C      end if GGL90diffTKEh .eq. 0.
391           ENDIF
392    #endif /* ALLOW_GGL90_HORIZDIFF */
393    
394           DO j=jMin,jMax
395            DO i=iMin,iMax
396    C     vertical shear term (dU/dz)^2+(dV/dz)^2
397             tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)
398         &                 -( uVel(i,j,k  ,bi,bj)+uVel(i+1,j,k  ,bi,bj)) )
399         &        *recip_drC(k)
400             tempV= .5 _d 0*( vVel(i,j,km1,bi,bj)+vVel(i,j+1,km1,bi,bj)
401         &                 -( vVel(i,j,k  ,bi,bj)+vVel(i,j+1,k  ,bi,bj)) )
402         &        *recip_drC(k)
403             verticalShear = tempU*tempU + tempV*tempV
404             RiNumber = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)
405    C     compute Prandtl number (always greater than 0)
406             prTemp = 1. _d 0
407             IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
408             TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
409    c         TKEPrandtlNumber(i,j,k) = 1. _d 0
410    
411    C     viscosity and diffusivity
412             KappaM = GGL90ck*GGL90mixingLength(i,j,k)*SQRTTKE(i,j,k)
413             GGL90visctmp(i,j,k) = MAX(KappaM,diffKrNrT(k))
414         &                            * maskC(i,j,k,bi,bj)
415    c        note: storing GGL90visctmp like this, and using it later to compute
416    c              GGL9rdiffKr etc. is robust in case of smoothing (e.g. see OPA)
417             KappaM = MAX(KappaM,viscArNr(k)) * maskC(i,j,k,bi,bj)
418             KappaH = KappaM/TKEPrandtlNumber(i,j,k)
419             KappaE(i,j,k) = GGL90alpha * KappaM * maskC(i,j,k,bi,bj)
420    
421    C     dissipation term
422             TKEdissipation = explDissFac*GGL90ceps
423         &        *SQRTTKE(i,j,k)*rMixingLength(i,j,k)
424         &        *GGL90TKE(i,j,k,bi,bj)
425    C     partial update with sum of explicit contributions
426             GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
427         &        + deltaTggl90*(
428         &        + KappaM*verticalShear
429         &        - KappaH*Nsquare(i,j,k)
430         &        - TKEdissipation
431         &        )
432            ENDDO
433         ENDDO         ENDDO
434  C     end if GGL90diffTKEh .eq. 0.  
435        ENDIF  #ifdef ALLOW_GGL90_HORIZDIFF
436           IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
437    C--    Add horiz. diffusion tendency
438            DO j=jMin,jMax
439             DO i=iMin,iMax
440              GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
441         &                          + gTKE(i,j)*deltaTggl90
442             ENDDO
443            ENDDO
444           ENDIF
445  #endif /* ALLOW_GGL90_HORIZDIFF */  #endif /* ALLOW_GGL90_HORIZDIFF */
446    
447    C--   end of k loop
448          ENDDO
449    
450  C     ============================================  C     ============================================
451  C     Implicit time step to update TKE for k=1,Nr;  C     Implicit time step to update TKE for k=1,Nr;
452  C     TKE(Nr+1)=0 by default  C     TKE(Nr+1)=0 by default
# Line 269  C     set up matrix Line 455  C     set up matrix
455  C--   Lower diagonal  C--   Lower diagonal
456        DO j=jMin,jMax        DO j=jMin,jMax
457         DO i=iMin,iMax         DO i=iMin,iMax
458           a(i,j,1) = 0. _d 0           a3d(i,j,1) = 0. _d 0
459         ENDDO         ENDDO
460        ENDDO        ENDDO
461        DO k=2,Nr        DO k=2,Nr
462         km1=MAX(1,k-1)         km1=MAX(2,k-1)
463         DO j=jMin,jMax         DO j=jMin,jMax
464          DO i=iMin,iMax          DO i=iMin,iMax
465           a(i,j,k) = -deltaTggl90  C-    We keep recip_hFacC in the diffusive flux calculation,
466       &        *recip_drF(km1)*recip_hFacI(i,j,k,bi,bj)  C-    but no hFacC in TKE volume control
467       &        *.5*(KappaE(i,j, k )+KappaE(i,j,km1))  C-    No need for maskC(k-1) with recip_hFacC(k-1)
468       &        *recip_drC(k)           a3d(i,j,k) = -deltaTggl90
469            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)
470         &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))
471         &        *recip_drC(k)*maskC(i,j,k,bi,bj)
472          ENDDO          ENDDO
473         ENDDO         ENDDO
474        ENDDO        ENDDO
475  C--   Upper diagonal  C--   Upper diagonal
476        DO j=jMin,jMax        DO j=jMin,jMax
477         DO i=iMin,iMax         DO i=iMin,iMax
478           c(i,j,1)  = 0. _d 0           c3d(i,j,1)  = 0. _d 0
          c(i,j,Nr) = 0. _d 0  
479         ENDDO         ENDDO
480        ENDDO        ENDDO
481  CML      DO k=1,Nr-1        DO k=2,Nr
       DO k=2,Nr-1  
        kp1=min(Nr,k+1)  
482         DO j=jMin,jMax         DO j=jMin,jMax
483          DO i=iMin,iMax          DO i=iMin,iMax
484            c(i,j,k) = -deltaTggl90            kp1=MAX(1,MIN(klowC(i,j,bi,bj),k+1))
485       &        *recip_drF( k )*recip_hFacI(i,j,k,bi,bj)  C-    We keep recip_hFacC in the diffusive flux calculation,
486       &               *.5*(KappaE(i,j,k)+KappaE(i,j,kp1))  C-    but no hFacC in TKE volume control
487       &        *recip_drC(k)  C-    No need for maskC(k) with recip_hFacC(k)
488            IF (recip_hFacI(i,j,kp1,bi,bj).EQ.0.) c(i,j,k)=0.            c3d(i,j,k) = -deltaTggl90
489         &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)
490         &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))
491         &        *recip_drC(k)*maskC(i,j,k-1,bi,bj)
492          ENDDO          ENDDO
493         ENDDO         ENDDO
494        ENDDO        ENDDO
495  C--   Center diagonal  C--   Center diagonal
496        DO k=1,Nr        DO k=1,Nr
497           km1 = MAX(k-1,1)
498         DO j=jMin,jMax         DO j=jMin,jMax
499          DO i=iMin,iMax          DO i=iMin,iMax
500            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)
501       &        + ab15*deltaTggl90*GGL90ceps*SQRT(GGL90TKE(I,J,K,bi,bj))       &        + implDissFac*deltaTggl90*GGL90ceps*SQRTTKE(i,j,k)
502       &        *rMixingLength(I,J,K)       &        * rMixingLength(i,j,k)
503         &        * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)
504           ENDDO           ENDDO
505         ENDDO         ENDDO
506        ENDDO        ENDDO
507  C     end set up matrix  C     end set up matrix
508    
 C  
509  C     Apply boundary condition  C     Apply boundary condition
510  C        kp1 = MIN(Nr,kSurf+1)
511        DO J=jMin,jMax        DO j=jMin,jMax
512         DO I=iMin,iMax         DO i=iMin,iMax
513  C     estimate friction velocity uStar from surface forcing  C     estimate friction velocity uStar from surface forcing
514          uStarSquare = SQRT(          uStarSquare = SQRT(
515       &         ( .5*( surfaceForcingU(I,  J,  bi,bj)       &    ( .5 _d 0*( surfaceForcingU(i,  j,  bi,bj)
516       &              + surfaceForcingU(I+1,J,  bi,bj) ) )**2       &              + surfaceForcingU(i+1,j,  bi,bj) ) )**2
517       &       + ( .5*( surfaceForcingV(I,  J,  bi,bj)       &  + ( .5 _d 0*( surfaceForcingV(i,  j,  bi,bj)
518       &              + surfaceForcingV(I,  J+1,bi,bj) ) )**2       &              + surfaceForcingV(i,  j+1,bi,bj) ) )**2
519       &                     )       &                     )
520  C     Dirichlet surface boundary condition for TKE  C     Dirichlet surface boundary condition for TKE
521          gTKE(I,J,kSurf) = MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)          GGL90TKE(i,j,kSurf,bi,bj) = maskC(i,j,kSurf,bi,bj)
522       &                     *maskC(I,J,kSurf,bi,bj)       &           *MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)
523            GGL90TKE(i,j,kp1,bi,bj) = GGL90TKE(i,j,kp1,bi,bj)
524         &               - a3d(i,j,kp1)*GGL90TKE(i,j,kSurf,bi,bj)
525            a3d(i,j,kp1) = 0. _d 0
526  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom
527          kBottom   = MIN(MAX(kLowC(I,J,bi,bj),1),Nr)          kBottom   = MAX(kLowC(i,j,bi,bj),1)
528          gTKE(I,J,kBottom) = gTKE(I,J,kBottom)          GGL90TKE(i,j,kBottom,bi,bj) = GGL90TKE(i,j,kBottom,bi,bj)
529       &       - GGL90TKEbottom*c(I,J,kBottom)       &                              - GGL90TKEbottom*c3d(i,j,kBottom)
530         ENDDO          c3d(i,j,kBottom) = 0. _d 0
531        ENDDO         ENDDO
532  C        ENDDO
533  C     solve tri-diagonal system, and store solution on gTKE (previously rhs)  
534  C  C     solve tri-diagonal system
535        CALL GGL90_SOLVE( bi, bj, iMin, iMax, jMin, jMax,        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
536       I     a, b, c,       I                        a3d, b3d, c3d,
537       U     gTKE,       U                        GGL90TKE(1-OLx,1-OLy,1,bi,bj),
538       I     myThid )       O                        errCode,
539  C       I                        bi, bj, myThid )
540  C     now update TKE  
541  C        DO k=1,Nr
542        DO K=1,Nr         DO j=jMin,jMax
543         DO J=jMin,jMax          DO i=iMin,iMax
         DO I=iMin,iMax  
544  C     impose minimum TKE to avoid numerical undershoots below zero  C     impose minimum TKE to avoid numerical undershoots below zero
545           GGL90TKE(I,J,K,bi,bj) = MAX( gTKE(I,J,K), GGL90TKEmin )           GGL90TKE(i,j,k,bi,bj) = maskC(i,j,k,bi,bj)
546       &        * maskC(I,J,K,bi,bj)       &                  *MAX( GGL90TKE(i,j,k,bi,bj), GGL90TKEmin )
547          ENDDO          ENDDO
548         ENDDO         ENDDO
549        ENDDO        ENDDO
550  C  
551  C     end of time step  C     end of time step
552  C     ===============================  C     ===============================
553  C     compute viscosity coefficients  
554  C        DO k=2,Nr
555        DO K=2,Nr         DO j=1,sNy
556         DO J=jMin,jMax          DO i=1,sNx
557          DO I=iMin,iMax  #ifdef ALLOW_GGL90_SMOOTH
558  C     Eq. (11), (18) and (21)           tmpVisc=
559           KappaM = GGL90ck*GGL90mixingLength(I,J,K)*       &  (
560       &                  SQRT( GGL90TKE(I,J,K,bi,bj) )       &   p4 *  GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
561           KappaH = KappaM/TKEPrandtlNumber(I,J,K)       &  +p8 *( GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj)
562  C     Set a minium (= background) value       &       + GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj)
563           KappaM = MAX(KappaM,viscAr)       &       + GGL90visctmp(i+1,j  ,k) * mskCor(i+1,j  ,bi,bj)
564           KappaH = MAX(KappaH,diffKrNrT(k))       &       + GGL90visctmp(i  ,j+1,k) * mskCor(i  ,j+1,bi,bj))
565  C     Set a maximum and mask land point       &  +p16*( GGL90visctmp(i+1,j+1,k) * mskCor(i+1,j+1,bi,bj)
566           GGL90viscAr(I,J,K,bi,bj) = MIN(KappaM,GGL90viscMax)       &       + GGL90visctmp(i+1,j-1,k) * mskCor(i+1,j-1,bi,bj)
567       &        * maskC(I,J,K,bi,bj)       &       + GGL90visctmp(i-1,j+1,k) * mskCor(i-1,j+1,bi,bj)
568           GGL90diffKr(I,J,K,bi,bj) = MIN(KappaH,GGL90diffMax)       &       + GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj))
569       &        * maskC(I,J,K,bi,bj)       &  )
570         & /(p4
571         &  +p8 *(       maskC(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)
572         &       +       maskC(i  ,j-1,k,bi,bj) * mskCor(i  ,j-1,bi,bj)
573         &       +       maskC(i+1,j  ,k,bi,bj) * mskCor(i+1,j  ,bi,bj)
574         &       +       maskC(i  ,j+1,k,bi,bj) * mskCor(i  ,j+1,bi,bj))
575         &  +p16*(       maskC(i+1,j+1,k,bi,bj) * mskCor(i+1,j+1,bi,bj)
576         &       +       maskC(i+1,j-1,k,bi,bj) * mskCor(i+1,j-1,bi,bj)
577         &       +       maskC(i-1,j+1,k,bi,bj) * mskCor(i-1,j+1,bi,bj)
578         &       +       maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj))
579         &  )*maskC(i,j,k,bi,bj)*mskCor(i,j,bi,bj)
580    #else
581             tmpVisc = GGL90visctmp(i,j,k)
582    #endif
583             tmpVisc = MIN(tmpVisc/TKEPrandtlNumber(i,j,k),GGL90diffMax)
584             GGL90diffKr(i,j,k,bi,bj)= MAX( tmpVisc , diffKrNrT(k) )
585            ENDDO
586           ENDDO
587          ENDDO
588    
589          DO k=2,Nr
590           DO j=1,sNy
591            DO i=1,sNx+1
592    #ifdef ALLOW_GGL90_SMOOTH
593            tmpVisc =
594         & (
595         &   p4 *(GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
596         &       +GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj))
597         &  +p8 *(GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj)
598         &       +GGL90visctmp(i-1,j+1,k) * mskCor(i-1,j+1,bi,bj)
599         &       +GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj)
600         &       +GGL90visctmp(i  ,j+1,k) * mskCor(i  ,j+1,bi,bj))
601         &  )
602         & /(p4 * 2. _d 0
603         &  +p8 *(      maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj)
604         &       +      maskC(i-1,j+1,k,bi,bj) * mskCor(i-1,j+1,bi,bj)
605         &       +      maskC(i  ,j-1,k,bi,bj) * mskCor(i  ,j-1,bi,bj)
606         &       +      maskC(i  ,j+1,k,bi,bj) * mskCor(i  ,j+1,bi,bj))
607         &  )
608         &  *maskC(i  ,j,k,bi,bj)*mskCor(i  ,j,bi,bj)
609         &  *maskC(i-1,j,k,bi,bj)*mskCor(i-1,j,bi,bj)
610    #else
611            tmpVisc = _maskW(i,j,k,bi,bj) *
612         &                   (.5 _d 0*(GGL90visctmp(i,j,k)
613         &                            +GGL90visctmp(i-1,j,k))
614         &                   )
615    #endif
616            tmpVisc = MIN( tmpVisc , GGL90viscMax )
617            GGL90viscArU(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
618          ENDDO          ENDDO
619         ENDDO         ENDDO
 C     end third k-loop  
620        ENDDO        ENDDO
621    
622          DO k=2,Nr
623           DO j=1,sNy+1
624            DO i=1,sNx
625    #ifdef ALLOW_GGL90_SMOOTH
626            tmpVisc =
627         & (
628         &   p4 *(GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
629         &       +GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj))
630         &  +p8 *(GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj)
631         &       +GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj)
632         &       +GGL90visctmp(i+1,j  ,k) * mskCor(i+1,j  ,bi,bj)
633         &       +GGL90visctmp(i+1,j-1,k) * mskCor(i+1,j-1,bi,bj))
634         &  )
635         & /(p4 * 2. _d 0
636         &  +p8 *(      maskC(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)
637         &       +      maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj)
638         &       +      maskC(i+1,j  ,k,bi,bj) * mskCor(i+1,j  ,bi,bj)
639         &       +      maskC(i+1,j-1,k,bi,bj) * mskCor(i+1,j-1,bi,bj))
640         &  )
641         &   *maskC(i,j  ,k,bi,bj)*mskCor(i,j  ,bi,bj)
642         &   *maskC(i,j-1,k,bi,bj)*mskCor(i,j-1,bi,bj)
643    #else
644            tmpVisc = _maskS(i,j,k,bi,bj) *
645         &                   (.5 _d 0*(GGL90visctmp(i,j,k)
646         &                            +GGL90visctmp(i,j-1,k))
647         &                   )
648    
649    #endif
650            tmpVisc = MIN( tmpVisc , GGL90viscMax )
651            GGL90viscArV(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
652            ENDDO
653           ENDDO
654          ENDDO
655    
656    #ifdef ALLOW_DIAGNOSTICS
657          IF ( useDiagnostics ) THEN
658             CALL DIAGNOSTICS_FILL( GGL90TKE   ,'GGL90TKE',
659         &                          0,Nr, 1, bi, bj, myThid )
660             CALL DIAGNOSTICS_FILL( GGL90viscArU,'GGL90ArU',
661         &                          0,Nr, 1, bi, bj, myThid )
662             CALL DIAGNOSTICS_FILL( GGL90viscArV,'GGL90ArV',
663         &                          0,Nr, 1, bi, bj, myThid )
664             CALL DIAGNOSTICS_FILL( GGL90diffKr,'GGL90Kr ',
665         &                          0,Nr, 1, bi, bj, myThid )
666             CALL DIAGNOSTICS_FILL( TKEPrandtlNumber ,'GGL90Prl',
667         &                          0,Nr, 2, bi, bj, myThid )
668             CALL DIAGNOSTICS_FILL( GGL90mixingLength,'GGL90Lmx',
669         &                          0,Nr, 2, bi, bj, myThid )
670          ENDIF
671    #endif
672    
673  #endif /* ALLOW_GGL90 */  #endif /* ALLOW_GGL90 */
674    
675        RETURN        RETURN
676        END        END
   
Legend:
Removed from v.1.8  
changed lines
  Added in v.1.26
  ViewVC Help
Powered by ViewVC 1.1.22