/[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.11 by dfer, Fri Jan 30 02:23:56 2009 UTC revision 1.16 by gforget, Fri Aug 6 18:37:05 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
# Line 86  c     _RL     SQRTTKE Line 84  c     _RL     SQRTTKE
84        _RL     RiNumber        _RL     RiNumber
85        _RL     TKEdissipation        _RL     TKEdissipation
86        _RL     tempU, tempV, prTemp        _RL     tempU, tempV, prTemp
87        _RL     MaxLength        _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  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
112        _RL p4, p8, p16, tmpdiffKrS        _RL p4, p8, p16
113        p4=0.25 _d 0        p4=0.25 _d 0
114        p8=0.125 _d 0        p8=0.125 _d 0
115        p16=0.0625 _d 0        p16=0.0625 _d 0
116  #endif  #endif
 CEOP  
117        iMin = 2-OLx        iMin = 2-OLx
118        iMax = sNx+OLx-1        iMax = sNx+OLx-1
119        jMin = 2-OLy        jMin = 2-OLy
# Line 121  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 135  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
              rMixingLength(I,J,K) = 0. _d 0  
140          ENDDO          ENDDO
141         ENDDO         ENDDO
142        ENDDO        ENDDO
143        DO J=1-Oly,sNy+Oly        DO J=1-Oly,sNy+Oly
144         DO I=1-Olx,sNx+Olx         DO I=1-Olx,sNx+Olx
145          rhoK    (I,J) = 0. _d 0          rhoK(I,J)          = 0. _d 0
146          rhoKm1  (I,J) = 0. _d 0          rhoKm1(I,J)        = 0. _d 0
147          totalDepth(I,J) = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)          totalDepth(I,J)    = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
148            rMixingLength(i,j,1) = 0. _d 0
149            mxLength_Dn(I,J,1) = GGL90mixingLengthMin
150            SQRTTKE(i,j,1) = SQRT( GGL90TKE(i,j,1,bi,bj) )
151         ENDDO         ENDDO
152        ENDDO        ENDDO
153    
# Line 167  c      Kp1 = MIN(Nr,K+1) Line 169  c      Kp1 = MIN(Nr,K+1)
169         DO J=jMin,jMax         DO J=jMin,jMax
170          DO I=iMin,iMax          DO I=iMin,iMax
171           SQRTTKE(i,j,k)=SQRT( GGL90TKE(I,J,K,bi,bj) )           SQRTTKE(i,j,k)=SQRT( GGL90TKE(I,J,K,bi,bj) )
172  C  
173  C     buoyancy frequency  C     buoyancy frequency
 C  
