/[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.18 by gforget, Wed Aug 11 03:32:29 2010 UTC revision 1.19 by jmc, Thu Aug 19 23:52:37 2010 UTC
# Line 7  CBOP Line 7  CBOP
7  C !ROUTINE: GGL90_CALC  C !ROUTINE: GGL90_CALC
8    
9  C !INTERFACE: ======================================================  C !INTERFACE: ======================================================
10        subroutine GGL90_CALC(        SUBROUTINE GGL90_CALC(
11       I     bi, bj, myTime, myThid )       I     bi, bj, myTime, myIter, myThid )
12    
13  C !DESCRIPTION: \bv  C !DESCRIPTION: \bv
14  C     *==========================================================*  C     *==========================================================*
# Line 43  C !INPUT PARAMETERS: =================== Line 43  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     myIter :: Current time-step number
47  C     myThid :: My Thread Id number  C     myThid :: My Thread Id number
48        INTEGER bi, bj        INTEGER bi, bj
49        _RL     myTime        _RL     myTime
50          INTEGER myIter
51        INTEGER myThid        INTEGER myThid
52  CEOP  CEOP
53    
# Line 53  CEOP Line 55  CEOP
55    
56  C !LOCAL VARIABLES: ====================================================  C !LOCAL VARIABLES: ====================================================
57  C Local constants  C Local constants
58  C     iMin, iMax, jMin, jMax, I, J - array computation indices  C     iMin,iMax,jMin,jMax :: index boundaries of computation domain
59  C     K, Kp1, km1, kSurf, kBottom  - vertical loop indices  C     i, j, k, kp1,km1 :: array computation indices
60  C     ab15, ab05       - weights for implicit timestepping  C     kSurf, kBottom   :: vertical indices of domain boundaries
61  C     uStarSquare      - square of friction velocity  C     explDissFac      :: explicit Dissipation Factor (in [0-1])
62  C     verticalShear    - (squared) vertical shear of horizontal velocity  C     implDissFac      :: implicit Dissipation Factor (in [0-1])
63  C     Nsquare          - squared buoyancy freqency  C     uStarSquare      :: square of friction velocity
64  C     RiNumber         - local Richardson number  C     verticalShear    :: (squared) vertical shear of horizontal velocity
65  C     KappaM           - (local) viscosity parameter (eq.10)  C     Nsquare          :: squared buoyancy freqency
66  C     KappaH           - (local) diffusivity parameter for temperature (eq.11)  C     RiNumber         :: local Richardson number
67  C     KappaE           - (local) diffusivity parameter for TKE (eq.15)  C     KappaM           :: (local) viscosity parameter (eq.10)
68  C     TKEdissipation   - dissipation of TKE  C     KappaH           :: (local) diffusivity parameter for temperature (eq.11)
69  C     GGL90mixingLength- mixing length of scheme following Banke+Delecuse  C     KappaE           :: (local) diffusivity parameter for TKE (eq.15)
70  C         rMixingLength- inverse of mixing length  C     TKEdissipation   :: dissipation of TKE
71  C     totalDepth       - thickness of water column (inverse of recip_Rcol)  C     GGL90mixingLength:: mixing length of scheme following Banke+Delecuse
72  C     TKEPrandtlNumber - here, an empirical function of the Richardson number  C         rMixingLength:: inverse of mixing length
73  C     rhoK, rhoKm1     - density at layer K and Km1 (relative to K)  C     totalDepth       :: thickness of water column (inverse of recip_Rcol)
74  C     gTKE             - right hand side of implicit equation  C     TKEPrandtlNumber :: here, an empirical function of the Richardson number
75    C     rhoK, rhoKm1     :: density at layer k and km1 (relative to k)
76        INTEGER iMin ,iMax ,jMin ,jMax        INTEGER iMin ,iMax ,jMin ,jMax
77        INTEGER I, J, K, Kp1, Km1, kSurf, kBottom        INTEGER i, j, k, kp1, km1, kSurf, kBottom
78        _RL     ab15, ab05        _RL     explDissFac, implDissFac
79        _RL     uStarSquare        _RL     uStarSquare
80        _RL     verticalShear        _RL     verticalShear
81        _RL     KappaM, KappaH        _RL     KappaM, KappaH
# Line 93  c     _RL     SQRTTKE Line 96  c     _RL     SQRTTKE
96        _RL     rhoK             (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     rhoK             (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
97        _RL     rhoKm1           (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     rhoKm1           (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98        _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL     gTKE             (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
99        _RL     GGL90visctmp     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     GGL90visctmp     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
100  C-    tri-diagonal matrix  C-    tri-diagonal matrix
101        _RL     a(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     a(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
# Line 101  C-    tri-diagonal matrix Line 103  C-    tri-diagonal matrix
103        _RL     c(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     c(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
104        INTEGER errCode        INTEGER errCode
105  #ifdef ALLOW_GGL90_HORIZDIFF  #ifdef ALLOW_GGL90_HORIZDIFF
106  C-    xA, yA   - area of lateral faces  C     xA, yA   :: area of lateral faces
107  C-    dfx, dfy - diffusive flux across lateral faces  C     dfx, dfy :: diffusive flux across lateral faces
108    C     gTKE     :: right hand side of diffusion equation
109        _RL     xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
110        _RL     yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
111        _RL     dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
112        _RL     dfy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     dfy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
113          _RL    gTKE(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
114  #endif /* ALLOW_GGL90_HORIZDIFF */  #endif /* ALLOW_GGL90_HORIZDIFF */
115  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
116        _RL p4, p8, p16        _RL p4, p8, p16
# Line 123  C     set separate time step (should be Line 127  C     set separate time step (should be
127        deltaTggl90 = dTtracerLev(1)        deltaTggl90 = dTtracerLev(1)
128    
129        kSurf = 1        kSurf = 1
130  C     implicit timestepping weights for dissipation  C     explicit/implicit timestepping weights for dissipation
131        ab15 =  1.5 _d 0        explDissFac = 0. _d 0
132        ab05 = -0.5 _d 0        implDissFac = 1. _d 0 - explDissFac
       ab15 =  1. _d 0  
       ab05 =  0. _d 0  
133    
134  C     Initialize local fields  C     Initialize local fields
135        DO K = 1, Nr        DO k = 1, Nr
136         DO J=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
137          DO I=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
138           gTKE(I,J,K)              = 0. _d 0           KappaE(i,j,k)            = 0. _d 0
139           KappaE(I,J,K)            = 0. _d 0           TKEPrandtlNumber(i,j,k)  = 1. _d 0
140           TKEPrandtlNumber(I,J,K)  = 1. _d 0           GGL90mixingLength(i,j,k) = GGL90mixingLengthMin
141           GGL90mixingLength(I,J,K) = GGL90mixingLengthMin           GGL90visctmp(i,j,k)      = 0. _d 0
          GGL90visctmp(I,J,K)      = 0. _d 0  
142          ENDDO          ENDDO
143         ENDDO         ENDDO
144        ENDDO        ENDDO
145        DO J=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
146         DO I=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
147          rhoK(I,J)          = 0. _d 0          rhoK(i,j)          = 0. _d 0
148          rhoKm1(I,J)        = 0. _d 0          rhoKm1(i,j)        = 0. _d 0
149          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)
150          rMixingLength(i,j,1) = 0. _d 0          rMixingLength(i,j,1) = 0. _d 0
151          mxLength_Dn(I,J,1) = GGL90mixingLengthMin          mxLength_Dn(i,j,1) = GGL90mixingLengthMin
152          SQRTTKE(i,j,1) = SQRT( GGL90TKE(i,j,1,bi,bj) )          SQRTTKE(i,j,1) = SQRT( GGL90TKE(i,j,1,bi,bj) )
153         ENDDO         ENDDO
154        ENDDO        ENDDO
155    
156  C     start k-loop  C     start k-loop
157        DO K = 2, Nr        DO k = 2, Nr
158         Km1 = K-1         km1 = k-1
159  c      Kp1 = MIN(Nr,K+1)  c      kp1 = MIN(Nr,k+1)
160         CALL FIND_RHO_2D(         CALL FIND_RHO_2D(
161       I      iMin, iMax, jMin, jMax, K,       I      iMin, iMax, jMin, jMax, k,
162       I      theta(1-OLx,1-OLy,Km1,bi,bj), salt(1-OLx,1-OLy,Km1,bi,bj),       I      theta(1-OLx,1-OLy,km1,bi,bj), salt(1-OLx,1-OLy,km1,bi,bj),
163       O      rhoKm1,       O      rhoKm1,
164       I      Km1, bi, bj, myThid )       I      km1, bi, bj, myThid )
165    
166         CALL FIND_RHO_2D(         CALL FIND_RHO_2D(
167       I      iMin, iMax, jMin, jMax, K,       I      iMin, iMax, jMin, jMax, k,
168       I      theta(1-OLx,1-OLy,K,bi,bj), salt(1-OLx,1-OLy,K,bi,bj),       I      theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
169       O      rhoK,       O      rhoK,
170       I      K, bi, bj, myThid )       I      k, bi, bj, myThid )
171         DO J=jMin,jMax         DO j=jMin,jMax
172          DO I=iMin,iMax          DO i=iMin,iMax
173           SQRTTKE(i,j,k)=SQRT( GGL90TKE(I,J,K,bi,bj) )           SQRTTKE(i,j,k)=SQRT( GGL90TKE(i,j,k,bi,bj) )
174    
175  C     buoyancy frequency  C     buoyancy frequency
176           Nsquare(i,j,k) = - gravity*recip_rhoConst*recip_drC(K)           Nsquare(i,j,k) = - gravity*recip_rhoConst*recip_drC(k)
177       &        * ( rhoKm1(I,J) - rhoK(I,J) )*maskC(I,J,K,bi,bj)       &        * ( rhoKm1(i,j) - rhoK(i,j) )*maskC(i,j,k,bi,bj)
178  cC     vertical shear term (dU/dz)^2+(dV/dz)^2  cC     vertical shear term (dU/dz)^2+(dV/dz)^2
179  c         tempU= .5 _d 0*( uVel(I,J,Km1,bi,bj)+uVel(I+1,J,Km1,bi,bj)  c         tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)
180  c     &                 -( uVel(I,J,K  ,bi,bj)+uVel(I+1,J,K  ,bi,bj)) )  c     &                 -( uVel(i,j,k  ,bi,bj)+uVel(i+1,j,k  ,bi,bj)) )
181  c     &        *recip_drC(K)  c     &        *recip_drC(k)
182  c         tempV= .5 _d 0*( vVel(I,J,Km1,bi,bj)+vVel(I,J+1,Km1,bi,bj)  c         tempV= .5 _d 0*( vVel(i,j,km1,bi,bj)+vVel(i,j+1,km1,bi,bj)
183  c     &                 -( vVel(I,J,K  ,bi,bj)+vVel(I,J+1,K  ,bi,bj)) )  c     &                 -( vVel(i,j,k  ,bi,bj)+vVel(i,j+1,k  ,bi,bj)) )
184  c     &        *recip_drC(K)  c     &        *recip_drC(k)
185  c         verticalShear = tempU*tempU + tempV*tempV  c         verticalShear = tempU*tempU + tempV*tempV
186  c         RiNumber   = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)  c         RiNumber   = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)
187  cC     compute Prandtl number (always greater than 0)  cC     compute Prandtl number (always greater than 0)
188  c         prTemp = 1. _d 0  c         prTemp = 1. _d 0
189  c         IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber  c         IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
190  c         TKEPrandtlNumber(I,J,K) = MIN(10. _d 0,prTemp)  c         TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
191  C     mixing length  C     mixing length
192           GGL90mixingLength(I,J,K) = SQRTTWO *           GGL90mixingLength(i,j,k) = SQRTTWO *
193       &        SQRTTKE(i,j,k)/SQRT( MAX(Nsquare(i,j,k),GGL90eps) )       &        SQRTTKE(i,j,k)/SQRT( MAX(Nsquare(i,j,k),GGL90eps) )
194          ENDDO          ENDDO
195         ENDDO         ENDDO
# Line 198  C- Impose upper and lower bound for mixi Line 199  C- Impose upper and lower bound for mixi
199        IF ( mxlMaxFlag .EQ. 0 ) THEN        IF ( mxlMaxFlag .EQ. 0 ) THEN
200  C-  C-
201         DO k=2,Nr         DO k=2,Nr
202          DO J=jMin,jMax          DO j=jMin,jMax
203           DO I=iMin,iMax           DO i=iMin,iMax
204            MaxLength=totalDepth(I,J)            MaxLength=totalDepth(i,j)
205            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
206       &                                   MaxLength)       &                                   MaxLength)
207           ENDDO           ENDDO
208          ENDDO          ENDDO
209         ENDDO         ENDDO
210    
211         DO k=2,Nr         DO k=2,Nr
212          DO J=jMin,jMax          DO j=jMin,jMax
213           DO I=iMin,iMax           DO i=iMin,iMax
214            GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
215       &                                   GGL90mixingLengthMin)       &                                   GGL90mixingLengthMin)
216            rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K)            rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
217           ENDDO           ENDDO
218          ENDDO          ENDDO
219         ENDDO         ENDDO
# Line 220  C- Line 221  C-
221        ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN
222  C-  C-
223         DO k=2,Nr         DO k=2,Nr
224          DO J=jMin,jMax          DO j=jMin,jMax
225           DO I=iMin,iMax           DO i=iMin,iMax
226            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))
227  c         MaxLength=MAX(MaxLength,20. _d 0)  c         MaxLength=MAX(MaxLength,20. _d 0)
228            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
229       &                                   MaxLength)       &                                   MaxLength)
230           ENDDO           ENDDO
231          ENDDO          ENDDO
232         ENDDO         ENDDO
233    
234         DO k=2,Nr         DO k=2,Nr
235          DO J=jMin,jMax          DO j=jMin,jMax
236           DO I=iMin,iMax           DO i=iMin,iMax
237            GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
238       &                                   GGL90mixingLengthMin)       &                                   GGL90mixingLengthMin)
239            rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K)            rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
240           ENDDO           ENDDO
241          ENDDO          ENDDO
242         ENDDO         ENDDO
# Line 243  c         MaxLength=MAX(MaxLength,20. _d Line 244  c         MaxLength=MAX(MaxLength,20. _d
244        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN
245  C-  C-
246  cgf ensure mixing between first and second level  cgf ensure mixing between first and second level
247  c      DO J=jMin,jMax  c      DO j=jMin,jMax
248  c        DO I=iMin,iMax  c        DO i=iMin,iMax
249  c         GGL90mixingLength(I,J,2)=drF(1)  c         GGL90mixingLength(i,j,2)=drF(1)
250  c        ENDDO  c        ENDDO
251  c      ENDDO  c      ENDDO
252  cgf  cgf
253         DO k=2,Nr         DO k=2,Nr
254          DO J=jMin,jMax          DO j=jMin,jMax
255           DO I=iMin,iMax           DO i=iMin,iMax
256            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
257       &        GGL90mixingLength(I,J,K-1)+drF(k-1))       &        GGL90mixingLength(i,j,k-1)+drF(k-1))
258           ENDDO           ENDDO
259          ENDDO          ENDDO
260         ENDDO         ENDDO
261         DO J=jMin,jMax         DO j=jMin,jMax
262          DO I=iMin,iMax          DO i=iMin,iMax
263            GGL90mixingLength(I,J,Nr) = MIN(GGL90mixingLength(I,J,Nr),            GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
264       &       GGL90mixingLengthMin+drF(Nr))       &       GGL90mixingLengthMin+drF(Nr))
265          ENDDO          ENDDO
266         ENDDO         ENDDO
267         DO k=Nr-1,2,-1         DO k=Nr-1,2,-1
268          DO J=jMin,jMax          DO j=jMin,jMax
269           DO I=iMin,iMax           DO i=iMin,iMax
270            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
271       &        GGL90mixingLength(I,J,K+1)+drF(k))       &        GGL90mixingLength(i,j,k+1)+drF(k))
272           ENDDO           ENDDO
273          ENDDO          ENDDO
274         ENDDO         ENDDO
275    
276         DO k=2,Nr         DO k=2,Nr
277          DO J=jMin,jMax          DO j=jMin,jMax
278           DO I=iMin,iMax           DO i=iMin,iMax
279            GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
280       &                                   GGL90mixingLengthMin)       &                                   GGL90mixingLengthMin)
281            rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K)            rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
282           ENDDO           ENDDO
283          ENDDO          ENDDO
284         ENDDO         ENDDO
# Line 285  cgf Line 286  cgf
286        ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN
287  C-  C-
288         DO k=2,Nr         DO k=2,Nr
289          DO J=jMin,jMax          DO j=jMin,jMax
290           DO I=iMin,iMax           DO i=iMin,iMax
291            mxLength_Dn(I,J,K) = MIN(GGL90mixingLength(I,J,K),            mxLength_Dn(i,j,k) = MIN(GGL90mixingLength(i,j,k),
292       &        mxLength_Dn(I,J,K-1)+drF(k-1))       &        mxLength_Dn(i,j,k-1)+drF(k-1))
293           ENDDO           ENDDO
294          ENDDO          ENDDO
295         ENDDO         ENDDO
296         DO J=jMin,jMax         DO j=jMin,jMax
297          DO I=iMin,iMax          DO i=iMin,iMax
298            GGL90mixingLength(I,J,Nr) = MIN(GGL90mixingLength(I,J,Nr),            GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
299       &       GGL90mixingLengthMin+drF(Nr))       &       GGL90mixingLengthMin+drF(Nr))
300          ENDDO          ENDDO
301         ENDDO         ENDDO
302         DO k=Nr-1,2,-1         DO k=Nr-1,2,-1
303          DO J=jMin,jMax          DO j=jMin,jMax
304           DO I=iMin,iMax           DO i=iMin,iMax
305            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
306       &        GGL90mixingLength(I,J,K+1)+drF(k))       &        GGL90mixingLength(i,j,k+1)+drF(k))
307           ENDDO           ENDDO
308          ENDDO          ENDDO
309         ENDDO         ENDDO
310    
311         DO k=2,Nr         DO k=2,Nr
312          DO J=jMin,jMax          DO j=jMin,jMax
313           DO I=iMin,iMax           DO i=iMin,iMax
314            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
315       &                                  mxLength_Dn(I,J,K))       &                                  mxLength_Dn(i,j,k))
316            tmpmlx = SQRT( GGL90mixingLength(I,J,K)*mxLength_Dn(I,J,K) )            tmpmlx = SQRT( GGL90mixingLength(i,j,k)*mxLength_Dn(i,j,k) )
317            tmpmlx = MAX( tmpmlx, GGL90mixingLengthMin)            tmpmlx = MAX( tmpmlx, GGL90mixingLengthMin)
318            rMixingLength(I,J,K) = 1. _d 0 / tmpmlx            rMixingLength(i,j,k) = 1. _d 0 / tmpmlx
319           ENDDO           ENDDO
320          ENDDO          ENDDO
321         ENDDO         ENDDO
# Line 325  C- Line 326  C-
326    
327  C- Impose minimum mixing length (to avoid division by zero)  C- Impose minimum mixing length (to avoid division by zero)
328  c      DO k=2,Nr  c      DO k=2,Nr
329  c      DO J=jMin,jMax  c      DO j=jMin,jMax
330  c       DO I=iMin,iMax  c       DO i=iMin,iMax
331  c        GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),  c        GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
332  c    &        GGL90mixingLengthMin)  c    &        GGL90mixingLengthMin)
333  c        rMixingLength(I,J,K) = 1. _d 0 /GGL90mixingLength(I,J,K)  c        rMixingLength(i,j,k) = 1. _d 0 /GGL90mixingLength(i,j,k)
334  c       ENDDO  c       ENDDO
335  c      ENDDO  c      ENDDO
336  c     ENDDO  c     ENDDO
337    
   
