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