/[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.9 by dfer, Mon Sep 29 13:47:23 2008 UTC revision 1.19 by jmc, Thu Aug 19 23:52:37 2010 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     a(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
102        _RL     b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
103        _RL     c(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     c(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
              rMixingLength(I,J,K) = 0. _d 0  
142          ENDDO          ENDDO
143         ENDDO         ENDDO
144        ENDDO        ENDDO
145        DO J=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
146         DO I=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
147          rhoK    (I,J) = 0. _d 0          rhoK(i,j)          = 0. _d 0
148          rhoKm1  (I,J) = 0. _d 0          rhoKm1(i,j)        = 0. _d 0
149          totalDepth(I,J)   = 0. _d 0          totalDepth(i,j)    = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
150          IF ( recip_Rcol(I,J,bi,bj) .NE. 0. )          rMixingLength(i,j,1) = 0. _d 0
151       &       totalDepth(I,J) = 1./recip_Rcol(I,J,bi,bj)          mxLength_Dn(i,j,1) = GGL90mixingLengthMin
152            SQRTTKE(i,j,1) = SQRT( GGL90TKE(i,j,1,bi,bj) )
153         ENDDO         ENDDO
154        ENDDO        ENDDO
155    
156  C     start k-loop  C     start k-loop
157        DO K = 2, Nr        DO k = 2, Nr
158         Km1 = K-1         km1 = k-1
159         Kp1 = MIN(Nr,K+1)  c      kp1 = MIN(Nr,k+1)
160         CALL FIND_RHO_2D(         CALL FIND_RHO_2D(
161       I      iMin, iMax, jMin, jMax, K,       I      iMin, iMax, jMin, jMax, k,
162       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),
163       O      rhoKm1,       O      rhoKm1,
164       I      Km1, bi, bj, myThid )       I      km1, bi, bj, myThid )
165    
166         CALL FIND_RHO_2D(         CALL FIND_RHO_2D(
167       I      iMin, iMax, jMin, jMax, K,       I      iMin, iMax, jMin, jMax, k,
168       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),
169       O      rhoK,       O      rhoK,
170       I      K, bi, bj, myThid )       I      k, bi, bj, myThid )
171         DO J=jMin,jMax         DO j=jMin,jMax
172          DO I=iMin,iMax          DO i=iMin,iMax
173           SQRTTKE=SQRT( GGL90TKE(I,J,K,bi,bj) )           SQRTTKE(i,j,k)=SQRT( GGL90TKE(i,j,k,bi,bj) )
174  C  
175  C     buoyancy frequency  C     buoyancy frequency
176  C           Nsquare(i,j,k) = - gravity*recip_rhoConst*recip_drC(k)
177           Nsquare = - gravity*recip_rhoConst*recip_drC(K)       &        * ( rhoKm1(i,j) - rhoK(i,j) )*maskC(i,j,k,bi,bj)
178       &        * ( rhoKm1(I,J) - rhoK(I,J) )*maskC(I,J,K,bi,bj)  cC     vertical shear term (dU/dz)^2+(dV/dz)^2
179  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)
180           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)) )
181       &             - (uVel(I,J,K  ,bi,bj)+uVel(I+1,J,K  ,bi,bj)) )  c     &        *recip_drC(k)
182       &        *recip_drC(K)  c         tempV= .5 _d 0*( vVel(i,j,km1,bi,bj)+vVel(i,j+1,km1,bi,bj)
183           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)) )
184       &             - (vVel(I,J,K  ,bi,bj)+vVel(I,J+1,K  ,bi,bj)) )  c     &        *recip_drC(k)
185       &        *recip_drC(K)  c         verticalShear = tempU*tempU + tempV*tempV
186           verticalShear = tempU*tempU + tempV*tempV  c         RiNumber   = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)
187           RiNumber   = MAX(Nsquare,0. _d 0)/(verticalShear+GGL90eps)  cC     compute Prandtl number (always greater than 0)
188  C     compute Prandtl number (always greater than 0)  c         prTemp = 1. _d 0
189           prTemp = 1. _d 0  c         IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
190           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)  
191  C     mixing length  C     mixing length
192           GGL90mixingLength(I,J,K) = SQRTTWO *           GGL90mixingLength(i,j,k) = SQRTTWO *
193       &        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  
      &        )  