338        DO k=2,Nr        DO k=2,Nr
339         Km1 = K-1         km1 = k-1
        DO J=jMin,jMax  
         DO I=iMin,iMax  
 C     vertical shear term (dU/dz)^2+(dV/dz)^2  
          tempU= .5 _d 0*( uVel(I,J,Km1,bi,bj)+uVel(I+1,J,Km1,bi,bj)  
      &                 -( uVel(I,J,K  ,bi,bj)+uVel(I+1,J,K  ,bi,bj)) )  
      &        *recip_drC(K)  
          tempV= .5 _d 0*( vVel(I,J,Km1,bi,bj)+vVel(I,J+1,Km1,bi,bj)  
      &                 -( vVel(I,J,K  ,bi,bj)+vVel(I,J+1,K  ,bi,bj)) )  
      &        *recip_drC(K)  
          verticalShear = tempU*tempU + tempV*tempV  
          RiNumber = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)  
 C     compute Prandtl number (always greater than 0)  
          prTemp = 1. _d 0  
          IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber  
          TKEPrandtlNumber(I,J,K) = MIN(10. _d 0,prTemp)  
 c         TKEPrandtlNumber(I,J,K) = 1. _d 0  
   
 C     viscosity and diffusivity  
          KappaM = GGL90ck*GGL90mixingLength(I,J,K)*SQRTTKE(i,j,k)  
          GGL90visctmp(I,J,K) = MAX(KappaM,diffKrNrT(k))  
      &                            * maskC(I,J,K,bi,bj)  
 c        note: storing GGL90visctmp like this, and using it later to compute  
 c              GGL9rdiffKr etc. is robust in case of smoothing (e.g. see OPA)  
          KappaM = MAX(KappaM,viscArNr(k)) * maskC(I,J,K,bi,bj)  
          KappaH = KappaM/TKEPrandtlNumber(I,J,K)  
          KappaE(I,J,K) = GGL90alpha * KappaM * maskC(I,J,K,bi,bj)  
   
 C     dissipation term  
          TKEdissipation = ab05*GGL90ceps  
      &        *SQRTTKE(i,j,k)*rMixingLength(I,J,K)  
      &        *GGL90TKE(I,J,K,bi,bj)  
 C     sum up contributions to form the right hand side  
          gTKE(I,J,K) = GGL90TKE(I,J,K,bi,bj)  
      &        + deltaTggl90*(  
      &        + KappaM*verticalShear  
      &        - KappaH*Nsquare(i,j,k)  
      &        - TKEdissipation  
      &        )  
         ENDDO  
        ENDDO  
       ENDDO  
