/[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.13 by dfer, Tue Oct 27 15:23:22 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 99  C-    tri-diagonal matrix Line 98  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
# Line 136  C     Initialize local fields Line 136  C     Initialize local fields
136           KappaE(I,J,K)            = 0. _d 0           KappaE(I,J,K)            = 0. _d 0
137           TKEPrandtlNumber(I,J,K)  = 1. _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
# Line 145  C     Initialize local fields Line 144  C     Initialize local fields
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          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 167  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 303  C- Line 303  C-
303           DO I=iMin,iMax           DO I=iMin,iMax
304            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),
305       &                                  mxLength_Dn(I,J,K))       &                                  mxLength_Dn(I,J,K))
306            tmpmlx = SQRT(GGL90mixingLength(I,J,K) * mxLength_Dn(I,J,K))            tmpmlx = SQRT( GGL90mixingLength(I,J,K)*mxLength_Dn(I,J,K) )
307            tmpmlx = MAX( tmpmlx, GGL90mixingLengthMin)            tmpmlx = MAX( tmpmlx, GGL90mixingLengthMin)
308            rMixingLength(I,J,K) = 1. _d 0 / tmpmlx            rMixingLength(I,J,K) = 1. _d 0 / tmpmlx
309           ENDDO           ENDDO
# Line 414  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 441  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 459  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 482  C--   Center diagonal Line 488  C--   Center diagonal
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 495  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 504  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.13  
changed lines
  Added in v.1.14
  ViewVC Help
Powered by ViewVC 1.1.22