/[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.12 by jmc, Thu Oct 8 20:07:18 2009 UTC revision 1.13 by dfer, Tue Oct 27 15:23:22 2009 UTC
# Line 85  c     _RL     SQRTTKE Line 85  c     _RL     SQRTTKE
85        _RL     RiNumber        _RL     RiNumber
86        _RL     TKEdissipation        _RL     TKEdissipation
87        _RL     tempU, tempV, prTemp        _RL     tempU, tempV, prTemp
88        _RL     MaxLength        _RL     MaxLength, tmpmlx
89        _RL     TKEPrandtlNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     TKEPrandtlNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
90        _RL     GGL90mixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     GGL90mixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
91        _RL         rMixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     rMixingLength    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
92          _RL     mxLength_Dn      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
93        _RL     KappaE           (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     KappaE           (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
94        _RL     rhoK             (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     rhoK             (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
95        _RL     rhoKm1           (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     rhoKm1           (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
96        _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
97        _RL     gTKE             (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     gTKE             (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
98  C     tri-diagonal matrix  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  #ifdef ALLOW_GGL90_HORIZDIFF  #ifdef ALLOW_GGL90_HORIZDIFF
103  C     xA, yA   - area of lateral faces  C-    xA, yA   - area of lateral faces
104  C     dfx, dfy - diffusive flux across lateral faces  C-    dfx, dfy - diffusive flux across lateral faces
105        _RL     xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
106        _RL     yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
107        _RL     dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 133  C     Initialize local fields Line 134  C     Initialize local fields
134          DO I=1-Olx,sNx+Olx          DO I=1-Olx,sNx+Olx
135           gTKE(I,J,K)              = 0. _d 0           gTKE(I,J,K)              = 0. _d 0
136           KappaE(I,J,K)            = 0. _d 0           KappaE(I,J,K)            = 0. _d 0
137           TKEPrandtlNumber(I,J,K)  = 0. _d 0           TKEPrandtlNumber(I,J,K)  = 1. _d 0
138           GGL90mixingLength(I,J,K) = GGL90mixingLengthMin           GGL90mixingLength(I,J,K) = GGL90mixingLengthMin
139               rMixingLength(I,J,K) = 0. _d 0               rMixingLength(I,J,K) = 0. _d 0
140          ENDDO          ENDDO
# Line 141  C     Initialize local fields Line 142  C     Initialize local fields
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            mxLength_Dn(I,J,1) = GGL90mixingLengthMin
149         ENDDO         ENDDO
150        ENDDO        ENDDO
151    
# Line 190  C     mixing length Line 192  C     mixing length
192         ENDDO         ENDDO
193        ENDDO        ENDDO
194    
195  C- Impose upper bound for mixing length (total depth)  C- Impose upper and lower bound for mixing length
196        IF ( mxlMaxFlag .EQ. 0 ) THEN        IF ( mxlMaxFlag .EQ. 0 ) THEN
197    C-
198         DO k=2,Nr         DO k=2,Nr
199          DO J=jMin,jMax          DO J=jMin,jMax
200           DO I=iMin,iMax           DO I=iMin,iMax
201            MaxLength=totalDepth(I,J)            MaxLength=totalDepth(I,J)
202            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),
203       &        MaxLength)       &                                   MaxLength)
204           ENDDO           ENDDO
205          ENDDO          ENDDO
206         ENDDO         ENDDO
207    
208           DO k=2,Nr
209            DO J=jMin,jMax
210             DO I=iMin,iMax
211              GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),
212         &                                   GGL90mixingLengthMin)
213              rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K)
214             ENDDO
215            ENDDO
216           ENDDO
217    
218        ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN
219    C-
220         DO k=2,Nr         DO k=2,Nr
221          DO J=jMin,jMax          DO J=jMin,jMax
222           DO I=iMin,iMax           DO I=iMin,iMax
223            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))
224  c         MaxLength=MAX(MaxLength,20. _d 0)  c         MaxLength=MAX(MaxLength,20. _d 0)
225            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),
226       &        MaxLength)       &                                   MaxLength)
227             ENDDO
228            ENDDO
229           ENDDO
230    
231           DO k=2,Nr
232            DO J=jMin,jMax
233             DO I=iMin,iMax
234              GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),
235         &                                   GGL90mixingLengthMin)
236              rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K)
237           ENDDO           ENDDO
238          ENDDO          ENDDO
239         ENDDO         ENDDO
240    
241        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN
242    C-
243         DO k=2,Nr         DO k=2,Nr
244          DO J=jMin,jMax          DO J=jMin,jMax
245           DO I=iMin,iMax           DO I=iMin,iMax
# Line 235  c         MaxLength=MAX(MaxLength,20. _d Line 262  c         MaxLength=MAX(MaxLength,20. _d
262           ENDDO           ENDDO
263          ENDDO          ENDDO
264         ENDDO         ENDDO
       ELSE  
        STOP 'GGL90_CALC: Wrong mxlMaxFlag (mixing lenght limit)'  
       ENDIF  