340    
341    #ifdef ALLOW_GGL90_HORIZDIFF
342           IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
343  C     horizontal diffusion of TKE (requires an exchange in  C     horizontal diffusion of TKE (requires an exchange in
344  C      do_fields_blocking_exchanges)  C      do_fields_blocking_exchanges)
 #ifdef ALLOW_GGL90_HORIZDIFF  
       IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN  
        DO K=2,Nr  
345  C     common factors  C     common factors
346          DO j=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
347           DO i=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
# Line 421  C     ... across y-faces Line 379  C     ... across y-faces
379  C     Compute divergence of fluxes  C     Compute divergence of fluxes
380          DO j=1-Oly,sNy+Oly-1          DO j=1-Oly,sNy+Oly-1
381           DO i=1-Olx,sNx+Olx-1           DO i=1-Olx,sNx+Olx-1
382            gTKE(i,j,k)=gTKE(i,j,k)            gTKE(i,j) =
383       &   -_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)
384       &         *( (dfx(i+1,j)-dfx(i,j))       &         *( (dfx(i+1,j)-dfx(i,j))
385       &           +(dfy(i,j+1)-dfy(i,j))       &           +(dfy(i,j+1)-dfy(i,j))
386       &           )*deltaTggl90       &          )
387           ENDDO           ENDDO
388          ENDDO          ENDDO
389  C       end of k-loop  C      end if GGL90diffTKEh .eq. 0.
390           ENDIF
391    #endif /* ALLOW_GGL90_HORIZDIFF */
392    
393           DO j=jMin,jMax
394            DO i=iMin,iMax
395    C     vertical shear term (dU/dz)^2+(dV/dz)^2
396             tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)
397         &                 -( uVel(i,j,k  ,bi,bj)+uVel(i+1,j,k  ,bi,bj)) )
398         &        *recip_drC(k)
399             tempV= .5 _d 0*( vVel(i,j,km1,bi,bj)+vVel(i,j+1,km1,bi,bj)
400         &                 -( vVel(i,j,k  ,bi,bj)+vVel(i,j+1,k  ,bi,bj)) )
401         &        *recip_drC(k)
402             verticalShear = tempU*tempU + tempV*tempV
403             RiNumber = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)
404    C     compute Prandtl number (always greater than 0)
405             prTemp = 1. _d 0
406             IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
407             TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
408    c         TKEPrandtlNumber(i,j,k) = 1. _d 0
409    
410    C     viscosity and diffusivity
411             KappaM = GGL90ck*GGL90mixingLength(i,j,k)*SQRTTKE(i,j,k)
412             GGL90visctmp(i,j,k) = MAX(KappaM,diffKrNrT(k))
413         &                            * maskC(i,j,k,bi,bj)
414    c        note: storing GGL90visctmp like this, and using it later to compute
415    c              GGL9rdiffKr etc. is robust in case of smoothing (e.g. see OPA)
416             KappaM = MAX(KappaM,viscArNr(k)) * maskC(i,j,k,bi,bj)
417             KappaH = KappaM/TKEPrandtlNumber(i,j,k)
418             KappaE(i,j,k) = GGL90alpha * KappaM * maskC(i,j,k,bi,bj)
419    
420    C     dissipation term
421             TKEdissipation = explDissFac*GGL90ceps
422         &        *SQRTTKE(i,j,k)*rMixingLength(i,j,k)
423         &        *GGL90TKE(i,j,k,bi,bj)
424    C     partial update with sum of explicit contributions
425             GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
426         &        + deltaTggl90*(
427         &        + KappaM*verticalShear
428         &        - KappaH*Nsquare(i,j,k)
429         &        - TKEdissipation
430         &        )
431            ENDDO
432         ENDDO         ENDDO
433  C     end if GGL90diffTKEh .eq. 0.  
434        ENDIF  #ifdef ALLOW_GGL90_HORIZDIFF
435           IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
436    C--    Add horiz. diffusion tendency
437            DO j=jMin,jMax
438             DO i=iMin,iMax
439              GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
440         &                          + gTKE(i,j)*deltaTggl90
441             ENDDO
442            ENDDO
443           ENDIF
444  #endif /* ALLOW_GGL90_HORIZDIFF */  #endif /* ALLOW_GGL90_HORIZDIFF */
445    
446    C--   end of k loop
447          ENDDO
448    
449  C     ============================================  C     ============================================
450  C     Implicit time step to update TKE for k=1,Nr;  C     Implicit time step to update TKE for k=1,Nr;
451  C     TKE(Nr+1)=0 by default  C     TKE(Nr+1)=0 by default
# Line 485  C--   Center diagonal Line 497  C--   Center diagonal
497         DO j=jMin,jMax         DO j=jMin,jMax
498          DO i=iMin,iMax          DO i=iMin,iMax
499            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)
500       &        + ab15*deltaTggl90*GGL90ceps*SQRTTKE(I,J,K)       &        + implDissFac*deltaTggl90*GGL90ceps*SQRTTKE(i,j,k)
501       &        * rMixingLength(I,J,K)       &        * rMixingLength(i,j,k)
502       &        * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)       &        * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)
503           ENDDO           ENDDO
504         ENDDO         ENDDO
# Line 495  C     end set up matrix Line 507  C     end set up matrix
507    
508  C     Apply boundary condition  C     Apply boundary condition
509        kp1 = MIN(Nr,kSurf+1)        kp1 = MIN(Nr,kSurf+1)
510        DO J=jMin,jMax        DO j=jMin,jMax
511         DO I=iMin,iMax         DO i=iMin,iMax
512  C     estimate friction velocity uStar from surface forcing  C     estimate friction velocity uStar from surface forcing
513          uStarSquare = SQRT(          uStarSquare = SQRT(
514       &    ( .5 _d 0*( surfaceForcingU(I,  J,  bi,bj)       &    ( .5 _d 0*( surfaceForcingU(i,  j,  bi,bj)
515       &              + surfaceForcingU(I+1,J,  bi,bj) ) )**2       &              + surfaceForcingU(i+1,j,  bi,bj) ) )**2
516       &  + ( .5 _d 0*( surfaceForcingV(I,  J,  bi,bj)       &  + ( .5 _d 0*( surfaceForcingV(i,  j,  bi,bj)
517       &              + surfaceForcingV(I,  J+1,bi,bj) ) )**2       &              + surfaceForcingV(i,  j+1,bi,bj) ) )**2
518       &                     )       &                     )
519  C     Dirichlet surface boundary condition for TKE  C     Dirichlet surface boundary condition for TKE
520          gTKE(I,J,kSurf) = MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)          GGL90TKE(i,j,kSurf,bi,bj) = maskC(i,j,kSurf,bi,bj)
521       &                     *maskC(I,J,kSurf,bi,bj)       &           *MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)
522          gTKE(i,j,kp1) = gTKE(i,j,kp1)          GGL90TKE(i,j,kp1,bi,bj) = GGL90TKE(i,j,kp1,bi,bj)
523       &                - a(i,j,kp1)*gTKE(i,j,kSurf)       &               - a(i,j,kp1)*GGL90TKE(i,j,kSurf,bi,bj)
524          a(i,j,kp1) = 0. _d 0          a(i,j,kp1) = 0. _d 0
525  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom
526          kBottom   = MAX(kLowC(I,J,bi,bj),1)          kBottom   = MAX(kLowC(i,j,bi,bj),1)
527          gTKE(I,J,kBottom) = gTKE(I,J,kBottom)          GGL90TKE(i,j,kBottom,bi,bj) = GGL90TKE(i,j,kBottom,bi,bj)
528       &                    - GGL90TKEbottom*c(I,J,kBottom)       &                              - GGL90TKEbottom*c(i,j,kBottom)
529          c(I,J,kBottom) = 0. _d 0          c(i,j,kBottom) = 0. _d 0
530         ENDDO         ENDDO
531        ENDDO        ENDDO
532    
533  C     solve tri-diagonal system, and store solution on gTKE (previously rhs)  C     solve tri-diagonal system
534        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
535       I                        a, b, c,       I                        a, b, c,
536       U                        gTKE,       U                        GGL90TKE,
537       O                        errCode,       O                        errCode,
538       I                        1, 1, myThid )       I                        bi, bj, myThid )
539    
540  C     now update TKE        DO k=1,Nr
541        DO K=1,Nr         DO j=jMin,jMax
542         DO J=jMin,jMax          DO i=iMin,iMax
         DO I=iMin,iMax  
