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