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