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