543  C     impose minimum TKE to avoid numerical undershoots below zero  C     impose minimum TKE to avoid numerical undershoots below zero
544           GGL90TKE(I,J,K,bi,bj) = MAX( gTKE(I,J,K), GGL90TKEmin )           GGL90TKE(i,j,k,bi,bj) = maskC(i,j,k,bi,bj)
545       &        * maskC(I,J,K,bi,bj)       &                  *MAX( GGL90TKE(i,j,k,bi,bj), GGL90TKEmin )
546          ENDDO          ENDDO
547         ENDDO         ENDDO
548        ENDDO        ENDDO
# Line 539  C     impose minimum TKE to avoid numeri Line 550  C     impose minimum TKE to avoid numeri
550  C     end of time step  C     end of time step
551  C     ===============================  C     ===============================
552    
553        DO K=2,Nr        DO k=2,Nr
554         DO J=1,sNy         DO j=1,sNy
555          DO I=1,sNx          DO i=1,sNx
556  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
557           tmpVisc=           tmpVisc=
558       &  (       &  (
# Line 566  C     =============================== Line 577  C     ===============================
577       &       +       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))
578       &  )*maskC(i,j,k,bi,bj)*mskCor(i,j,bi,bj)       &  )*maskC(i,j,k,bi,bj)*mskCor(i,j,bi,bj)
579  #else  #else
580           tmpVisc = GGL90visctmp(I,J,K)           tmpVisc = GGL90visctmp(i,j,k)
581  #endif  #endif
582           tmpVisc = MIN(tmpVisc/TKEPrandtlNumber(i,j,k),GGL90diffMax)           tmpVisc = MIN(tmpVisc/TKEPrandtlNumber(i,j,k),GGL90diffMax)
583           GGL90diffKr(I,J,K,bi,bj)= MAX( tmpVisc , diffKrNrT(k) )           GGL90diffKr(i,j,k,bi,bj)= MAX( tmpVisc , diffKrNrT(k) )
584          ENDDO          ENDDO
585         ENDDO         ENDDO
586        ENDDO        ENDDO
587    
588          DO k=2,Nr
589           DO j=1,sNy
590        DO K=2,Nr          DO i=1,sNx+1
        DO J=1,sNy  
         DO I=1,sNx+1  
