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