194          ENDDO          ENDDO
195         ENDDO         ENDDO
196        ENDDO        ENDDO
197    
198    C- Impose upper and lower bound for mixing length
199          IF ( mxlMaxFlag .EQ. 0 ) THEN
200    C-
201           DO k=2,Nr
202            DO j=jMin,jMax
203             DO i=iMin,iMax
204              MaxLength=totalDepth(i,j)
205              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
206         &                                   MaxLength)
207             ENDDO
208            ENDDO
209           ENDDO
210    
211           DO k=2,Nr
212            DO j=jMin,jMax
213             DO i=iMin,iMax
214              GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
215         &                                   GGL90mixingLengthMin)
216              rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
217             ENDDO
218            ENDDO
219           ENDDO
220    
221          ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN
222    C-
223           DO k=2,Nr
224            DO j=jMin,jMax
225             DO i=iMin,iMax
226              MaxLength=MIN(Ro_surf(i,j,bi,bj)-rF(k),rF(k)-R_low(i,j,bi,bj))
227    c         MaxLength=MAX(MaxLength,20. _d 0)
228              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
229         &                                   MaxLength)
230             ENDDO
231            ENDDO
232           ENDDO
233    
234           DO k=2,Nr
235            DO j=jMin,jMax
236             DO i=iMin,iMax
237              GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
238         &                                   GGL90mixingLengthMin)
239              rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
240             ENDDO
241            ENDDO
242           ENDDO
243    
244          ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN
245    C-
246    cgf ensure mixing between first and second level
247    c      DO j=jMin,jMax
248    c        DO i=iMin,iMax
249    c         GGL90mixingLength(i,j,2)=drF(1)
250    c        ENDDO
251    c      ENDDO
252    cgf
253           DO k=2,Nr
254            DO j=jMin,jMax
255             DO i=iMin,iMax
256              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
257         &        GGL90mixingLength(i,j,k-1)+drF(k-1))
258             ENDDO
259            ENDDO
260           ENDDO
261           DO j=jMin,jMax
262            DO i=iMin,iMax
263              GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
264         &       GGL90mixingLengthMin+drF(Nr))
265            ENDDO
266           ENDDO
267           DO k=Nr-1,2,-1
268            DO j=jMin,jMax
269             DO i=iMin,iMax
270              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
271         &        GGL90mixingLength(i,j,k+1)+drF(k))
272             ENDDO
273            ENDDO
274           ENDDO
275    
276           DO k=2,Nr
277            DO j=jMin,jMax
278             DO i=iMin,iMax
279              GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
280         &                                   GGL90mixingLengthMin)
281              rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
282             ENDDO
283            ENDDO
284           ENDDO
285    
286          ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN
287    C-
288           DO k=2,Nr
289            DO j=jMin,jMax
290             DO i=iMin,iMax
291              mxLength_Dn(i,j,k) = MIN(GGL90mixingLength(i,j,k),
292         &        mxLength_Dn(i,j,k-1)+drF(k-1))
293             ENDDO
294            ENDDO
295           ENDDO
296           DO j=jMin,jMax
297            DO i=iMin,iMax
298              GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
299         &       GGL90mixingLengthMin+drF(Nr))
300            ENDDO
301           ENDDO
302           DO k=Nr-1,2,-1
303            DO j=jMin,jMax
304             DO i=iMin,iMax
305              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
306         &        GGL90mixingLength(i,j,k+1)+drF(k))
307             ENDDO
308            ENDDO
309           ENDDO
310    
311           DO k=2,Nr
312            DO j=jMin,jMax
313             DO i=iMin,iMax
314              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
315         &                                  mxLength_Dn(i,j,k))
316              tmpmlx = SQRT( GGL90mixingLength(i,j,k)*mxLength_Dn(i,j,k) )
317              tmpmlx = MAX( tmpmlx, GGL90mixingLengthMin)
318              rMixingLength(i,j,k) = 1. _d 0 / tmpmlx
319             ENDDO
320            ENDDO
321           ENDDO
322    
323          ELSE
324           STOP 'GGL90_CALC: Wrong mxlMaxFlag (mixing length limit)'
325          ENDIF
326    
327    C- Impose minimum mixing length (to avoid division by zero)
328    c      DO k=2,Nr
329    c      DO j=jMin,jMax
330    c       DO i=iMin,iMax
331    c        GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
332    c    &        GGL90mixingLengthMin)
333    c        rMixingLength(i,j,k) = 1. _d 0 /GGL90mixingLength(i,j,k)
334    c       ENDDO
335    c      ENDDO
336    c     ENDDO
337    
338          DO k=2,Nr
339           km1 = k-1
340    
341    #ifdef ALLOW_GGL90_HORIZDIFF
342           IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
343  C     horizontal diffusion of TKE (requires an exchange in  C     horizontal diffusion of TKE (requires an exchange in
344  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  
345  C     common factors  C     common factors
346          DO j=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
347           DO i=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
# Line 223  C     common factors Line 354  C     common factors
354  C     Compute diffusive fluxes  C     Compute diffusive fluxes
355  C     ... across x-faces  C     ... across x-faces
356          DO j=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
357           dfx(1-Olx,j)=0.           dfx(1-Olx,j)=0. _d 0
358           DO i=1-Olx+1,sNx+Olx           DO i=1-Olx+1,sNx+Olx
359            dfx(i,j) = -GGL90diffTKEh*xA(i,j)            dfx(i,j) = -GGL90diffTKEh*xA(i,j)
360       &      *_recip_dxC(i,j,bi,bj)       &      *_recip_dxC(i,j,bi,bj)
# Line 233  C     ... across x-faces Line 364  C     ... across x-faces
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.           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
# Line 248  C     ... across y-faces Line 379  C     ... across y-faces
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)            gTKE(i,j) =
383       &   -_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)
384       &         *( (dfx(i+1,j)-dfx(i,j))       &         *( (dfx(i+1,j)-dfx(i,j))
385       &           +(dfy(i,j+1)-dfy(i,j))       &           +(dfy(i,j+1)-dfy(i,j))
386       &           )       &          )
387           ENDDO           ENDDO
388          ENDDO          ENDDO
389  C       end of k-loop  C      end if GGL90diffTKEh .eq. 0.
390           ENDIF
391    #endif /* ALLOW_GGL90_HORIZDIFF */
392    
393           DO j=jMin,jMax
394            DO i=iMin,iMax
395    C     vertical shear term (dU/dz)^2+(dV/dz)^2
396             tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)
397         &                 -( uVel(i,j,k  ,bi,bj)+uVel(i+1,j,k  ,bi,bj)) )
398         &        *recip_drC(k)
399             tempV= .5 _d 0*( vVel(i,j,km1,bi,bj)+vVel(i,j+1,km1,bi,bj)
400         &                 -( vVel(i,j,k  ,bi,bj)+vVel(i,j+1,k  ,bi,bj)) )
401         &        *recip_drC(k)
402             verticalShear = tempU*tempU + tempV*tempV
403             RiNumber = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)
404    C     compute Prandtl number (always greater than 0)
405             prTemp = 1. _d 0
406             IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
407             TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
408    c         TKEPrandtlNumber(i,j,k) = 1. _d 0
409    
410    C     viscosity and diffusivity
411             KappaM = GGL90ck*GGL90mixingLength(i,j,k)*SQRTTKE(i,j,k)
412             GGL90visctmp(i,j,k) = MAX(KappaM,diffKrNrT(k))
413         &                            * maskC(i,j,k,bi,bj)
414    c        note: storing GGL90visctmp like this, and using it later to compute
415    c              GGL9rdiffKr etc. is robust in case of smoothing (e.g. see OPA)
416             KappaM = MAX(KappaM,viscArNr(k)) * maskC(i,j,k,bi,bj)
417             KappaH = KappaM/TKEPrandtlNumber(i,j,k)
418             KappaE(i,j,k) = GGL90alpha * KappaM * maskC(i,j,k,bi,bj)
419    
420    C     dissipation term
421             TKEdissipation = explDissFac*GGL90ceps
422         &        *SQRTTKE(i,j,k)*rMixingLength(i,j,k)
423         &        *GGL90TKE(i,j,k,bi,bj)
424    C     partial update with sum of explicit contributions
425             GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
426         &        + deltaTggl90*(
427         &        + KappaM*verticalShear
428         &        - KappaH*Nsquare(i,j,k)
429         &        - TKEdissipation
430         &        )
431            ENDDO
432         ENDDO         ENDDO
433  C     end if GGL90diffTKEh .eq. 0.  
434        ENDIF  #ifdef ALLOW_GGL90_HORIZDIFF
435           IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
436    C--    Add horiz. diffusion tendency
437            DO j=jMin,jMax
438             DO i=iMin,iMax
439              GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
440         &                          + gTKE(i,j)*deltaTggl90
441             ENDDO
442            ENDDO
443           ENDIF
444  #endif /* ALLOW_GGL90_HORIZDIFF */  #endif /* ALLOW_GGL90_HORIZDIFF */
445    
446    C--   end of k loop
447          ENDDO
448    
449  C     ============================================  C     ============================================
450  C     Implicit time step to update TKE for k=1,Nr;  C     Implicit time step to update TKE for k=1,Nr;
451  C     TKE(Nr+1)=0 by default  C     TKE(Nr+1)=0 by default
# Line 273  C--   Lower diagonal Line 458  C--   Lower diagonal
458         ENDDO         ENDDO
459        ENDDO        ENDDO
460        DO k=2,Nr        DO k=2,Nr
461         km1=MAX(1,k-1)         km1=MAX(2,k-1)
462         DO j=jMin,jMax         DO j=jMin,jMax
463          DO i=iMin,iMax          DO i=iMin,iMax
464    C-    We keep recip_hFacC in the diffusive flux calculation,
465    C-    but no hFacC in TKE volume control
466    C-    No need for maskC(k-1) with recip_hFacC(k-1)
467           a(i,j,k) = -deltaTggl90           a(i,j,k) = -deltaTggl90
468       &        *recip_drF(km1)*recip_hFacI(i,j,k,bi,bj)       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)
469       &        *.5*(KappaE(i,j, k )+KappaE(i,j,km1))       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))
470       &        *recip_drC(k)       &        *recip_drC(k)*maskC(i,j,k,bi,bj)
           IF (recip_hFacI(i,j,km1,bi,bj).EQ.0.) a(i,j,k)=0.  
