/[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.19 by jmc, Thu Aug 19 23:52:37 2010 UTC revision 1.29 by jmc, Sat Feb 21 17:13:20 2015 UTC
# Line 8  C !ROUTINE: GGL90_CALC Line 8  C !ROUTINE: GGL90_CALC
8    
9  C !INTERFACE: ======================================================  C !INTERFACE: ======================================================
10        SUBROUTINE GGL90_CALC(        SUBROUTINE GGL90_CALC(
11       I     bi, bj, myTime, myIter, myThid )       I                 bi, bj, sigmaR, myTime, myIter, myThid )
12    
13  C !DESCRIPTION: \bv  C !DESCRIPTION: \bv
14  C     *==========================================================*  C     *==========================================================*
# Line 41  C !USES: =============================== Line 41  C !USES: ===============================
41    
42  C !INPUT PARAMETERS: ===================================================  C !INPUT PARAMETERS: ===================================================
43  C Routine arguments  C Routine arguments
44  C     bi, bj :: array indices on which to apply calculations  C     bi, bj :: Current tile indices
45    C     sigmaR :: Vertical gradient of iso-neutral density
46  C     myTime :: Current time in simulation  C     myTime :: Current time in simulation
47  C     myIter :: Current time-step number  C     myIter :: Current time-step number
48  C     myThid :: My Thread Id number  C     myThid :: My Thread Id number
49        INTEGER bi, bj        INTEGER bi, bj
50          _RL     sigmaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
51        _RL     myTime        _RL     myTime
52        INTEGER myIter        INTEGER myIter
53        INTEGER myThid        INTEGER myThid
 CEOP  
54    
55  #ifdef ALLOW_GGL90  #ifdef ALLOW_GGL90
56    
# Line 72  C     GGL90mixingLength:: mixing length Line 73  C     GGL90mixingLength:: mixing length
73  C         rMixingLength:: inverse of mixing length  C         rMixingLength:: inverse of mixing length
74  C     totalDepth       :: thickness of water column (inverse of recip_Rcol)  C     totalDepth       :: thickness of water column (inverse of recip_Rcol)
75  C     TKEPrandtlNumber :: here, an empirical function of the Richardson number  C     TKEPrandtlNumber :: here, an empirical function of the Richardson number
 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     explDissFac, implDissFac        _RL     explDissFac, implDissFac
# Line 85  c     _RL     Nsquare Line 85  c     _RL     Nsquare
85  c     _RL     SQRTTKE  c     _RL     SQRTTKE
86        _RL     SQRTTKE(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     SQRTTKE(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
87        _RL     RiNumber        _RL     RiNumber
88    #ifdef ALLOW_GGL90_IDEMIX
89          _RL     IDEMIX_RiNumber
90    #endif
91        _RL     TKEdissipation        _RL     TKEdissipation
92        _RL     tempU, tempV, prTemp        _RL     tempU, tempV, prTemp
93        _RL     MaxLength, tmpmlx, tmpVisc        _RL     MaxLength, tmpmlx, tmpVisc
# Line 93  c     _RL     SQRTTKE Line 96  c     _RL     SQRTTKE
96        _RL     rMixingLength    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     rMixingLength    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
97        _RL     mxLength_Dn      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     mxLength_Dn      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
98        _RL     KappaE           (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     KappaE           (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
       _RL     rhoK             (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL     rhoKm1           (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
99        _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
100        _RL     GGL90visctmp     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     GGL90visctmp     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
101    #ifdef ALLOW_DIAGNOSTICS
102          _RL     surf_flx_tke     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
103    #endif /* ALLOW_DIAGNOSTICS */
104  C-    tri-diagonal matrix  C-    tri-diagonal matrix
105        _RL     a(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     a3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
106        _RL     b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     b3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
107        _RL     c(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     c3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
108        INTEGER errCode        INTEGER errCode
109  #ifdef ALLOW_GGL90_HORIZDIFF  #ifdef ALLOW_GGL90_HORIZDIFF
110    C     hFac     :: fractional thickness of W-cell
111  C     xA, yA   :: area of lateral faces  C     xA, yA   :: area of lateral faces
112  C     dfx, dfy :: diffusive flux across lateral faces  C     dfx, dfy :: diffusive flux across lateral faces
113  C     gTKE     :: right hand side of diffusion equation  C     gTKE     :: right hand side of diffusion equation
114          _RL     hFac
115        _RL     xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
116        _RL     yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
117        _RL     dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 114  C     gTKE     :: right hand side of dif Line 120  C     gTKE     :: right hand side of dif
120  #endif /* ALLOW_GGL90_HORIZDIFF */  #endif /* ALLOW_GGL90_HORIZDIFF */
121  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
122        _RL p4, p8, p16        _RL p4, p8, p16
123    CEOP
124        p4=0.25 _d 0        p4=0.25 _d 0
125        p8=0.125 _d 0        p8=0.125 _d 0
126        p16=0.0625 _d 0        p16=0.0625 _d 0
# Line 133  C     explicit/implicit timestepping wei Line 140  C     explicit/implicit timestepping wei
140    
141  C     Initialize local fields  C     Initialize local fields
142        DO k = 1, Nr        DO k = 1, Nr
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             rMixingLength(i,j,k)     = 0. _d 0
146             mxLength_Dn(i,j,k)       = 0. _d 0
147             GGL90visctmp(i,j,k)      = 0. _d 0
148           KappaE(i,j,k)            = 0. _d 0           KappaE(i,j,k)            = 0. _d 0
149           TKEPrandtlNumber(i,j,k)  = 1. _d 0           TKEPrandtlNumber(i,j,k)  = 1. _d 0
150           GGL90mixingLength(i,j,k) = GGL90mixingLengthMin           GGL90mixingLength(i,j,k) = GGL90mixingLengthMin
151           GGL90visctmp(i,j,k)      = 0. _d 0           GGL90visctmp(i,j,k)      = 0. _d 0
152    #ifndef SOLVE_DIAGONAL_LOWMEMORY
153             a3d(i,j,k) = 0. _d 0
154             b3d(i,j,k) = 1. _d 0
155             c3d(i,j,k) = 0. _d 0
156    #endif
157             Nsquare(i,j,k) = 0. _d 0
158             SQRTTKE(i,j,k) = 0. _d 0
159          ENDDO          ENDDO
160         ENDDO         ENDDO
161        ENDDO        ENDDO
162        DO j=1-Oly,sNy+Oly        DO j=1-OLy,sNy+OLy
163         DO i=1-Olx,sNx+Olx         DO i=1-OLx,sNx+OLx
         rhoK(i,j)          = 0. _d 0  
         rhoKm1(i,j)        = 0. _d 0  
164          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)
165          rMixingLength(i,j,1) = 0. _d 0          rMixingLength(i,j,1) = 0. _d 0
166          mxLength_Dn(i,j,1) = GGL90mixingLengthMin          mxLength_Dn(i,j,1) = GGL90mixingLengthMin
167          SQRTTKE(i,j,1) = SQRT( GGL90TKE(i,j,1,bi,bj) )          SQRTTKE(i,j,1) = SQRT( GGL90TKE(i,j,1,bi,bj) )
168    #ifdef ALLOW_GGL90_HORIZDIFF
169            xA(i,j)  = 0. _d 0
170            yA(i,j)  = 0. _d 0
171            dfx(i,j) = 0. _d 0
172            dfy(i,j) = 0. _d 0
173            gTKE(i,j) = 0. _d 0
174    #endif /* ALLOW_GGL90_HORIZDIFF */
175         ENDDO         ENDDO
176        ENDDO        ENDDO
177    
178    #ifdef ALLOW_GGL90_IDEMIX
179          IF ( useIDEMIX) CALL GGL90_IDEMIX(
180         &   bi, bj, sigmaR, myTime, myIter, myThid )
181    #endif /* ALLOW_GGL90_IDEMIX */
182    
183  C     start k-loop  C     start k-loop
184        DO k = 2, Nr        DO k = 2, Nr
185         km1 = k-1  c      km1 = k-1
186  c      kp1 = MIN(Nr,k+1)  c      kp1 = MIN(Nr,k+1)
        CALL FIND_RHO_2D(  
      I      iMin, iMax, jMin, jMax, k,  
      I      theta(1-OLx,1-OLy,km1,bi,bj), salt(1-OLx,1-OLy,km1,bi,bj),  
      O      rhoKm1,  
      I      km1, bi, bj, myThid )  
   
        CALL FIND_RHO_2D(  
      I      iMin, iMax, jMin, jMax, k,  
      I      theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),  
      O      rhoK,  
      I      k, bi, bj, myThid )  
187         DO j=jMin,jMax         DO j=jMin,jMax
188          DO i=iMin,iMax          DO i=iMin,iMax
189           SQRTTKE(i,j,k)=SQRT( GGL90TKE(i,j,k,bi,bj) )           SQRTTKE(i,j,k)=SQRT( GGL90TKE(i,j,k,bi,bj) )
190    
191  C     buoyancy frequency  C     buoyancy frequency
192           Nsquare(i,j,k) = - gravity*recip_rhoConst*recip_drC(k)           Nsquare(i,j,k) = gravity*gravitySign*recip_rhoConst
193       &        * ( rhoKm1(i,j) - rhoK(i,j) )*maskC(i,j,k,bi,bj)       &                  * sigmaR(i,j,k)
194  cC     vertical shear term (dU/dz)^2+(dV/dz)^2  cC     vertical shear term (dU/dz)^2+(dV/dz)^2
195  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)
196  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)) )
# Line 195  C     mixing length Line 211  C     mixing length
211         ENDDO         ENDDO
212        ENDDO        ENDDO
213    
214    C- ensure mixing between first and second level
215          IF (mxlSurfFlag) THEN
216           DO j=jMin,jMax
217            DO i=iMin,iMax
218             GGL90mixingLength(i,j,2)=drF(1)
219            ENDDO
220           ENDDO
221          ENDIF
222    
223  C- Impose upper and lower bound for mixing length  C- Impose upper and lower bound for mixing length
224        IF ( mxlMaxFlag .EQ. 0 ) THEN        IF ( mxlMaxFlag .EQ. 0 ) THEN
225  C-  
226         DO k=2,Nr         DO k=2,Nr
227          DO j=jMin,jMax          DO j=jMin,jMax
228           DO i=iMin,iMax           DO i=iMin,iMax
# Line 219  C- Line 244  C-
244         ENDDO         ENDDO
245    
246        ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN
247  C-  
248         DO k=2,Nr         DO k=2,Nr
249          DO j=jMin,jMax          DO j=jMin,jMax
250           DO i=iMin,iMax           DO i=iMin,iMax
# Line 242  c         MaxLength=MAX(MaxLength,20. _d Line 267  c         MaxLength=MAX(MaxLength,20. _d
267         ENDDO         ENDDO
268    
269        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN
270  C-  
 cgf ensure mixing between first and second level  
 c      DO j=jMin,jMax  
 c        DO i=iMin,iMax  
 c         GGL90mixingLength(i,j,2)=drF(1)  
 c        ENDDO  
 c      ENDDO  
 cgf  
271         DO k=2,Nr         DO k=2,Nr
272          DO j=jMin,jMax          DO j=jMin,jMax
273           DO i=iMin,iMax           DO i=iMin,iMax
# Line 284  cgf Line 302  cgf
302         ENDDO         ENDDO
303    
304        ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN
305  C-  
306         DO k=2,Nr         DO k=2,Nr
307          DO j=jMin,jMax          DO j=jMin,jMax
308           DO i=iMin,iMax           DO i=iMin,iMax
# Line 339  c     ENDDO Line 357  c     ENDDO
357         km1 = k-1         km1 = k-1
358    
359  #ifdef ALLOW_GGL90_HORIZDIFF  #ifdef ALLOW_GGL90_HORIZDIFF
360         IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN        IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
361  C     horizontal diffusion of TKE (requires an exchange in  C     horizontal diffusion of TKE (requires an exchange in
362  C      do_fields_blocking_exchanges)  C      do_fields_blocking_exchanges)
363  C     common factors  C     common factors
364          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
365           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
366            xA(i,j) = _dyG(i,j,bi,bj)            xA(i,j) = _dyG(i,j,bi,bj)*drC(k)*
367       &         *drF(k)*_hFacW(i,j,k,bi,bj)       &                 (min(.5 _d 0,_hFacW(i,j,k-1,bi,bj) ) +
368            yA(i,j) = _dxG(i,j,bi,bj)       &                  min(.5 _d 0,_hFacW(i,j,k  ,bi,bj) ) )
369       &         *drF(k)*_hFacS(i,j,k,bi,bj)            yA(i,j) = _dxG(i,j,bi,bj)*drC(k)*
370         &                 (min(.5 _d 0,_hFacS(i,j,k-1,bi,bj) ) +
371         &                  min(.5 _d 0,_hFacS(i,j,k  ,bi,bj) ) )
372           ENDDO           ENDDO
373          ENDDO          ENDDO
374  C     Compute diffusive fluxes  C     Compute diffusive fluxes
375  C     ... across x-faces  C     ... across x-faces
376          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
377           dfx(1-Olx,j)=0. _d 0           dfx(1-OLx,j)=0. _d 0
378           DO i=1-Olx+1,sNx+Olx           DO i=1-OLx+1,sNx+OLx
379            dfx(i,j) = -GGL90diffTKEh*xA(i,j)            dfx(i,j) = -GGL90diffTKEh*xA(i,j)
380       &      *_recip_dxC(i,j,bi,bj)       &      *_recip_dxC(i,j,bi,bj)
381       &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i-1,j,k,bi,bj))       &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i-1,j,k,bi,bj))
382    #ifdef ISOTROPIC_COS_SCALING
383       &      *CosFacU(j,bi,bj)       &      *CosFacU(j,bi,bj)
384    #endif /* ISOTROPIC_COS_SCALING */
385           ENDDO           ENDDO
386          ENDDO          ENDDO
387  C     ... across y-faces  C     ... across y-faces
388          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
389           dfy(i,1-Oly)=0. _d 0           dfy(i,1-OLy)=0. _d 0
390          ENDDO          ENDDO
391          DO j=1-Oly+1,sNy+Oly          DO j=1-OLy+1,sNy+OLy
392           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
393            dfy(i,j) = -GGL90diffTKEh*yA(i,j)            dfy(i,j) = -GGL90diffTKEh*yA(i,j)
394       &      *_recip_dyC(i,j,bi,bj)       &      *_recip_dyC(i,j,bi,bj)
395       &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i,j-1,k,bi,bj))       &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i,j-1,k,bi,bj))
# Line 377  C     ... across y-faces Line 399  C     ... across y-faces
399           ENDDO           ENDDO
400          ENDDO          ENDDO
401  C     Compute divergence of fluxes  C     Compute divergence of fluxes
402          DO j=1-Oly,sNy+Oly-1          DO j=1-OLy,sNy+OLy-1
403           DO i=1-Olx,sNx+Olx-1           DO i=1-OLx,sNx+OLx-1
404            gTKE(i,j) =  #ifdef ALLOW_GGL90_IDEMIX
405       &    -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj)            gTKE(i,j) = -recip_drC(k)*recip_rA(i,j,bi,bj)
406       &         *( (dfx(i+1,j)-dfx(i,j))       &         *recip_hFacI(i,j,k,bi,bj)
407       &           +(dfy(i,j+1)-dfy(i,j))  #else
408       &          )            hFac = MIN(.5 _d 0,_hFacC(i,j,k-1,bi,bj) ) +
409         &           MIN(.5 _d 0,_hFacC(i,j,k  ,bi,bj) )
410              gTKE(i,j) = 0.0
411              IF ( hFac .ne. 0.0 )
412         &      gTKE(i,j) = -recip_drC(k)*recip_rA(i,j,bi,bj)/hFac
413    #endif
414         &         *((dfx(i+1,j)-dfx(i,j))
415         &          +(dfy(i,j+1)-dfy(i,j)) )
416           ENDDO           ENDDO
417          ENDDO          ENDDO
418  C      end if GGL90diffTKEh .eq. 0.  C      end if GGL90diffTKEh .eq. 0.
# Line 396  C     vertical shear term (dU/dz)^2+(dV/ Line 425  C     vertical shear term (dU/dz)^2+(dV/
425           tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)           tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)
426       &                 -( uVel(i,j,k  ,bi,bj)+uVel(i+1,j,k  ,bi,bj)) )       &                 -( uVel(i,j,k  ,bi,bj)+uVel(i+1,j,k  ,bi,bj)) )
427       &        *recip_drC(k)       &        *recip_drC(k)
428    #ifdef ALLOW_GGL90_IDEMIX
429         &        *recip_hFacI(i,j,k,bi,bj)
430    #endif
431           tempV= .5 _d 0*( vVel(i,j,km1,bi,bj)+vVel(i,j+1,km1,bi,bj)           tempV= .5 _d 0*( vVel(i,j,km1,bi,bj)+vVel(i,j+1,km1,bi,bj)
432       &                 -( vVel(i,j,k  ,bi,bj)+vVel(i,j+1,k  ,bi,bj)) )       &                 -( vVel(i,j,k  ,bi,bj)+vVel(i,j+1,k  ,bi,bj)) )
433       &        *recip_drC(k)       &        *recip_drC(k)
434    #ifdef ALLOW_GGL90_IDEMIX
435         &        *recip_hFacI(i,j,k,bi,bj)
436    #endif
437           verticalShear = tempU*tempU + tempV*tempV           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  
438    
439  C     viscosity and diffusivity  C     viscosity and diffusivity
440           KappaM = GGL90ck*GGL90mixingLength(i,j,k)*SQRTTKE(i,j,k)           KappaM = GGL90ck*GGL90mixingLength(i,j,k)*SQRTTKE(i,j,k)
441           GGL90visctmp(i,j,k) = MAX(KappaM,diffKrNrT(k))           GGL90visctmp(i,j,k) = MAX(KappaM,diffKrNrT(k))
442       &                            * maskC(i,j,k,bi,bj)       &                            * maskC(i,j,k,bi,bj)
443  c        note: storing GGL90visctmp like this, and using it later to compute  C        note: storing GGL90visctmp like this, and using it later to compute
444  c              GGL9rdiffKr etc. is robust in case of smoothing (e.g. see OPA)  C              GGL9rdiffKr etc. is robust in case of smoothing (e.g. see OPA)
445           KappaM = MAX(KappaM,viscArNr(k)) * maskC(i,j,k,bi,bj)           KappaM = MAX(KappaM,viscArNr(k)) * maskC(i,j,k,bi,bj)
446    
447    C     compute Prandtl number (always greater than 0)
448             RiNumber = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)
449    #ifdef ALLOW_GGL90_IDEMIX
450    CML         IDEMIX_RiNumber = 1./GGL90eps
451             IDEMIX_RiNumber = MAX( KappaM*Nsquare(i,j,k), 0. _d 0)/
452         &    (GGL90eps+IDEMIX_tau_d(i,j,k,bi,bj)*IDEMIX_E(i,j,k,bi,bj)**2)
453             prTemp         = MIN(5.*RiNumber, 6.6*IDEMIX_RiNumber)
454    #else
455             prTemp = 1. _d 0
456             IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
457    #endif /* ALLOW_GGL90_IDEMIX */
458             TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
459             TKEPrandtlNumber(i,j,k) = MAX( 1. _d 0,TKEPrandtlNumber(i,j,k))
460    
461    C        diffusivity
462           KappaH = KappaM/TKEPrandtlNumber(i,j,k)           KappaH = KappaM/TKEPrandtlNumber(i,j,k)
463           KappaE(i,j,k) = GGL90alpha * KappaM * maskC(i,j,k,bi,bj)           KappaE(i,j,k) = GGL90alpha * KappaM * maskC(i,j,k,bi,bj)
464    
# Line 427  C     partial update with sum of explici Line 472  C     partial update with sum of explici
472       &        + KappaM*verticalShear       &        + KappaM*verticalShear
473       &        - KappaH*Nsquare(i,j,k)       &        - KappaH*Nsquare(i,j,k)
474       &        - TKEdissipation       &        - TKEdissipation
475    #ifdef ALLOW_GGL90_IDEMIX
476         &        + IDEMIX_tau_d(i,j,k,bi,bj)*IDEMIX_E(i,j,k,bi,bj)**2
477    #endif
478       &        )       &        )
479          ENDDO          ENDDO
480         ENDDO         ENDDO
# Line 454  C     set up matrix Line 502  C     set up matrix
502  C--   Lower diagonal  C--   Lower diagonal
503        DO j=jMin,jMax        DO j=jMin,jMax
504         DO i=iMin,iMax         DO i=iMin,iMax
505           a(i,j,1) = 0. _d 0           a3d(i,j,1) = 0. _d 0
506         ENDDO         ENDDO
507        ENDDO        ENDDO
508        DO k=2,Nr        DO k=2,Nr
# Line 464  C--   Lower diagonal Line 512  C--   Lower diagonal
512  C-    We keep recip_hFacC in the diffusive flux calculation,  C-    We keep recip_hFacC in the diffusive flux calculation,
513  C-    but no hFacC in TKE volume control  C-    but no hFacC in TKE volume control
514  C-    No need for maskC(k-1) with recip_hFacC(k-1)  C-    No need for maskC(k-1) with recip_hFacC(k-1)
515           a(i,j,k) = -deltaTggl90           a3d(i,j,k) = -deltaTggl90
516       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)
517       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))
518       &        *recip_drC(k)*maskC(i,j,k,bi,bj)       &        *recip_drC(k)*maskC(i,j,k,bi,bj)
519    #ifdef ALLOW_GGL90_IDEMIX
520         &        *recip_hFacI(i,j,k,bi,bj)
521    #endif
522          ENDDO          ENDDO
523         ENDDO         ENDDO
524        ENDDO        ENDDO
525  C--   Upper diagonal  C--   Upper diagonal
526        DO j=jMin,jMax        DO j=jMin,jMax
527         DO i=iMin,iMax         DO i=iMin,iMax
528           c(i,j,1)  = 0. _d 0           c3d(i,j,1)  = 0. _d 0
529         ENDDO         ENDDO
530        ENDDO        ENDDO
531        DO k=2,Nr        DO k=2,Nr
# Line 484  C--   Upper diagonal Line 535  C--   Upper diagonal
535  C-    We keep recip_hFacC in the diffusive flux calculation,  C-    We keep recip_hFacC in the diffusive flux calculation,
536  C-    but no hFacC in TKE volume control  C-    but no hFacC in TKE volume control
537  C-    No need for maskC(k) with recip_hFacC(k)  C-    No need for maskC(k) with recip_hFacC(k)
538            c(i,j,k) = -deltaTggl90            c3d(i,j,k) = -deltaTggl90
539       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)
540       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))
541       &        *recip_drC(k)*maskC(i,j,k-1,bi,bj)       &        *recip_drC(k)*maskC(i,j,k-1,bi,bj)
542    #ifdef ALLOW_GGL90_IDEMIX
543         &        *recip_hFacI(i,j,k,bi,bj)
544    #endif
545          ENDDO          ENDDO
546         ENDDO         ENDDO
547        ENDDO        ENDDO
548    
549          IF (.NOT.GGL90_dirichlet) THEN
550    C      Neumann bottom boundary condition for TKE: no flux from bottom
551           DO j=jMin,jMax
552            DO i=iMin,iMax
553             kBottom   = MAX(kLowC(i,j,bi,bj),1)
554             c3d(i,j,kBottom) = 0. _d 0
555            ENDDO
556           ENDDO
557          ENDIF
558    
559  C--   Center diagonal  C--   Center diagonal
560        DO k=1,Nr        DO k=1,Nr
561         km1 = MAX(k-1,1)         km1 = MAX(k-1,1)
562         DO j=jMin,jMax         DO j=jMin,jMax
563          DO i=iMin,iMax          DO i=iMin,iMax
564            b(i,j,k) = 1. _d 0 - c(i,j,k) - a(i,j,k)            b3d(i,j,k) = 1. _d 0 - c3d(i,j,k) - a3d(i,j,k)
565       &        + implDissFac*deltaTggl90*GGL90ceps*SQRTTKE(i,j,k)       &        + implDissFac*deltaTggl90*GGL90ceps*SQRTTKE(i,j,k)
566       &        * rMixingLength(i,j,k)       &        * rMixingLength(i,j,k)
567       &        * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)       &        * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)
# Line 520  C     Dirichlet surface boundary conditi Line 585  C     Dirichlet surface boundary conditi
585          GGL90TKE(i,j,kSurf,bi,bj) = maskC(i,j,kSurf,bi,bj)          GGL90TKE(i,j,kSurf,bi,bj) = maskC(i,j,kSurf,bi,bj)
586       &           *MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)       &           *MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)
587          GGL90TKE(i,j,kp1,bi,bj) = GGL90TKE(i,j,kp1,bi,bj)          GGL90TKE(i,j,kp1,bi,bj) = GGL90TKE(i,j,kp1,bi,bj)
588       &               - a(i,j,kp1)*GGL90TKE(i,j,kSurf,bi,bj)       &               - a3d(i,j,kp1)*GGL90TKE(i,j,kSurf,bi,bj)
589          a(i,j,kp1) = 0. _d 0          a3d(i,j,kp1) = 0. _d 0
 C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom  
         kBottom   = MAX(kLowC(i,j,bi,bj),1)  
         GGL90TKE(i,j,kBottom,bi,bj) = GGL90TKE(i,j,kBottom,bi,bj)  
      &                              - GGL90TKEbottom*c(i,j,kBottom)  
         c(i,j,kBottom) = 0. _d 0  