174           Nsquare(i,j,k) = - gravity*recip_rhoConst*recip_drC(K)           Nsquare(i,j,k) = - gravity*recip_rhoConst*recip_drC(K)
175       &        * ( rhoKm1(I,J) - rhoK(I,J) )*maskC(I,J,K,bi,bj)       &        * ( rhoKm1(I,J) - rhoK(I,J) )*maskC(I,J,K,bi,bj)
176  cC     vertical shear term (dU/dz)^2+(dV/dz)^2  cC     vertical shear term (dU/dz)^2+(dV/dz)^2
# Line 192  C     mixing length Line 193  C     mixing length
193         ENDDO         ENDDO
194        ENDDO        ENDDO
195    
196  C- Impose upper bound for mixing length (total depth)  C- Impose upper and lower bound for mixing length
197        IF ( mxlMaxFlag .EQ. 0 ) THEN        IF ( mxlMaxFlag .EQ. 0 ) THEN
198    C-
199         DO k=2,Nr         DO k=2,Nr
200          DO J=jMin,jMax          DO J=jMin,jMax
201           DO I=iMin,iMax           DO I=iMin,iMax
202            MaxLength=totalDepth(I,J)            MaxLength=totalDepth(I,J)
203            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),
204       &        MaxLength)       &                                   MaxLength)
205             ENDDO
206            ENDDO
207           ENDDO
208    
209           DO k=2,Nr
210            DO J=jMin,jMax
211             DO I=iMin,iMax
212              GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),
213         &                                   GGL90mixingLengthMin)
214              rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K)
215           ENDDO           ENDDO
216          ENDDO          ENDDO
217         ENDDO         ENDDO
218    
219        ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN
220    C-
221         DO k=2,Nr         DO k=2,Nr
222          DO J=jMin,jMax          DO J=jMin,jMax
223           DO I=iMin,iMax           DO I=iMin,iMax
224            MaxLength=MIN(Ro_surf(I,J,bi,bj)-rF(k),rF(k)-R_low(I,J,bi,bj))            MaxLength=MIN(Ro_surf(I,J,bi,bj)-rF(k),rF(k)-R_low(I,J,bi,bj))
225  c         MaxLength=MAX(MaxLength,20. _d 0)  c         MaxLength=MAX(MaxLength,20. _d 0)
226            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),
227       &        MaxLength)       &                                   MaxLength)
228           ENDDO           ENDDO
229          ENDDO          ENDDO
230         ENDDO         ENDDO
231    
232           DO k=2,Nr
233            DO J=jMin,jMax
234             DO I=iMin,iMax
235              GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),
236         &                                   GGL90mixingLengthMin)
237              rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K)
238             ENDDO
239            ENDDO
240           ENDDO
241    
242        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN
243    C-
244    cgf ensure mixing between first and second level
245    c      DO J=jMin,jMax
246    c        DO I=iMin,iMax
247    c         GGL90mixingLength(I,J,2)=drF(1)
248    c        ENDDO
249    c      ENDDO
250    cgf
251         DO k=2,Nr         DO k=2,Nr
252          DO J=jMin,jMax          DO J=jMin,jMax
253           DO I=iMin,iMax           DO I=iMin,iMax
# Line 234  c         MaxLength=MAX(MaxLength,20. _d Line 267  c         MaxLength=MAX(MaxLength,20. _d
267           DO I=iMin,iMax           DO I=iMin,iMax
268            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),
269       &        GGL90mixingLength(I,J,K+1)+drF(k))       &        GGL90mixingLength(I,J,K+1)+drF(k))
270           ENDDO           ENDDO
271          ENDDO            ENDDO
272         ENDDO         ENDDO
       ELSE  
        STOP 'GGL90_CALC: Wrong mxlMaxFlag (mixing lenght limit)'  
       ENDIF  
273    
274  C- Impose minimum mixing length (to avoid division by zero)         DO k=2,Nr
275        DO k=2,Nr          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. 3 ) THEN
285    C-
286           DO k=2,Nr
287            DO J=jMin,jMax
288             DO I=iMin,iMax
289              mxLength_Dn(I,J,K) = MIN(GGL90mixingLength(I,J,K),
290         &        mxLength_Dn(I,J,K-1)+drF(k-1))
291             ENDDO
292            ENDDO
293           ENDDO
294         DO J=jMin,jMax         DO J=jMin,jMax
295          DO I=iMin,iMax          DO I=iMin,iMax
296           GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),            GGL90mixingLength(I,J,Nr) = MIN(GGL90mixingLength(I,J,Nr),
297       &        GGL90mixingLengthMin)       &       GGL90mixingLengthMin+drF(Nr))
          rMixingLength(I,J,K) = 1. _d 0 /GGL90mixingLength(I,J,K)  