471          ENDDO          ENDDO
472         ENDDO         ENDDO
473        ENDDO        ENDDO
# Line 288  C--   Upper diagonal Line 475  C--   Upper diagonal
475        DO j=jMin,jMax        DO j=jMin,jMax
476         DO i=iMin,iMax         DO i=iMin,iMax
477           c(i,j,1)  = 0. _d 0           c(i,j,1)  = 0. _d 0
          c(i,j,Nr) = 0. _d 0  
478         ENDDO         ENDDO
479        ENDDO        ENDDO
480  CML      DO k=1,Nr-1        DO k=2,Nr
       DO k=2,Nr-1  
        kp1=min(Nr,k+1)  
481         DO j=jMin,jMax         DO j=jMin,jMax
482          DO i=iMin,iMax          DO i=iMin,iMax
483              kp1=MAX(1,MIN(klowC(i,j,bi,bj),k+1))
484    C-    We keep recip_hFacC in the diffusive flux calculation,
485    C-    but no hFacC in TKE volume control
486    C-    No need for maskC(k) with recip_hFacC(k)
487            c(i,j,k) = -deltaTggl90            c(i,j,k) = -deltaTggl90
488       &        *recip_drF( k )*recip_hFacI(i,j,k,bi,bj)       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)
489       &               *.5*(KappaE(i,j,k)+KappaE(i,j,kp1))       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))
490       &        *recip_drC(k)       &        *recip_drC(k)*maskC(i,j,k-1,bi,bj)
           IF (recip_hFacI(i,j,kp1,bi,bj).EQ.0.) c(i,j,k)=0.  
