/[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.14 by jmc, Sat Oct 31 03:18:50 2009 UTC
# Line 27  C global parameters updated by ggl90_cal Line 27  C global parameters updated by ggl90_cal
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: ============================================================
# Line 85  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
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  C-    tri-diagonal matrix
98        _RL     a(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     a(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
99        _RL     b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
100        _RL     c(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     c(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
101    c     INTEGER errCode
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
              rMixingLength(I,J,K) = 0. _d 0  
139          ENDDO          ENDDO
140         ENDDO         ENDDO
141        ENDDO        ENDDO
142        DO J=1-Oly,sNy+Oly        DO J=1-Oly,sNy+Oly
143         DO I=1-Olx,sNx+Olx         DO I=1-Olx,sNx+Olx
144          rhoK    (I,J) = 0. _d 0          rhoK(I,J)          = 0. _d 0
145          rhoKm1  (I,J) = 0. _d 0          rhoKm1(I,J)        = 0. _d 0
146          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)
147            rMixingLength(i,j,1) = 0. _d 0
148            mxLength_Dn(I,J,1) = GGL90mixingLengthMin
149            SQRTTKE(i,j,1) = SQRT( GGL90TKE(i,j,1,bi,bj) )
150         ENDDO         ENDDO
151        ENDDO        ENDDO
152    
# Line 165  c      Kp1 = MIN(Nr,K+1) Line 168  c      Kp1 = MIN(Nr,K+1)
168         DO J=jMin,jMax         DO J=jMin,jMax
169          DO I=iMin,iMax          DO I=iMin,iMax
170           SQRTTKE(i,j,k)=SQRT( GGL90TKE(I,J,K,bi,bj) )           SQRTTKE(i,j,k)=SQRT( GGL90TKE(I,J,K,bi,bj) )
171  C  
172  C     buoyancy frequency  C     buoyancy frequency
 C  
173           Nsquare(i,j,k) = - gravity*recip_rhoConst*recip_drC(K)           Nsquare(i,j,k) = - gravity*recip_rhoConst*recip_drC(K)
174       &        * ( rhoKm1(I,J) - rhoK(I,J) )*maskC(I,J,K,bi,bj)       &        * ( rhoKm1(I,J) - rhoK(I,J) )*maskC(I,J,K,bi,bj)
175  cC     vertical shear term (dU/dz)^2+(dV/dz)^2  cC     vertical shear term (dU/dz)^2+(dV/dz)^2
# 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
205            ENDDO
206           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           ENDDO
215          ENDDO          ENDDO
216         ENDDO         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))
290           rMixingLength(I,J,K) = 1. _d 0 /GGL90mixingLength(I,J,K)          ENDDO
291           ENDDO
292           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          ENDDO
299         ENDDO         ENDDO
300        ENDDO  
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 338  C     ... across y-faces Line 414  C     ... across y-faces
414           ENDDO           ENDDO
415          ENDDO          ENDDO
416  C     Compute divergence of fluxes  C     Compute divergence of fluxes
417    C- jmc: concerned about missing a deltaT since gTKE is already the future TKE.
418          DO j=1-Oly,sNy+Oly-1          DO j=1-Oly,sNy+Oly-1
419           DO i=1-Olx,sNx+Olx-1           DO i=1-Olx,sNx+Olx-1
420            gTKE(i,j,k)=gTKE(i,j,k)            gTKE(i,j,k)=gTKE(i,j,k)
# Line 365  C--   Lower diagonal Line 442  C--   Lower diagonal
442         ENDDO         ENDDO
443        ENDDO        ENDDO
444        DO k=2,Nr        DO k=2,Nr
445    C- jmc: concerned that a(k=2) should always be zero
446    C       and would be better without recip_hFacC
447         km1=max(2,k-1)         km1=max(2,k-1)
448         DO j=jMin,jMax         DO j=jMin,jMax
449          DO i=iMin,iMax          DO i=iMin,iMax
# Line 383  C--   Upper diagonal Line 462  C--   Upper diagonal
462         ENDDO         ENDDO
463        ENDDO        ENDDO
464        DO k=2,Nr        DO k=2,Nr
465    C- jmc: concerned that c(k) from k=kLow to k=Nr should always be zero
466         DO j=jMin,jMax         DO j=jMin,jMax
467          DO i=iMin,iMax          DO i=iMin,iMax
468    C- jmc: this is dangerous since klowC could be zero if land column
469    C       and would be better without recip_hFacC
470            kp1=min(klowC(i,j,bi,bj),k+1)            kp1=min(klowC(i,j,bi,bj),k+1)
471            c(i,j,k) = -deltaTggl90            c(i,j,k) = -deltaTggl90
472  c     &        *recip_drF( k )*recip_hFacI(i,j,k,bi,bj)  c     &        *recip_drF( k )*recip_hFacI(i,j,k,bi,bj)
# Line 399  C--   Center diagonal Line 481  C--   Center diagonal
481         DO j=jMin,jMax         DO j=jMin,jMax
482          DO i=iMin,iMax          DO i=iMin,iMax
483            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)
484       &        + ab15*deltaTggl90*GGL90ceps*SQRT(GGL90TKE(I,J,K,bi,bj))       &        + ab15*deltaTggl90*GGL90ceps*SQRTTKE(I,J,K)
485       &        *rMixingLength(I,J,K)*maskC(i,j,k,bi,bj)       &        *rMixingLength(I,J,K)*maskC(i,j,k,bi,bj)
486           ENDDO           ENDDO
487         ENDDO         ENDDO
488        ENDDO        ENDDO
489  C     end set up matrix  C     end set up matrix
490    
 C  
491  C     Apply boundary condition  C     Apply boundary condition
492  C  C- jmc: concerned about conservation when a or c are changed after computing b
493        DO J=jMin,jMax        DO J=jMin,jMax
494         DO I=iMin,iMax         DO I=iMin,iMax
495  C     estimate friction velocity uStar from surface forcing  C     estimate friction velocity uStar from surface forcing
# Line 419  C     estimate friction velocity uStar f Line 500  C     estimate friction velocity uStar f
500       &              + surfaceForcingV(I,  J+1,bi,bj) ) )**2       &              + surfaceForcingV(I,  J+1,bi,bj) ) )**2
501       &                     )       &                     )
502  C     Dirichlet surface boundary condition for TKE  C     Dirichlet surface boundary condition for TKE
503    C- jmc: would be much better to update the provisional TKE (i.e. gTKE) at k=2
504          gTKE(I,J,kSurf) = MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)          gTKE(I,J,kSurf) = MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)
505       &                     *maskC(I,J,kSurf,bi,bj)       &                     *maskC(I,J,kSurf,bi,bj)
506  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom
# Line 428  C     Dirichlet bottom boundary conditio Line 510  C     Dirichlet bottom boundary conditio
510          c(I,J,kBottom) = 0. _d 0          c(I,J,kBottom) = 0. _d 0
511         ENDDO         ENDDO
512        ENDDO        ENDDO
513  C  
514  C     solve tri-diagonal system, and store solution on gTKE (previously rhs)  C     solve tri-diagonal system, and store solution on gTKE (previously rhs)
515  C  c     CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
516    c    I                        a, b, c,
517    c    U                        gTKE,
518    c    O                        errCode,
519    c    I                        bi, bj, myThid )
520        CALL GGL90_SOLVE( bi, bj, iMin, iMax, jMin, jMax,        CALL GGL90_SOLVE( bi, bj, iMin, iMax, jMin, jMax,
521       I     a, b, c,       I     a, b, c,
522       U     gTKE,       U     gTKE,
523       I     myThid )       I     myThid )
524  C  
525  C     now update TKE  C     now update TKE
 C  
526        DO K=1,Nr        DO K=1,Nr
527         DO J=jMin,jMax         DO J=jMin,jMax
528          DO I=iMin,iMax          DO I=iMin,iMax
Legend:
Removed from v.1.12  
changed lines
  Added in v.1.14
  ViewVC Help
Powered by ViewVC 1.1.22