298          ENDDO          ENDDO
299         ENDDO         ENDDO
300        ENDDO         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) = MIN(GGL90mixingLength(I,J,K),
313         &                                  mxLength_Dn(I,J,K))
314              tmpmlx = SQRT( GGL90mixingLength(I,J,K)*mxLength_Dn(I,J,K) )
315              tmpmlx = MAX( tmpmlx, GGL90mixingLengthMin)
316              rMixingLength(I,J,K) = 1. _d 0 / tmpmlx
317             ENDDO
318            ENDDO
319           ENDDO
320    
321          ELSE
322           STOP 'GGL90_CALC: Wrong mxlMaxFlag (mixing length limit)'
323          ENDIF
324    
325    C- Impose minimum mixing length (to avoid division by zero)
326    c      DO k=2,Nr
327    c      DO J=jMin,jMax
328    c       DO I=iMin,iMax
329    c        GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),
330    c    &        GGL90mixingLengthMin)
331    c        rMixingLength(I,J,K) = 1. _d 0 /GGL90mixingLength(I,J,K)
332    c       ENDDO
333    c      ENDDO
334    c     ENDDO
335    
336    
337        DO k=2,Nr        DO k=2,Nr
338         Km1 = K-1         Km1 = K-1
# Line 264  C     vertical shear term (dU/dz)^2+(dV/ Line 346  C     vertical shear term (dU/dz)^2+(dV/
346       &                 -( 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)) )
347       &        *recip_drC(K)       &        *recip_drC(K)
348           verticalShear = tempU*tempU + tempV*tempV           verticalShear = tempU*tempU + tempV*tempV
349           RiNumber   = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)           RiNumber = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)
350  C     compute Prandtl number (always greater than 0)  C     compute Prandtl number (always greater than 0)
351           prTemp = 1. _d 0           prTemp = 1. _d 0
352           IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber           IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
353           TKEPrandtlNumber(I,J,K) = MIN(10. _d 0,prTemp)           TKEPrandtlNumber(I,J,K) = MIN(10. _d 0,prTemp)
354    c         TKEPrandtlNumber(I,J,K) = 1. _d 0
355    
356  C     viscosity and diffusivity  C     viscosity and diffusivity
357           KappaM = GGL90ck*GGL90mixingLength(I,J,K)*SQRTTKE(i,j,k)           KappaM = GGL90ck*GGL90mixingLength(I,J,K)*SQRTTKE(i,j,k)
358             GGL90visctmp(I,J,K) = MAX(KappaM,diffKrNrT(k))
359         &                            * maskC(I,J,K,bi,bj)
360    c        note: storing GGL90visctmp like this, and using it later to compute
361    c              GGL9rdiffKr etc. is robust in case of smoothing (e.g. see OPA)
362             KappaM = MAX(KappaM,viscArNr(k)) * maskC(I,J,K,bi,bj)
363           KappaH = KappaM/TKEPrandtlNumber(I,J,K)           KappaH = KappaM/TKEPrandtlNumber(I,J,K)
364             KappaE(I,J,K) = GGL90alpha * KappaM * maskC(I,J,K,bi,bj)
 C     Set a minium (= background) and maximum value  
          KappaM = MAX(KappaM,viscAr)  
          KappaH = MAX(KappaH,diffKrNrT(k))  
          KappaM = MIN(KappaM,GGL90viscMax)  
          KappaH = MIN(KappaH,GGL90diffMax)  
   
 C     Mask land points and storage  
          GGL90viscAr(I,J,K,bi,bj) = KappaM * maskC(I,J,K,bi,bj)  
          GGL90diffKr(I,J,K,bi,bj) = KappaH * maskC(I,J,K,bi,bj)  
          KappaE(I,J,K) = GGL90alpha * GGL90viscAr(I,J,K,bi,bj)  
365    
366  C     dissipation term  C     dissipation term
367           TKEdissipation = ab05*GGL90ceps           TKEdissipation = ab05*GGL90ceps
# Line 346  C     Compute divergence of fluxes Line 424  C     Compute divergence of fluxes
424       &   -_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)
425       &         *( (dfx(i+1,j)-dfx(i,j))       &         *( (dfx(i+1,j)-dfx(i,j))
426       &           +(dfy(i,j+1)-dfy(i,j))       &           +(dfy(i,j+1)-dfy(i,j))
427       &           )       &           )*deltaTggl90
428           ENDDO           ENDDO
429          ENDDO          ENDDO
430  C       end of k-loop  C       end of k-loop
# Line 367  C--   Lower diagonal Line 445  C--   Lower diagonal
445         ENDDO         ENDDO
446        ENDDO        ENDDO
447        DO k=2,Nr        DO k=2,Nr
448         km1=max(2,k-1)         km1=MAX(2,k-1)
449         DO j=jMin,jMax         DO j=jMin,jMax
450          DO i=iMin,iMax          DO i=iMin,iMax
451    C-    We keep recip_hFacC in the diffusive flux calculation,
452    C-    but no hFacC in TKE volume control
453    C-    No need for maskC(k-1) with recip_hFacC(k-1)
454           a(i,j,k) = -deltaTggl90           a(i,j,k) = -deltaTggl90
 c     &        *recip_drF(km1)*recip_hFacI(i,j,k,bi,bj)  