491          ENDDO          ENDDO
492         ENDDO         ENDDO
493        ENDDO        ENDDO
494  C--   Center diagonal  C--   Center diagonal
495        DO k=1,Nr        DO k=1,Nr
496           km1 = MAX(k-1,1)
497         DO j=jMin,jMax         DO j=jMin,jMax
498          DO i=iMin,iMax          DO i=iMin,iMax
499            b(i,j,k) = 1. _d 0 - c(i,j,k) - a(i,j,k)            b(i,j,k) = 1. _d 0 - c(i,j,k) - a(i,j,k)
500       &        + ab15*deltaTggl90*GGL90ceps*SQRT(GGL90TKE(I,J,K,bi,bj))       &        + implDissFac*deltaTggl90*GGL90ceps*SQRTTKE(i,j,k)
501       &        *rMixingLength(I,J,K)       &        * rMixingLength(i,j,k)
502         &        * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)
503           ENDDO           ENDDO
504         ENDDO         ENDDO
505        ENDDO        ENDDO
506  C     end set up matrix  C     end set up matrix
507    
 C  
508  C     Apply boundary condition  C     Apply boundary condition
509  C        kp1 = MIN(Nr,kSurf+1)
510        DO J=jMin,jMax        DO j=jMin,jMax
511         DO I=iMin,iMax         DO i=iMin,iMax
512  C     estimate friction velocity uStar from surface forcing  C     estimate friction velocity uStar from surface forcing
513          uStarSquare = SQRT(          uStarSquare = SQRT(
514       &         ( .5*( surfaceForcingU(I,  J,  bi,bj)       &    ( .5 _d 0*( surfaceForcingU(i,  j,  bi,bj)
515       &              + surfaceForcingU(I+1,J,  bi,bj) ) )**2       &              + surfaceForcingU(i+1,j,  bi,bj) ) )**2
516       &       + ( .5*( surfaceForcingV(I,  J,  bi,bj)       &  + ( .5 _d 0*( surfaceForcingV(i,  j,  bi,bj)
517       &              + surfaceForcingV(I,  J+1,bi,bj) ) )**2       &              + surfaceForcingV(i,  j+1,bi,bj) ) )**2
518       &                     )       &                     )
519  C     Dirichlet surface boundary condition for TKE  C     Dirichlet surface boundary condition for TKE
520          gTKE(I,J,kSurf) = MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)          GGL90TKE(i,j,kSurf,bi,bj) = maskC(i,j,kSurf,bi,bj)
521       &                     *maskC(I,J,kSurf,bi,bj)       &           *MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)
522            GGL90TKE(i,j,kp1,bi,bj) = GGL90TKE(i,j,kp1,bi,bj)
523         &               - a(i,j,kp1)*GGL90TKE(i,j,kSurf,bi,bj)
524            a(i,j,kp1) = 0. _d 0
525  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom
526          kBottom   = MIN(MAX(kLowC(I,J,bi,bj),1),Nr)          kBottom   = MAX(kLowC(i,j,bi,bj),1)
527          gTKE(I,J,kBottom) = gTKE(I,J,kBottom)          GGL90TKE(i,j,kBottom,bi,bj) = GGL90TKE(i,j,kBottom,bi,bj)
528       &       - GGL90TKEbottom*c(I,J,kBottom)       &                              - GGL90TKEbottom*c(i,j,kBottom)
529         ENDDO          c(i,j,kBottom) = 0. _d 0
530        ENDDO         ENDDO
531  C        ENDDO
532  C     solve tri-diagonal system, and store solution on gTKE (previously rhs)  
533  C  C     solve tri-diagonal system
534        CALL GGL90_SOLVE( bi, bj, iMin, iMax, jMin, jMax,        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
535       I     a, b, c,       I                        a, b, c,
536       U     gTKE,       U                        GGL90TKE,
537       I     myThid )       O                        errCode,
538  C       I                        bi, bj, myThid )
539  C     now update TKE  
540  C        DO k=1,Nr
541        DO K=1,Nr         DO j=jMin,jMax
542         DO J=jMin,jMax          DO i=iMin,iMax
         DO I=iMin,iMax  