590         ENDDO         ENDDO
591        ENDDO        ENDDO
592    
593          IF (GGL90_dirichlet) THEN
594    C      Dirichlet bottom boundary condition for TKE = GGL90TKEbottom
595           DO j=jMin,jMax
596            DO i=iMin,iMax
597             kBottom   = MAX(kLowC(i,j,bi,bj),1)
598             GGL90TKE(i,j,kBottom,bi,bj) = GGL90TKE(i,j,kBottom,bi,bj)
599         &                              - GGL90TKEbottom*c3d(i,j,kBottom)
600             c3d(i,j,kBottom) = 0. _d 0
601            ENDDO
602           ENDDO
603          ENDIF
604    
605  C     solve tri-diagonal system  C     solve tri-diagonal system
606        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
607       I                        a, b, c,       I                        a3d, b3d, c3d,
608       U                        GGL90TKE,       U                        GGL90TKE(1-OLx,1-OLy,1,bi,bj),
609       O                        errCode,       O                        errCode,
610       I                        bi, bj, myThid )       I                        bi, bj, myThid )
611    
# Line 654  C     =============================== Line 726  C     ===============================
726    
727  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
728        IF ( useDiagnostics ) THEN        IF ( useDiagnostics ) THEN
729           CALL DIAGNOSTICS_FILL( GGL90TKE   ,'GGL90TKE',          CALL DIAGNOSTICS_FILL( GGL90TKE   ,'GGL90TKE',
730       &                          0,Nr, 1, bi, bj, myThid )       &                         0,Nr, 1, bi, bj, myThid )
731           CALL DIAGNOSTICS_FILL( GGL90viscArU,'GGL90ArU',          CALL DIAGNOSTICS_FILL( GGL90viscArU,'GGL90ArU',
732       &                          0,Nr, 1, bi, bj, myThid )       &                         0,Nr, 1, bi, bj, myThid )
733           CALL DIAGNOSTICS_FILL( GGL90viscArV,'GGL90ArV',          CALL DIAGNOSTICS_FILL( GGL90viscArV,'GGL90ArV',
734       &                          0,Nr, 1, bi, bj, myThid )       &                         0,Nr, 1, bi, bj, myThid )
735           CALL DIAGNOSTICS_FILL( GGL90diffKr,'GGL90Kr ',          CALL DIAGNOSTICS_FILL( GGL90diffKr,'GGL90Kr ',
736       &                          0,Nr, 1, bi, bj, myThid )       &                         0,Nr, 1, bi, bj, myThid )
737           CALL DIAGNOSTICS_FILL( TKEPrandtlNumber ,'GGL90Prl',          CALL DIAGNOSTICS_FILL( TKEPrandtlNumber ,'GGL90Prl',
738       &                          0,Nr, 2, bi, bj, myThid )       &                         0,Nr, 2, bi, bj, myThid )
739           CALL DIAGNOSTICS_FILL( GGL90mixingLength,'GGL90Lmx',          CALL DIAGNOSTICS_FILL( GGL90mixingLength,'GGL90Lmx',
740       &                          0,Nr, 2, bi, bj, myThid )       &                         0,Nr, 2, bi, bj, myThid )
741    
742            kp1 = MIN(Nr,kSurf+1)
743            DO j=jMin,jMax
744             DO i=iMin,iMax
745    C     diagnose surface flux of TKE
746              surf_flx_tke(i,j) =(GGL90TKE(i,j,kSurf,bi,bj)-
747         &                        GGL90TKE(i,j,kp1,bi,bj))
748         &        *recip_drF(kSurf)*recip_hFacC(i,j,kSurf,bi,bj)
749         &        *KappaE(i,j,kp1)
750             ENDDO
751            ENDDO
752            CALL DIAGNOSTICS_FILL( surf_flx_tke,'GGL90flx',
753         &                         0, 1, 2, bi, bj, myThid )
754    
755            k=kSurf
756            DO j=jMin,jMax
757             DO i=iMin,iMax
758    C     diagnose work done by the wind
759              surf_flx_tke(i,j) =
760         &      halfRL*( surfaceForcingU(i,  j,bi,bj)*uVel(i  ,j,k,bi,bj)
761         &              +surfaceForcingU(i+1,j,bi,bj)*uVel(i+1,j,k,bi,bj))
762         &    + halfRL*( surfaceForcingV(i,j,  bi,bj)*vVel(i,j  ,k,bi,bj)
763         &              +surfaceForcingV(i,j+1,bi,bj)*vVel(i,j+1,k,bi,bj))
764             ENDDO
765            ENDDO
766            CALL DIAGNOSTICS_FILL( surf_flx_tke,'GGL90tau',
767         &                         0, 1, 2, bi, bj, myThid )
768    
769        ENDIF        ENDIF
770  #endif  #endif /* ALLOW_DIAGNOSTICS */
771    
772  #endif /* ALLOW_GGL90 */  #endif /* ALLOW_GGL90 */
773    
Legend:
Removed from v.1.19  
changed lines
  Added in v.1.29
  ViewVC Help
Powered by ViewVC 1.1.22