455       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)
456       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))
457       &        *recip_drC(k)*maskC(i,j,k,bi,bj)*maskC(i,j,k-1,bi,bj)       &        *recip_drC(k)*maskC(i,j,k,bi,bj)
458          ENDDO          ENDDO
459         ENDDO         ENDDO
460        ENDDO        ENDDO
# Line 387  C--   Upper diagonal Line 467  C--   Upper diagonal
467        DO k=2,Nr        DO k=2,Nr
468         DO j=jMin,jMax         DO j=jMin,jMax
469          DO i=iMin,iMax          DO i=iMin,iMax
470            kp1=min(klowC(i,j,bi,bj),k+1)            kp1=MAX(1,MIN(klowC(i,j,bi,bj),k+1))
471    C-    We keep recip_hFacC in the diffusive flux calculation,
472    C-    but no hFacC in TKE volume control
473    C-    No need for maskC(k) with recip_hFacC(k)
474            c(i,j,k) = -deltaTggl90            c(i,j,k) = -deltaTggl90
 c     &        *recip_drF( k )*recip_hFacI(i,j,k,bi,bj)  
475       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)
476       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))
477       &        *recip_drC(k)*maskC(i,j,k,bi,bj)*maskC(i,j,k-1,bi,bj)       &        *recip_drC(k)*maskC(i,j,k-1,bi,bj)
478          ENDDO          ENDDO
479         ENDDO         ENDDO
480        ENDDO        ENDDO
481  C--   Center diagonal  C--   Center diagonal
482        DO k=1,Nr        DO k=1,Nr
483           km1 = MAX(k-1,1)
484         DO j=jMin,jMax         DO j=jMin,jMax
485          DO i=iMin,iMax          DO i=iMin,iMax
486            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)
487       &        + ab15*deltaTggl90*GGL90ceps*SQRT(GGL90TKE(I,J,K,bi,bj))       &        + ab15*deltaTggl90*GGL90ceps*SQRTTKE(I,J,K)
488       &        *rMixingLength(I,J,K)*maskC(i,j,k,bi,bj)       &        * rMixingLength(I,J,K)
489         &        * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)
490           ENDDO           ENDDO
491         ENDDO         ENDDO
492        ENDDO        ENDDO
493  C     end set up matrix  C     end set up matrix
494    
 C  
495  C     Apply boundary condition  C     Apply boundary condition
496  C        kp1 = MIN(Nr,kSurf+1)
497        DO J=jMin,jMax        DO J=jMin,jMax
498         DO I=iMin,iMax         DO I=iMin,iMax
499  C     estimate friction velocity uStar from surface forcing  C     estimate friction velocity uStar from surface forcing
# Line 423  C     estimate friction velocity uStar f Line 506  C     estimate friction velocity uStar f
506  C     Dirichlet surface boundary condition for TKE  C     Dirichlet surface boundary condition for TKE
507          gTKE(I,J,kSurf) = MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)          gTKE(I,J,kSurf) = MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)
508       &                     *maskC(I,J,kSurf,bi,bj)       &                     *maskC(I,J,kSurf,bi,bj)
509            gTKE(i,j,kp1) = gTKE(i,j,kp1)
510         &                - a(i,j,kp1)*gTKE(i,j,kSurf)
511            a(i,j,kp1) = 0. _d 0
512  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom
513          kBottom   = MAX(kLowC(I,J,bi,bj),1)          kBottom   = MAX(kLowC(I,J,bi,bj),1)
514          gTKE(I,J,kBottom) = gTKE(I,J,kBottom)          gTKE(I,J,kBottom) = gTKE(I,J,kBottom)
# Line 430  C     Dirichlet bottom boundary conditio Line 516  C     Dirichlet bottom boundary conditio
516          c(I,J,kBottom) = 0. _d 0          c(I,J,kBottom) = 0. _d 0
517         ENDDO         ENDDO
518        ENDDO        ENDDO
519  C  
520  C     solve tri-diagonal system, and store solution on gTKE (previously rhs)  C     solve tri-diagonal system, and store solution on gTKE (previously rhs)
521  C        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
522        CALL GGL90_SOLVE( bi, bj, iMin, iMax, jMin, jMax,       I                        a, b, c,
523       I     a, b, c,       U                        gTKE,
524       U     gTKE,       O                        errCode,
525       I     myThid )       I                        bi, bj, myThid )
526  C  
527  C     now update TKE  C     now update TKE
 C  