543  C     impose minimum TKE to avoid numerical undershoots below zero  C     impose minimum TKE to avoid numerical undershoots below zero
544           GGL90TKE(I,J,K,bi,bj) = MAX( gTKE(I,J,K), GGL90TKEmin )           GGL90TKE(i,j,k,bi,bj) = maskC(i,j,k,bi,bj)
545       &        * maskC(I,J,K,bi,bj)       &                  *MAX( GGL90TKE(i,j,k,bi,bj), GGL90TKEmin )
546          ENDDO          ENDDO
547         ENDDO         ENDDO
548        ENDDO        ENDDO
549  C  
550  C     end of time step  C     end of time step
551  C     ===============================  C     ===============================
552  C     compute viscosity coefficients  
553  C        DO k=2,Nr
554        DO K=2,Nr         DO j=1,sNy
555         DO J=jMin,jMax          DO i=1,sNx
556          DO I=iMin,iMax  #ifdef ALLOW_GGL90_SMOOTH
557  C     Eq. (11), (18) and (21)           tmpVisc=
558           KappaM = GGL90ck*GGL90mixingLength(I,J,K)*       &  (
559       &                  SQRT( GGL90TKE(I,J,K,bi,bj) )       &   p4 *  GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
560           KappaH = KappaM/TKEPrandtlNumber(I,J,K)       &  +p8 *( GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj)
561  C     Set a minium (= background) value       &       + GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj)
562           KappaM = MAX(KappaM,viscAr)       &       + GGL90visctmp(i+1,j  ,k) * mskCor(i+1,j  ,bi,bj)
563           KappaH = MAX(KappaH,diffKrNrT(k))       &       + GGL90visctmp(i  ,j+1,k) * mskCor(i  ,j+1,bi,bj))
564  C     Set a maximum and mask land point       &  +p16*( GGL90visctmp(i+1,j+1,k) * mskCor(i+1,j+1,bi,bj)
565           GGL90viscAr(I,J,K,bi,bj) = MIN(KappaM,GGL90viscMax)       &       + GGL90visctmp(i+1,j-1,k) * mskCor(i+1,j-1,bi,bj)
566       &        * maskC(I,J,K,bi,bj)       &       + GGL90visctmp(i-1,j+1,k) * mskCor(i-1,j+1,bi,bj)
567           GGL90diffKr(I,J,K,bi,bj) = MIN(KappaH,GGL90diffMax)       &       + GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj))
568       &        * maskC(I,J,K,bi,bj)       &  )
569         & /(p4
570         &  +p8 *(       maskC(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)
571         &       +       maskC(i  ,j-1,k,bi,bj) * mskCor(i  ,j-1,bi,bj)
572         &       +       maskC(i+1,j  ,k,bi,bj) * mskCor(i+1,j  ,bi,bj)
573         &       +       maskC(i  ,j+1,k,bi,bj) * mskCor(i  ,j+1,bi,bj))
574         &  +p16*(       maskC(i+1,j+1,k,bi,bj) * mskCor(i+1,j+1,bi,bj)
575         &       +       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,j,k,bi,bj)*mskCor(i,j,bi,bj)
579    #else
580             tmpVisc = GGL90visctmp(i,j,k)
581    #endif
582             tmpVisc = MIN(tmpVisc/TKEPrandtlNumber(i,j,k),GGL90diffMax)
583             GGL90diffKr(i,j,k,bi,bj)= MAX( tmpVisc , diffKrNrT(k) )
584            ENDDO
585           ENDDO
586          ENDDO
587    
588          DO k=2,Nr
589           DO j=1,sNy
590            DO i=1,sNx+1
591    #ifdef ALLOW_GGL90_SMOOTH
592            tmpVisc =
593         & (
594         &   p4 *(GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
595         &       +GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj))
596         &  +p8 *(GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj)
597         &       +GGL90visctmp(i-1,j+1,k) * mskCor(i-1,j+1,bi,bj)
598         &       +GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj)
599         &       +GGL90visctmp(i  ,j+1,k) * mskCor(i  ,j+1,bi,bj))
600         &  )
601         & /(p4 * 2. _d 0
602         &  +p8 *(      maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj)
603         &       +      maskC(i-1,j+1,k,bi,bj) * mskCor(i-1,j+1,bi,bj)
604         &       +      maskC(i  ,j-1,k,bi,bj) * mskCor(i  ,j-1,bi,bj)
605         &       +      maskC(i  ,j+1,k,bi,bj) * mskCor(i  ,j+1,bi,bj))
606         &  )
607         &  *maskC(i  ,j,k,bi,bj)*mskCor(i  ,j,bi,bj)
608         &  *maskC(i-1,j,k,bi,bj)*mskCor(i-1,j,bi,bj)
609    #else
610            tmpVisc = _maskW(i,j,k,bi,bj) *
611         &                   (.5 _d 0*(GGL90visctmp(i,j,k)
612         &                            +GGL90visctmp(i-1,j,k))
613         &                   )
614    #endif
615            tmpVisc = MIN( tmpVisc , GGL90viscMax )
616            GGL90viscArU(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
617          ENDDO          ENDDO
618         ENDDO         ENDDO
 C     end third k-loop  
619        ENDDO        ENDDO
620    
621          DO k=2,Nr
622           DO j=1,sNy+1
623            DO i=1,sNx
624    #ifdef ALLOW_GGL90_SMOOTH
625            tmpVisc =
626         & (
627         &   p4 *(GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
628         &       +GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj))
629         &  +p8 *(GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj)
630         &       +GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj)
631         &       +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         &  )
634         & /(p4 * 2. _d 0
635         &  +p8 *(      maskC(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)
636         &       +      maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj)
637         &       +      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         &  )
640         &   *maskC(i,j  ,k,bi,bj)*mskCor(i,j  ,bi,bj)
641         &   *maskC(i,j-1,k,bi,bj)*mskCor(i,j-1,bi,bj)
642    #else
643            tmpVisc = _maskS(i,j,k,bi,bj) *
644         &                   (.5 _d 0*(GGL90visctmp(i,j,k)
645         &                            +GGL90visctmp(i,j-1,k))
646         &                   )
647    
648    #endif
649            tmpVisc = MIN( tmpVisc , GGL90viscMax )
650            GGL90viscArV(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
651            ENDDO
652           ENDDO
653          ENDDO
654    
655    #ifdef ALLOW_DIAGNOSTICS
656          IF ( useDiagnostics ) THEN
657             CALL DIAGNOSTICS_FILL( GGL90TKE   ,'GGL90TKE',
658         &                          0,Nr, 1, bi, bj, myThid )
659             CALL DIAGNOSTICS_FILL( GGL90viscArU,'GGL90ArU',
660         &                          0,Nr, 1, bi, bj, myThid )
661             CALL DIAGNOSTICS_FILL( GGL90viscArV,'GGL90ArV',
662         &                          0,Nr, 1, bi, bj, myThid )
663             CALL DIAGNOSTICS_FILL( GGL90diffKr,'GGL90Kr ',
664         &                          0,Nr, 1, bi, bj, myThid )
665             CALL DIAGNOSTICS_FILL( TKEPrandtlNumber ,'GGL90Prl',
666         &                          0,Nr, 2, bi, bj, myThid )
667             CALL DIAGNOSTICS_FILL( GGL90mixingLength,'GGL90Lmx',
668         &                          0,Nr, 2, bi, bj, myThid )
669          ENDIF
670    #endif
671    
672  #endif /* ALLOW_GGL90 */  #endif /* ALLOW_GGL90 */
673    
674        RETURN        RETURN
675        END        END
   
Legend:
Removed from v.1.9  
changed lines
  Added in v.1.19
  ViewVC Help
Powered by ViewVC 1.1.22