591  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
592          tmpVisc =          tmpVisc =
593       & (       & (
594       &   p4 *(GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)       &   p4 *(GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
595       &       +GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj))       &       +GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj))
# Line 603  C     =============================== Line 612  C     ===============================
612       &                            +GGL90visctmp(i-1,j,k))       &                            +GGL90visctmp(i-1,j,k))
613       &                   )       &                   )
614  #endif  #endif
615          tmpVisc = MIN( tmpVisc , GGL90viscMax )          tmpVisc = MIN( tmpVisc , GGL90viscMax )
616          GGL90viscArU(i,j,k,bi,bj) = MAX( tmpVisc , viscArNr(k)  )          GGL90viscArU(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
617          ENDDO          ENDDO
618         ENDDO         ENDDO
619        ENDDO        ENDDO
620    
621          DO k=2,Nr
622        DO K=2,Nr         DO j=1,sNy+1
623         DO J=1,sNy+1          DO i=1,sNx
         DO I=1,sNx  
624  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
625          tmpVisc =          tmpVisc =
626       & (       & (
# Line 639  C     =============================== Line 647  C     ===============================
647    
648  #endif  #endif
649          tmpVisc = MIN( tmpVisc , GGL90viscMax )          tmpVisc = MIN( tmpVisc , GGL90viscMax )
650          GGL90viscArV(i,j,k,bi,bj) = MAX( tmpVisc , viscArNr(k)  )          GGL90viscArV(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
651          ENDDO          ENDDO
652         ENDDO         ENDDO
653        ENDDO        ENDDO
Legend:
Removed from v.1.18  
changed lines
  Added in v.1.19
  ViewVC Help
Powered by ViewVC 1.1.22