265    
266  C- Impose minimum mixing length (to avoid division by zero)         DO k=2,Nr
267        DO k=2,Nr          DO J=jMin,jMax
268             DO I=iMin,iMax
269              GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),
270         &                                   GGL90mixingLengthMin)
271              rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K)
272             ENDDO
273            ENDDO
274           ENDDO
275    
276          ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN
277    C-
278           DO k=2,Nr
279            DO J=jMin,jMax
280             DO I=iMin,iMax
281              mxLength_Dn(I,J,K) = MIN(GGL90mixingLength(I,J,K),
282         &        mxLength_Dn(I,J,K-1)+drF(k-1))
283             ENDDO
284            ENDDO
285           ENDDO
286         DO J=jMin,jMax         DO J=jMin,jMax
287          DO I=iMin,iMax          DO I=iMin,iMax
288           GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),            GGL90mixingLength(I,J,Nr) = MIN(GGL90mixingLength(I,J,Nr),
289       &        GGL90mixingLengthMin)       &       GGL90mixingLengthMin+drF(Nr))
          rMixingLength(I,J,K) = 1. _d 0 /GGL90mixingLength(I,J,K)  
290          ENDDO          ENDDO
291         ENDDO         ENDDO
292        ENDDO         DO k=Nr-1,2,-1
293            DO J=jMin,jMax
294             DO I=iMin,iMax
295              GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),
296         &        GGL90mixingLength(I,J,K+1)+drF(k))
297             ENDDO
298            ENDDO
299           ENDDO
300    
301           DO k=2,Nr
302            DO J=jMin,jMax
303             DO I=iMin,iMax
304              GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),
305         &                                  mxLength_Dn(I,J,K))
306              tmpmlx = SQRT(GGL90mixingLength(I,J,K) * mxLength_Dn(I,J,K))
307              tmpmlx = MAX( tmpmlx, GGL90mixingLengthMin)
308              rMixingLength(I,J,K) = 1. _d 0 / tmpmlx
309             ENDDO
310            ENDDO
311           ENDDO
312    
313          ELSE
314           STOP 'GGL90_CALC: Wrong mxlMaxFlag (mixing length limit)'
315          ENDIF
316    
317    C- Impose minimum mixing length (to avoid division by zero)
318    c      DO k=2,Nr
319    c      DO J=jMin,jMax
320    c       DO I=iMin,iMax
321    c        GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),
322    c    &        GGL90mixingLengthMin)
323    c        rMixingLength(I,J,K) = 1. _d 0 /GGL90mixingLength(I,J,K)
324    c       ENDDO
325    c      ENDDO
326    c     ENDDO
327    
328        DO k=2,Nr        DO k=2,Nr
329         Km1 = K-1         Km1 = K-1
# Line 267  C     compute Prandtl number (always gre Line 342  C     compute Prandtl number (always gre
342           prTemp = 1. _d 0           prTemp = 1. _d 0
343           IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber           IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
344           TKEPrandtlNumber(I,J,K) = MIN(10. _d 0,prTemp)           TKEPrandtlNumber(I,J,K) = MIN(10. _d 0,prTemp)
345    c         TKEPrandtlNumber(I,J,K) = 1. _d 0
346    
347  C     viscosity and diffusivity  C     viscosity and diffusivity
348           KappaM = GGL90ck*GGL90mixingLength(I,J,K)*SQRTTKE(i,j,k)           KappaM = GGL90ck*GGL90mixingLength(I,J,K)*SQRTTKE(i,j,k)
# Line 399  C--   Center diagonal Line 475  C--   Center diagonal
475         DO j=jMin,jMax         DO j=jMin,jMax
476          DO i=iMin,iMax          DO i=iMin,iMax
477            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)
478       &        + ab15*deltaTggl90*GGL90ceps*SQRT(GGL90TKE(I,J,K,bi,bj))       &        + ab15*deltaTggl90*GGL90ceps*SQRTTKE(I,J,K)
479       &        *rMixingLength(I,J,K)*maskC(i,j,k,bi,bj)       &        *rMixingLength(I,J,K)*maskC(i,j,k,bi,bj)
480           ENDDO           ENDDO
481         ENDDO         ENDDO
Legend:
Removed from v.1.12  
changed lines
  Added in v.1.13
  ViewVC Help
Powered by ViewVC 1.1.22