/[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.2 by mlosch, Mon Sep 27 08:02:04 2004 UTC revision 1.7 by mlosch, Tue Jun 6 22:15:19 2006 UTC
# Line 68  C     KappaE           - (local) diffusi Line 68  C     KappaE           - (local) diffusi
68  C     buoyFreq         - buoyancy freqency  C     buoyFreq         - buoyancy freqency
69  C     TKEdissipation   - dissipation of TKE  C     TKEdissipation   - dissipation of TKE
70  C     GGL90mixingLength- mixing length of scheme following Banke+Delecuse  C     GGL90mixingLength- mixing length of scheme following Banke+Delecuse
71    C         rMixingLength- inverse of mixing length
72  C     totalDepth       - thickness of water column (inverse of recip_Rcol)  C     totalDepth       - thickness of water column (inverse of recip_Rcol)
73  C     TKEPrandtlNumber - here, an empirical function of the Richardson number  C     TKEPrandtlNumber - here, an empirical function of the Richardson number
74  C     rhoK, rhoKm1     - density at layer K and Km1 (relative to K)  C     rhoK, rhoKm1     - density at layer K and Km1 (relative to K)
# Line 86  C     gTKE             - right hand side Line 87  C     gTKE             - right hand side
87        _RL     tempU, tempV, prTemp        _RL     tempU, tempV, prTemp
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)
91        _RL     KappaE           (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     KappaE           (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
92        _RL     rhoK             (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     rhoK             (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
93        _RL     rhoKm1           (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     rhoKm1           (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 110  CEOP Line 112  CEOP
112        jMax = sNy+OLy-1        jMax = sNy+OLy-1
113    
114  C     set separate time step (should be deltaTtracer)  C     set separate time step (should be deltaTtracer)
115        deltaTggl90 = deltaTtracer        deltaTggl90 = dTtracerLev(1)
116  C      C    
117        kSurf = 1        kSurf = 1
118  C     implicit timestepping weights for dissipation  C     implicit timestepping weights for dissipation
# Line 126  C     Initialize local fields Line 128  C     Initialize local fields
128           gTKE(I,J,K)              = 0. _d 0           gTKE(I,J,K)              = 0. _d 0
129           KappaE(I,J,K)            = 0. _d 0           KappaE(I,J,K)            = 0. _d 0
130           TKEPrandtlNumber(I,J,K)  = 0. _d 0           TKEPrandtlNumber(I,J,K)  = 0. _d 0
131           GGL90mixingLength(I,J,K) = 0. _d 0           GGL90mixingLength(I,J,K) = GGL90mixingLengthMin
132                 rMixingLength(I,J,K) = 0. _d 0
133          ENDDO          ENDDO
134         ENDDO             ENDDO    
135        ENDDO        ENDDO
# Line 174  C     vertical shear term (dU/dz)^2+(dV/ Line 177  C     vertical shear term (dU/dz)^2+(dV/
177  C     compute Prandtl number (always greater than 0)  C     compute Prandtl number (always greater than 0)
178           prTemp = 1. _d 0           prTemp = 1. _d 0
179           IF ( RiNumber .GE. 0.2 ) prTemp = 5.0 * RiNumber           IF ( RiNumber .GE. 0.2 ) prTemp = 5.0 * RiNumber
180           TKEPrandtlNumber(I,J,K) = MIN(10.,prTemp)           TKEPrandtlNumber(I,J,K) = MIN(10.0 _d 0,prTemp)
181  C     mixing length  C     mixing length
182           GGL90mixingLength(I,J,K) =           GGL90mixingLength(I,J,K) =
183       &        SQRTTKE/SQRT( MAX(Nsquare,GGL90eps) )       &        SQRTTKE/SQRT( MAX(Nsquare,GGL90eps) )
# Line 184  C     impose upper bound for mixing leng Line 187  C     impose upper bound for mixing leng
187  C     impose minimum mixing length (to avoid division by zero)  C     impose minimum mixing length (to avoid division by zero)
188           GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),           GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),
189       &        GGL90mixingLengthMin)       &        GGL90mixingLengthMin)
190             rMixingLength(I,J,K) = 1. _d 0 /GGL90mixingLength(I,J,K)
191  C     viscosity of last timestep  C     viscosity of last timestep
192           KappaM = GGL90ck*GGL90mixingLength(I,J,K)*SQRTTKE           KappaM = GGL90ck*GGL90mixingLength(I,J,K)*SQRTTKE
193           KappaE(I,J,K) = KappaM*GGL90alpha           KappaE(I,J,K) = KappaM*GGL90alpha
194  C     dissipation term  C     dissipation term
195           TKEdissipation = ab05*GGL90ceps           TKEdissipation = ab05*GGL90ceps
196       &        *SQRTTKE/GGL90mixingLength(I,J,K)       &        *SQRTTKE*rMixingLength(I,J,K)
197       &        *GGL90TKE(I,J,K,bi,bj)             &        *GGL90TKE(I,J,K,bi,bj)      
198  C     sum up contributions to form the right hand side  C     sum up contributions to form the right hand side
199           gTKE(I,J,K) = GGL90TKE(I,J,K,bi,bj)           gTKE(I,J,K) = GGL90TKE(I,J,K,bi,bj)
# Line 305  C--   Center diagonal Line 309  C--   Center diagonal
309          DO i=iMin,iMax          DO i=iMin,iMax
310            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)
311       &        + ab15*deltaTggl90*GGL90ceps*SQRT(GGL90TKE(I,J,K,bi,bj))       &        + ab15*deltaTggl90*GGL90ceps*SQRT(GGL90TKE(I,J,K,bi,bj))
312       &        /GGL90mixingLength(I,J,K)       &        *rMixingLength(I,J,K)
313           ENDDO           ENDDO
314         ENDDO         ENDDO
315        ENDDO        ENDDO
# Line 324  C     estimate friction velocity uStar f Line 328  C     estimate friction velocity uStar f
328       &              + surfaceForcingV(I,  J+1,bi,bj) ) )**2       &              + surfaceForcingV(I,  J+1,bi,bj) ) )**2
329       &                     )       &                     )
330  C     Dirichlet surface boundary condition for TKE  C     Dirichlet surface boundary condition for TKE
331          gTKE(I,J,kSurf) = MAX(GGL90TKEmin,GGL90m2*uStarSquare)          gTKE(I,J,kSurf) = MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)
332       &                     *maskC(I,J,kSurf,bi,bj)       &                     *maskC(I,J,kSurf,bi,bj)
333  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom
334          kBottom   = MIN(MAX(kLowC(I,J,bi,bj),1),Nr)          kBottom   = MIN(MAX(kLowC(I,J,bi,bj),1),Nr)
# Line 365  C     Eq. (11), (18) and (21) Line 369  C     Eq. (11), (18) and (21)
369           KappaH = KappaM/TKEPrandtlNumber(I,J,K)           KappaH = KappaM/TKEPrandtlNumber(I,J,K)
370  C     Set a minium (= background) value  C     Set a minium (= background) value
371           KappaM = MAX(KappaM,viscAr)           KappaM = MAX(KappaM,viscAr)
372           KappaH = MAX(KappaH,diffKrT)           KappaH = MAX(KappaH,diffKrNrT(k))
373  C     Set a maximum and mask land point  C     Set a maximum and mask land point
374           GGL90viscAr(I,J,K,bi,bj) = MIN(KappaM,GGL90viscMax)           GGL90viscAr(I,J,K,bi,bj) = MIN(KappaM,GGL90viscMax)
375       &        * maskC(I,J,K,bi,bj)       &        * maskC(I,J,K,bi,bj)
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.7
  ViewVC Help
Powered by ViewVC 1.1.22