528        DO K=1,Nr        DO K=1,Nr
529         DO J=jMin,jMax         DO J=jMin,jMax
530          DO I=iMin,iMax          DO I=iMin,iMax
# Line 453  C     impose minimum TKE to avoid numeri Line 538  C     impose minimum TKE to avoid numeri
538  C     end of time step  C     end of time step
539  C     ===============================  C     ===============================
540    
 #ifdef ALLOW_GGL90_SMOOTH  
541        DO K=1,Nr        DO K=1,Nr
542         DO J=jMin,jMax         DO J=jMin,jMax
543          DO I=iMin,iMax          DO I=iMin,iMax
544           tmpdiffKrS=  #ifdef ALLOW_GGL90_SMOOTH
545             tmpVisc=
546       &  (       &  (
547       &   p4 *  GGL90viscAr(i  ,j  ,k,bi,bj) * mskCor(i  ,j  ,bi,bj)       &   p4 *  GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
548       &  +p8 *( GGL90viscAr(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)       &  +p8 *( GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj)
549       &       + GGL90viscAr(i  ,j-1,k,bi,bj) * mskCor(i  ,j-1,bi,bj)       &       + GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj)
550       &       + GGL90viscAr(i+1,j  ,k,bi,bj) * mskCor(i+1,j  ,bi,bj)       &       + GGL90visctmp(i+1,j  ,k) * mskCor(i+1,j  ,bi,bj)
551       &       + GGL90viscAr(i  ,j+1,k,bi,bj) * mskCor(i  ,j+1,bi,bj))       &       + GGL90visctmp(i  ,j+1,k) * mskCor(i  ,j+1,bi,bj))
552       &  +p16*( GGL90viscAr(i+1,j+1,k,bi,bj) * mskCor(i+1,j+1,bi,bj)       &  +p16*( GGL90visctmp(i+1,j+1,k) * mskCor(i+1,j+1,bi,bj)
553       &       + GGL90viscAr(i+1,j-1,k,bi,bj) * mskCor(i+1,j-1,bi,bj)       &       + GGL90visctmp(i+1,j-1,k) * mskCor(i+1,j-1,bi,bj)
554       &       + GGL90viscAr(i-1,j+1,k,bi,bj) * mskCor(i-1,j+1,bi,bj)       &       + GGL90visctmp(i-1,j+1,k) * mskCor(i-1,j+1,bi,bj)
555       &       + GGL90viscAr(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj))       &       + GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj))
556       &  )       &  )
557       & /(p4       & /(p4
558       &  +p8 *(       maskC(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)       &  +p8 *(       maskC(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)
# Line 479  C     =============================== Line 564  C     ===============================
564       &       +       maskC(i-1,j+1,k,bi,bj) * mskCor(i-1,j+1,bi,bj)       &       +       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))       &       +       maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj))
566       &  )*maskC(i,j,k,bi,bj)*mskCor(i,j,bi,bj)       &  )*maskC(i,j,k,bi,bj)*mskCor(i,j,bi,bj)
567       &   /TKEPrandtlNumber(i,j,k)  #else
568           GGL90diffKrS(I,J,K,bi,bj)= MAX( tmpdiffKrS , diffKrNrT(k) )           tmpVisc = GGL90visctmp(I,J,K)
569    #endif
570             tmpVisc = MIN(tmpVisc/TKEPrandtlNumber(i,j,k),GGL90diffMax)
571             GGL90diffKr(I,J,K,bi,bj)= MAX( tmpVisc , diffKrNrT(k) )
572          ENDDO          ENDDO
573         ENDDO         ENDDO
574        ENDDO        ENDDO
575    
576    
577    
578          DO K=1,Nr
579           DO J=jMin,jMax
580            DO I=iMin,iMax
581    #ifdef ALLOW_GGL90_SMOOTH
582            tmpVisc =
583         & (
584         &   p4 *(GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
585         &       +GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj))
586         &  +p8 *(GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj)
587         &       +GGL90visctmp(i-1,j+1,k) * mskCor(i-1,j+1,bi,bj)
588         &       +GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj)
589         &       +GGL90visctmp(i  ,j+1,k) * mskCor(i  ,j+1,bi,bj))
590         &  )
591         & /(p4 * 2. _d 0
592         &  +p8 *(      maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj)
593         &       +      maskC(i-1,j+1,k,bi,bj) * mskCor(i-1,j+1,bi,bj)
594         &       +      maskC(i  ,j-1,k,bi,bj) * mskCor(i  ,j-1,bi,bj)
595         &       +      maskC(i  ,j+1,k,bi,bj) * mskCor(i  ,j+1,bi,bj))
596         &  )
597         &  *maskC(i  ,j,k,bi,bj)*mskCor(i  ,j,bi,bj)
598         &  *maskC(i-1,j,k,bi,bj)*mskCor(i-1,j,bi,bj)
599    #else
600            tmpVisc = _maskW(i,j,k,bi,bj) *
601         &                   (.5 _d 0*(GGL90visctmp(i,j,k)
602         &                            +GGL90visctmp(i-1,j,k))
603         &                   )
604  #endif  #endif
605            tmpVisc = MIN( tmpVisc , GGL90viscMax )
606            GGL90viscArU(i,j,k,bi,bj) = MAX( tmpVisc , viscArNr(k)  )
607            ENDDO
608           ENDDO
609          ENDDO
610    
611    
612          DO K=1,Nr
613           DO J=jMin,jMax
614            DO I=iMin,iMax
615    #ifdef ALLOW_GGL90_SMOOTH
616            tmpVisc =
617         & (
618         &   p4 *(GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
619         &       +GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj))
620         &  +p8 *(GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj)
621         &       +GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj)
622         &       +GGL90visctmp(i+1,j  ,k) * mskCor(i+1,j  ,bi,bj)
623         &       +GGL90visctmp(i+1,j-1,k) * mskCor(i+1,j-1,bi,bj))
624         &  )
625         & /(p4 * 2. _d 0
626         &  +p8 *(      maskC(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)
627         &       +      maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj)
628         &       +      maskC(i+1,j  ,k,bi,bj) * mskCor(i+1,j  ,bi,bj)
629         &       +      maskC(i+1,j-1,k,bi,bj) * mskCor(i+1,j-1,bi,bj))
630         &  )
631         &   *maskC(i,j  ,k,bi,bj)*mskCor(i,j  ,bi,bj)
632         &   *maskC(i,j-1,k,bi,bj)*mskCor(i,j-1,bi,bj)
633    #else
634            tmpVisc = _maskS(i,j,k,bi,bj) *
635         &                   (.5 _d 0*(GGL90visctmp(i,j,k)
636         &                            +GGL90visctmp(i,j-1,k))
637         &                   )
638    
639    #endif
640            tmpVisc = MIN( tmpVisc , GGL90viscMax )
641            GGL90viscArV(i,j,k,bi,bj) = MAX( tmpVisc , viscArNr(k)  )
642            ENDDO
643           ENDDO
644          ENDDO
645    
646  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
647        IF ( useDiagnostics ) THEN        IF ( useDiagnostics ) THEN
648           CALL DIAGNOSTICS_FILL( GGL90TKE   ,'GGL90TKE',           CALL DIAGNOSTICS_FILL( GGL90TKE   ,'GGL90TKE',
649       &                          0,Nr, 1, bi, bj, myThid )       &                          0,Nr, 1, bi, bj, myThid )
650           CALL DIAGNOSTICS_FILL( GGL90viscAr,'GGL90Ar ',           CALL DIAGNOSTICS_FILL( GGL90viscArU,'GGL90ArU',
651         &                          0,Nr, 1, bi, bj, myThid )
652             CALL DIAGNOSTICS_FILL( GGL90viscArV,'GGL90ArV',
653       &                          0,Nr, 1, bi, bj, myThid )       &                          0,Nr, 1, bi, bj, myThid )
654           CALL DIAGNOSTICS_FILL( GGL90diffKr,'GGL90Kr ',           CALL DIAGNOSTICS_FILL( GGL90diffKr,'GGL90Kr ',
655       &                          0,Nr, 1, bi, bj, myThid )       &                          0,Nr, 1, bi, bj, myThid )
# Line 505  C     =============================== Line 664  C     ===============================
664    
665        RETURN        RETURN
666        END        END
   
Legend:
Removed from v.1.11  
changed lines
  Added in v.1.16
  ViewVC Help
Powered by ViewVC 1.1.22