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