/[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.22 by jmc, Thu Jun 28 15:35:52 2012 UTC revision 1.32 by mlosch, Thu Feb 26 16:45:24 2015 UTC
# Line 10  C !INTERFACE: ========================== Line 10  C !INTERFACE: ==========================
10        SUBROUTINE GGL90_CALC(        SUBROUTINE GGL90_CALC(
11       I                 bi, bj, sigmaR, myTime, myIter, myThid )       I                 bi, bj, sigmaR, myTime, myIter, myThid )
12    
   
13  C !DESCRIPTION: \bv  C !DESCRIPTION: \bv
14  C     *==========================================================*  C     *==========================================================*
15  C     | SUBROUTINE GGL90_CALC                                    |  C     | SUBROUTINE GGL90_CALC                                    |
# Line 52  C     myThid :: My Thread Id number Line 51  C     myThid :: My Thread Id number
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 61  C Local constants Line 59  C Local constants
59  C     iMin,iMax,jMin,jMax :: index boundaries of computation domain  C     iMin,iMax,jMin,jMax :: index boundaries of computation domain
60  C     i, j, k, kp1,km1 :: array computation indices  C     i, j, k, kp1,km1 :: array computation indices
61  C     kSurf, kBottom   :: vertical indices of domain boundaries  C     kSurf, kBottom   :: vertical indices of domain boundaries
62    C     hFac/hFacI       :: fractional thickness of W-cell
63  C     explDissFac      :: explicit Dissipation Factor (in [0-1])  C     explDissFac      :: explicit Dissipation Factor (in [0-1])
64  C     implDissFac      :: implicit Dissipation Factor (in [0-1])  C     implDissFac      :: implicit Dissipation Factor (in [0-1])
65  C     uStarSquare      :: square of friction velocity  C     uStarSquare      :: square of friction velocity
# Line 79  C     TKEPrandtlNumber :: here, an empir Line 78  C     TKEPrandtlNumber :: here, an empir
78        INTEGER i, j, k, kp1, km1, kSurf, kBottom        INTEGER i, j, k, kp1, km1, kSurf, kBottom
79        _RL     explDissFac, implDissFac        _RL     explDissFac, implDissFac
80        _RL     uStarSquare        _RL     uStarSquare
81        _RL     verticalShear        _RL     verticalShear(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
82        _RL     KappaM, KappaH        _RL     KappaM(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
83          _RL     KappaH
84  c     _RL     Nsquare  c     _RL     Nsquare
85        _RL     Nsquare(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     Nsquare(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
86        _RL     deltaTggl90        _RL     deltaTggl90
87  c     _RL     SQRTTKE  c     _RL     SQRTTKE
88        _RL     SQRTTKE(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     SQRTTKE(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
89        _RL     RiNumber        _RL     RiNumber
90    #ifdef ALLOW_GGL90_IDEMIX
91          _RL     IDEMIX_RiNumber
92    #endif
93        _RL     TKEdissipation        _RL     TKEdissipation
94        _RL     tempU, tempV, prTemp        _RL     tempU, tempV, prTemp
95        _RL     MaxLength, tmpmlx, tmpVisc        _RL     MaxLength, tmpmlx, tmpVisc
# Line 97  c     _RL     SQRTTKE Line 100  c     _RL     SQRTTKE
100        _RL     KappaE           (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     KappaE           (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
101        _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102        _RL     GGL90visctmp     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     GGL90visctmp     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
103    #ifdef ALLOW_DIAGNOSTICS
104          _RL     surf_flx_tke     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105    #endif /* ALLOW_DIAGNOSTICS */
106    C     hFac(I)  :: fractional thickness of W-cell
107          _RL       hFac
108    #ifdef ALLOW_GGL90_IDEMIX
109          _RL       hFacI(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
110    #endif /* ALLOW_GGL90_IDEMIX */
111          _RL recip_hFacI(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
112  C-    tri-diagonal matrix  C-    tri-diagonal matrix
113        _RL     a3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     a3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
114        _RL     b3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     b3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
# Line 114  C     gTKE     :: right hand side of dif Line 126  C     gTKE     :: right hand side of dif
126  #endif /* ALLOW_GGL90_HORIZDIFF */  #endif /* ALLOW_GGL90_HORIZDIFF */
127  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
128        _RL p4, p8, p16        _RL p4, p8, p16
129    CEOP
130        p4=0.25 _d 0        p4=0.25 _d 0
131        p8=0.125 _d 0        p8=0.125 _d 0
132        p16=0.0625 _d 0        p16=0.0625 _d 0
# Line 131  C     explicit/implicit timestepping wei Line 144  C     explicit/implicit timestepping wei
144        explDissFac = 0. _d 0        explDissFac = 0. _d 0
145        implDissFac = 1. _d 0 - explDissFac        implDissFac = 1. _d 0 - explDissFac
146    
147    C     For nonlinear free surface and especially with r*-coordinates, the
148    C     hFacs change every timestep, so we need to update them here in the
149    C     case of using IDEMIX.
150           DO K=1,Nr
151            Km1 = MAX(K-1,1)
152            DO j=1-OLy,sNy+OLy
153             DO i=1-OLx,sNx+OLx
154              hFac =
155         &         MIN(.5 _d 0,_hFacC(i,j,km1,bi,bj) ) +
156         &         MIN(.5 _d 0,_hFacC(i,j,k  ,bi,bj) )
157              recip_hFacI(I,J,K)=0. _d 0
158              IF ( hFac .NE. 0. _d 0 )
159         &         recip_hFacI(I,J,K)=1. _d 0/hFac
160    #ifdef ALLOW_GGL90_IDEMIX
161              hFacI(i,j,k) = hFac
162    #endif /* ALLOW_GGL90_IDEMIX */
163             ENDDO
164            ENDDO
165           ENDDO
166    
167  C     Initialize local fields  C     Initialize local fields
168        DO k = 1, Nr        DO k = 1, Nr
169         DO j=1-OLy,sNy+OLy         DO j=1-OLy,sNy+OLy
170          DO i=1-OLx,sNx+OLx          DO i=1-OLx,sNx+OLx
171             rMixingLength(i,j,k)     = 0. _d 0
172             mxLength_Dn(i,j,k)       = 0. _d 0
173             GGL90visctmp(i,j,k)      = 0. _d 0
174           KappaE(i,j,k)            = 0. _d 0           KappaE(i,j,k)            = 0. _d 0
175           TKEPrandtlNumber(i,j,k)  = 1. _d 0           TKEPrandtlNumber(i,j,k)  = 1. _d 0
176           GGL90mixingLength(i,j,k) = GGL90mixingLengthMin           GGL90mixingLength(i,j,k) = GGL90mixingLengthMin
# Line 144  C     Initialize local fields Line 180  C     Initialize local fields
180           b3d(i,j,k) = 1. _d 0           b3d(i,j,k) = 1. _d 0
181           c3d(i,j,k) = 0. _d 0           c3d(i,j,k) = 0. _d 0
182  #endif  #endif
183             Nsquare(i,j,k) = 0. _d 0
184             SQRTTKE(i,j,k) = 0. _d 0
185          ENDDO          ENDDO
186         ENDDO         ENDDO
187        ENDDO        ENDDO
188        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
189         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
190            KappaM(i,j)        = 0. _d 0
191            verticalShear(i,j) = 0. _d 0
192          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)
193          rMixingLength(i,j,1) = 0. _d 0          rMixingLength(i,j,1) = 0. _d 0
194          mxLength_Dn(i,j,1) = GGL90mixingLengthMin          mxLength_Dn(i,j,1) = GGL90mixingLengthMin
195          SQRTTKE(i,j,1) = SQRT( GGL90TKE(i,j,1,bi,bj) )          SQRTTKE(i,j,1) = SQRT( GGL90TKE(i,j,1,bi,bj) )
196    #ifdef ALLOW_GGL90_HORIZDIFF
197            xA(i,j)  = 0. _d 0
198            yA(i,j)  = 0. _d 0
199            dfx(i,j) = 0. _d 0
200            dfy(i,j) = 0. _d 0
201            gTKE(i,j) = 0. _d 0
202    #endif /* ALLOW_GGL90_HORIZDIFF */
203         ENDDO         ENDDO
204        ENDDO        ENDDO
205    
206  C     start k-loop  #ifdef ALLOW_GGL90_IDEMIX
207          IF ( useIDEMIX) CALL GGL90_IDEMIX(
208         &   bi, bj, hFacI, recip_hFacI, sigmaR, myTime, myIter, myThid )
209    #endif /* ALLOW_GGL90_IDEMIX */
210    
211        DO k = 2, Nr        DO k = 2, Nr
 c      km1 = k-1  
 c      kp1 = MIN(Nr,k+1)  
212         DO j=jMin,jMax         DO j=jMin,jMax
213          DO i=iMin,iMax          DO i=iMin,iMax
214           SQRTTKE(i,j,k)=SQRT( GGL90TKE(i,j,k,bi,bj) )           SQRTTKE(i,j,k)=SQRT( GGL90TKE(i,j,k,bi,bj) )
# Line 167  c      kp1 = MIN(Nr,k+1) Line 216  c      kp1 = MIN(Nr,k+1)
216  C     buoyancy frequency  C     buoyancy frequency
217           Nsquare(i,j,k) = gravity*gravitySign*recip_rhoConst           Nsquare(i,j,k) = gravity*gravitySign*recip_rhoConst
218       &                  * sigmaR(i,j,k)       &                  * sigmaR(i,j,k)
219  cC     vertical shear term (dU/dz)^2+(dV/dz)^2  C     vertical shear term (dU/dz)^2+(dV/dz)^2 is computed later
220  c         tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)  C     to save some memory
 c     &                 -( uVel(i,j,k  ,bi,bj)+uVel(i+1,j,k  ,bi,bj)) )  
 c     &        *recip_drC(k)  
 c         tempV= .5 _d 0*( vVel(i,j,km1,bi,bj)+vVel(i,j+1,km1,bi,bj)  
 c     &                 -( vVel(i,j,k  ,bi,bj)+vVel(i,j+1,k  ,bi,bj)) )  
 c     &        *recip_drC(k)  
 c         verticalShear = tempU*tempU + tempV*tempV  
 c         RiNumber   = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)  
 cC     compute Prandtl number (always greater than 0)  
 c         prTemp = 1. _d 0  
 c         IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber  
 c         TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)  
221  C     mixing length  C     mixing length
222           GGL90mixingLength(i,j,k) = SQRTTWO *           GGL90mixingLength(i,j,k) = SQRTTWO *
223       &        SQRTTKE(i,j,k)/SQRT( MAX(Nsquare(i,j,k),GGL90eps) )       &        SQRTTKE(i,j,k)/SQRT( MAX(Nsquare(i,j,k),GGL90eps) )
# Line 187  C     mixing length Line 225  C     mixing length
225         ENDDO         ENDDO
226        ENDDO        ENDDO
227    
228  C- Impose upper and lower bound for mixing length  C- ensure mixing between first and second level
229          IF (mxlSurfFlag) THEN
230           DO j=jMin,jMax
231            DO i=iMin,iMax
232             GGL90mixingLength(i,j,2)=drF(1)
233            ENDDO
234           ENDDO
235          ENDIF
236    
237    C--   Impose upper and lower bound for mixing length
238    C--   Impose minimum mixing length to avoid division by zero
239        IF ( mxlMaxFlag .EQ. 0 ) THEN        IF ( mxlMaxFlag .EQ. 0 ) THEN
240    
241         DO k=2,Nr         DO k=2,Nr
# Line 235  c         MaxLength=MAX(MaxLength,20. _d Line 283  c         MaxLength=MAX(MaxLength,20. _d
283    
284        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN
285    
 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  
286         DO k=2,Nr         DO k=2,Nr
287          DO j=jMin,jMax          DO j=jMin,jMax
288           DO i=iMin,iMax           DO i=iMin,iMax
# Line 316  cgf Line 357  cgf
357         STOP 'GGL90_CALC: Wrong mxlMaxFlag (mixing length limit)'         STOP 'GGL90_CALC: Wrong mxlMaxFlag (mixing length limit)'
358        ENDIF        ENDIF
359    
360  C- Impose minimum mixing length (to avoid division by zero)  C     start "proper" k-loop (the code above was moved out and up to
361  c      DO k=2,Nr  C     implemement various mixing parameters efficiently)
 c      DO j=jMin,jMax  
 c       DO i=iMin,iMax  
 c        GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),  
 c    &        GGL90mixingLengthMin)  
 c        rMixingLength(i,j,k) = 1. _d 0 /GGL90mixingLength(i,j,k)  
 c       ENDDO  
 c      ENDDO  
 c     ENDDO  
   
362        DO k=2,Nr        DO k=2,Nr
363         km1 = k-1         km1 = k-1
364    
365  #ifdef ALLOW_GGL90_HORIZDIFF  #ifdef ALLOW_GGL90_HORIZDIFF
366         IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN        IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
367  C     horizontal diffusion of TKE (requires an exchange in  C     horizontal diffusion of TKE (requires an exchange in
368  C      do_fields_blocking_exchanges)  C      do_fields_blocking_exchanges)
369  C     common factors  C     common factors
370          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
371           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
372            xA(i,j) = _dyG(i,j,bi,bj)            xA(i,j) = _dyG(i,j,bi,bj)*drC(k)*
373       &         *drF(k)*_hFacW(i,j,k,bi,bj)       &                 (min(.5 _d 0,_hFacW(i,j,k-1,bi,bj) ) +
374            yA(i,j) = _dxG(i,j,bi,bj)       &                  min(.5 _d 0,_hFacW(i,j,k  ,bi,bj) ) )
375       &         *drF(k)*_hFacS(i,j,k,bi,bj)            yA(i,j) = _dxG(i,j,bi,bj)*drC(k)*
376         &                 (min(.5 _d 0,_hFacS(i,j,k-1,bi,bj) ) +
377         &                  min(.5 _d 0,_hFacS(i,j,k  ,bi,bj) ) )
378           ENDDO           ENDDO
379          ENDDO          ENDDO
380  C     Compute diffusive fluxes  C     Compute diffusive fluxes
# Line 351  C     ... across x-faces Line 385  C     ... across x-faces
385            dfx(i,j) = -GGL90diffTKEh*xA(i,j)            dfx(i,j) = -GGL90diffTKEh*xA(i,j)
386       &      *_recip_dxC(i,j,bi,bj)       &      *_recip_dxC(i,j,bi,bj)
387       &      *(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))
388    #ifdef ISOTROPIC_COS_SCALING
389       &      *CosFacU(j,bi,bj)       &      *CosFacU(j,bi,bj)
390    #endif /* ISOTROPIC_COS_SCALING */
391           ENDDO           ENDDO
392          ENDDO          ENDDO
393  C     ... across y-faces  C     ... across y-faces
# Line 371  C     ... across y-faces Line 407  C     ... across y-faces
407  C     Compute divergence of fluxes  C     Compute divergence of fluxes
408          DO j=1-OLy,sNy+OLy-1          DO j=1-OLy,sNy+OLy-1
409           DO i=1-OLx,sNx+OLx-1           DO i=1-OLx,sNx+OLx-1
410            gTKE(i,j) =            gTKE(i,j) = -recip_drC(k)*recip_rA(i,j,bi,bj)
411       &    -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj)       &         *recip_hFacI(i,j,k)
412       &         *( (dfx(i+1,j)-dfx(i,j))       &         *((dfx(i+1,j)-dfx(i,j))
413       &           +(dfy(i,j+1)-dfy(i,j))       &         + (dfy(i,j+1)-dfy(i,j)) )
      &          )  
414           ENDDO           ENDDO
415          ENDDO          ENDDO
416  C      end if GGL90diffTKEh .eq. 0.  C     end if GGL90diffTKEh .eq. 0.
417         ENDIF         ENDIF
418  #endif /* ALLOW_GGL90_HORIZDIFF */  #endif /* ALLOW_GGL90_HORIZDIFF */
419    
420    C     viscosity and diffusivity
421           DO j=jMin,jMax
422            DO i=iMin,iMax
423             KappaM(i,j) = GGL90ck*GGL90mixingLength(i,j,k)*SQRTTKE(i,j,k)
424             GGL90visctmp(i,j,k) = MAX(KappaM(i,j),diffKrNrS(k))
425         &                            * maskC(i,j,k,bi,bj)
426    C        note: storing GGL90visctmp like this, and using it later to compute
427    C              GGL9rdiffKr etc. is robust in case of smoothing (e.g. see OPA)
428             KappaM(i,j) = MAX(KappaM(i,j),viscArNr(k)) * maskC(i,j,k,bi,bj)
429            ENDDO
430           ENDDO
431    
432    C     compute Prandtl number (always greater than 0)
433    #ifdef ALLOW_GGL90_IDEMIX
434           IF ( useIDEMIX ) THEN
435         DO j=jMin,jMax         DO j=jMin,jMax
436          DO i=iMin,iMax          DO i=iMin,iMax
437  C     vertical shear term (dU/dz)^2+(dV/dz)^2  C     vertical shear term (dU/dz)^2+(dV/dz)^2
438           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)
439       &                 -( 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)) )
440         &        *recip_drC(k)*recip_hFacI(i,j,k)
441             tempV= .5 _d 0*( vVel(i,j,km1,bi,bj)+vVel(i,j+1,km1,bi,bj)
442         &                 -( vVel(i,j,k  ,bi,bj)+vVel(i,j+1,k  ,bi,bj)) )
443         &        *recip_drC(k)*recip_hFacI(i,j,k)
444             verticalShear(i,j) = tempU*tempU + tempV*tempV
445             RiNumber = MAX(Nsquare(i,j,k),0. _d 0)
446         &        /(verticalShear(i,j)+GGL90eps)
447    CML         IDEMIX_RiNumber = 1./GGL90eps
448             IDEMIX_RiNumber = MAX( KappaM(i,j)*Nsquare(i,j,k), 0. _d 0)/
449         &    (GGL90eps+IDEMIX_tau_d(i,j,k,bi,bj)*IDEMIX_E(i,j,k,bi,bj)**2)
450             prTemp         = MIN(5.*RiNumber, 6.6 _d 0*IDEMIX_RiNumber)
451             TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
452             TKEPrandtlNumber(i,j,k) = MAX( 1. _d 0,TKEPrandtlNumber(i,j,k))
453            ENDDO
454           ENDDO
455           ELSE
456    #else /* ndef ALLOW_GGL90_IDEMIX */
457           IF (.TRUE.) THEN
458    #endif /* ALLOW_GGL90_IDEMIX */
459           DO j=jMin,jMax
460            DO i=iMin,iMax
461             tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)
462         &                 -( uVel(i,j,k  ,bi,bj)+uVel(i+1,j,k  ,bi,bj)) )
463       &        *recip_drC(k)       &        *recip_drC(k)
464           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)
465       &                 -( 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)) )
466       &        *recip_drC(k)       &        *recip_drC(k)
467           verticalShear = tempU*tempU + tempV*tempV           verticalShear(i,j) = tempU*tempU + tempV*tempV
468           RiNumber = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)           RiNumber = MAX(Nsquare(i,j,k),0. _d 0)
469  C     compute Prandtl number (always greater than 0)       &        /(verticalShear(i,j)+GGL90eps)
470           prTemp = 1. _d 0           prTemp = 1. _d 0
471           IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber           IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
472           TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)           TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
473  c         TKEPrandtlNumber(i,j,k) = 1. _d 0          ENDDO
474           ENDDO
475           ENDIF
476    
477  C     viscosity and diffusivity         DO j=jMin,jMax
478           KappaM = GGL90ck*GGL90mixingLength(i,j,k)*SQRTTKE(i,j,k)          DO i=iMin,iMax
479           GGL90visctmp(i,j,k) = MAX(KappaM,diffKrNrT(k))  C        diffusivity
480       &                            * maskC(i,j,k,bi,bj)           KappaH = KappaM(i,j)/TKEPrandtlNumber(i,j,k)
481  c        note: storing GGL90visctmp like this, and using it later to compute           KappaE(i,j,k) = GGL90alpha * KappaM(i,j) * maskC(i,j,k,bi,bj)
 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)  
482    
483  C     dissipation term  C     dissipation term
484           TKEdissipation = explDissFac*GGL90ceps           TKEdissipation = explDissFac*GGL90ceps
# Line 416  C     dissipation term Line 487  C     dissipation term
487  C     partial update with sum of explicit contributions  C     partial update with sum of explicit contributions
488           GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)           GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
489       &        + deltaTggl90*(       &        + deltaTggl90*(
490       &        + KappaM*verticalShear       &        + KappaM(i,j)*verticalShear(i,j)
491       &        - KappaH*Nsquare(i,j,k)       &        - KappaH*Nsquare(i,j,k)
492       &        - TKEdissipation       &        - TKEdissipation
493       &        )       &        )
494          ENDDO          ENDDO
495         ENDDO         ENDDO
496    
497    #ifdef ALLOW_GGL90_IDEMIX
498           IF ( useIDEMIX ) THEN
499    C     add IDEMIX contribution to the turbulent kinetic energy
500            DO j=jMin,jMax
501             DO i=iMin,iMax
502              GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
503         &         + deltaTggl90*(
504         &         + IDEMIX_tau_d(i,j,k,bi,bj)*IDEMIX_E(i,j,k,bi,bj)**2
505         &         )
506             ENDDO
507            ENDDO
508           ENDIF
509    #endif /* ALLOW_GGL90_IDEMIX */
510    
511  #ifdef ALLOW_GGL90_HORIZDIFF  #ifdef ALLOW_GGL90_HORIZDIFF
512         IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN         IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
513  C--    Add horiz. diffusion tendency  C--    Add horiz. diffusion tendency
# Line 472  C--   Upper diagonal Line 557  C--   Upper diagonal
557        DO k=2,Nr        DO k=2,Nr
558         DO j=jMin,jMax         DO j=jMin,jMax
559          DO i=iMin,iMax          DO i=iMin,iMax
560            kp1=MAX(1,MIN(klowC(i,j,bi,bj),k+1))           kp1=MAX(1,MIN(klowC(i,j,bi,bj),k+1))
561  C-    We keep recip_hFacC in the diffusive flux calculation,  C-    We keep recip_hFacC in the diffusive flux calculation,
562  C-    but no hFacC in TKE volume control  C-    but no hFacC in TKE volume control
563  C-    No need for maskC(k) with recip_hFacC(k)  C-    No need for maskC(k) with recip_hFacC(k)
564            c3d(i,j,k) = -deltaTggl90           c3d(i,j,k) = -deltaTggl90
565       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)
566       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))
567       &        *recip_drC(k)*maskC(i,j,k-1,bi,bj)       &        *recip_drC(k)*maskC(i,j,k-1,bi,bj)
568          ENDDO          ENDDO
569         ENDDO         ENDDO
570        ENDDO        ENDDO
571    
572    #ifdef ALLOW_GGL90_IDEMIX
573          IF ( useIDEMIX ) THEN
574           DO k=2,Nr
575            DO j=jMin,jMax
576             DO i=iMin,iMax
577              a3d(i,j,k) = a3d(i,j,k)*recip_hFacI(i,j,k)
578              c3d(i,j,k) = c3d(i,j,k)*recip_hFacI(i,j,k)
579             ENDDO
580            ENDDO
581           ENDDO
582          ENDIF
583    #endif /* ALLOW_GGL90_IDEMIX */
584    
585          IF (.NOT.GGL90_dirichlet) THEN
586    C      Neumann bottom boundary condition for TKE: no flux from bottom
587           DO j=jMin,jMax
588            DO i=iMin,iMax
589             kBottom   = MAX(kLowC(i,j,bi,bj),1)
590             c3d(i,j,kBottom) = 0. _d 0
591            ENDDO
592           ENDDO
593          ENDIF
594    
595  C--   Center diagonal  C--   Center diagonal
596        DO k=1,Nr        DO k=1,Nr
597         km1 = MAX(k-1,1)         km1 = MAX(k-1,1)
598         DO j=jMin,jMax         DO j=jMin,jMax
599          DO i=iMin,iMax          DO i=iMin,iMax
600            b3d(i,j,k) = 1. _d 0 - c3d(i,j,k) - a3d(i,j,k)           b3d(i,j,k) = 1. _d 0 - c3d(i,j,k) - a3d(i,j,k)
601       &        + implDissFac*deltaTggl90*GGL90ceps*SQRTTKE(i,j,k)       &        + implDissFac*deltaTggl90*GGL90ceps*SQRTTKE(i,j,k)
602       &        * rMixingLength(i,j,k)       &        * rMixingLength(i,j,k)
603       &        * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)       &        * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)
604           ENDDO          ENDDO
605         ENDDO         ENDDO
606        ENDDO        ENDDO
607  C     end set up matrix  C     end set up matrix
# Line 514  C     Dirichlet surface boundary conditi Line 623  C     Dirichlet surface boundary conditi
623          GGL90TKE(i,j,kp1,bi,bj) = GGL90TKE(i,j,kp1,bi,bj)          GGL90TKE(i,j,kp1,bi,bj) = GGL90TKE(i,j,kp1,bi,bj)
624       &               - a3d(i,j,kp1)*GGL90TKE(i,j,kSurf,bi,bj)       &               - a3d(i,j,kp1)*GGL90TKE(i,j,kSurf,bi,bj)
625          a3d(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*c3d(i,j,kBottom)  
         c3d(i,j,kBottom) = 0. _d 0  
626         ENDDO         ENDDO
627        ENDDO        ENDDO
628    
629          IF (GGL90_dirichlet) THEN
630    C      Dirichlet bottom boundary condition for TKE = GGL90TKEbottom
631           DO j=jMin,jMax
632            DO i=iMin,iMax
633             kBottom   = MAX(kLowC(i,j,bi,bj),1)
634             GGL90TKE(i,j,kBottom,bi,bj) = GGL90TKE(i,j,kBottom,bi,bj)
635         &                              - GGL90TKEbottom*c3d(i,j,kBottom)
636             c3d(i,j,kBottom) = 0. _d 0
637            ENDDO
638           ENDDO
639          ENDIF
640    
641  C     solve tri-diagonal system  C     solve tri-diagonal system
642        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
643       I                        a3d, b3d, c3d,       I                        a3d, b3d, c3d,
644       U                        GGL90TKE,       U                        GGL90TKE(1-OLx,1-OLy,1,bi,bj),
645       O                        errCode,       O                        errCode,
646       I                        bi, bj, myThid )       I                        bi, bj, myThid )
647    
# Line 572  C     =============================== Line 688  C     ===============================
688           tmpVisc = GGL90visctmp(i,j,k)           tmpVisc = GGL90visctmp(i,j,k)
689  #endif  #endif
690           tmpVisc = MIN(tmpVisc/TKEPrandtlNumber(i,j,k),GGL90diffMax)           tmpVisc = MIN(tmpVisc/TKEPrandtlNumber(i,j,k),GGL90diffMax)
691           GGL90diffKr(i,j,k,bi,bj)= MAX( tmpVisc , diffKrNrT(k) )           GGL90diffKr(i,j,k,bi,bj)= MAX( tmpVisc , diffKrNrS(k) )
692          ENDDO          ENDDO
693         ENDDO         ENDDO
694        ENDDO        ENDDO
# Line 646  C     =============================== Line 762  C     ===============================
762    
763  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
764        IF ( useDiagnostics ) THEN        IF ( useDiagnostics ) THEN
765           CALL DIAGNOSTICS_FILL( GGL90TKE   ,'GGL90TKE',          CALL DIAGNOSTICS_FILL( GGL90TKE   ,'GGL90TKE',
766       &                          0,Nr, 1, bi, bj, myThid )       &                         0,Nr, 1, bi, bj, myThid )
767           CALL DIAGNOSTICS_FILL( GGL90viscArU,'GGL90ArU',          CALL DIAGNOSTICS_FILL( GGL90viscArU,'GGL90ArU',
768       &                          0,Nr, 1, bi, bj, myThid )       &                         0,Nr, 1, bi, bj, myThid )
769           CALL DIAGNOSTICS_FILL( GGL90viscArV,'GGL90ArV',          CALL DIAGNOSTICS_FILL( GGL90viscArV,'GGL90ArV',
770       &                          0,Nr, 1, bi, bj, myThid )       &                         0,Nr, 1, bi, bj, myThid )
771           CALL DIAGNOSTICS_FILL( GGL90diffKr,'GGL90Kr ',          CALL DIAGNOSTICS_FILL( GGL90diffKr,'GGL90Kr ',
772       &                          0,Nr, 1, bi, bj, myThid )       &                         0,Nr, 1, bi, bj, myThid )
773           CALL DIAGNOSTICS_FILL( TKEPrandtlNumber ,'GGL90Prl',          CALL DIAGNOSTICS_FILL( TKEPrandtlNumber ,'GGL90Prl',
774       &                          0,Nr, 2, bi, bj, myThid )       &                         0,Nr, 2, bi, bj, myThid )
775           CALL DIAGNOSTICS_FILL( GGL90mixingLength,'GGL90Lmx',          CALL DIAGNOSTICS_FILL( GGL90mixingLength,'GGL90Lmx',
776       &                          0,Nr, 2, bi, bj, myThid )       &                         0,Nr, 2, bi, bj, myThid )
777    
778            kp1 = MIN(Nr,kSurf+1)
779            DO j=jMin,jMax
780             DO i=iMin,iMax
781    C     diagnose surface flux of TKE
782              surf_flx_tke(i,j) =(GGL90TKE(i,j,kSurf,bi,bj)-
783         &                        GGL90TKE(i,j,kp1,bi,bj))
784         &        *recip_drF(kSurf)*recip_hFacC(i,j,kSurf,bi,bj)
785         &        *KappaE(i,j,kp1)
786             ENDDO
787            ENDDO
788            CALL DIAGNOSTICS_FILL( surf_flx_tke,'GGL90flx',
789         &                         0, 1, 2, bi, bj, myThid )
790    
791            k=kSurf
792            DO j=jMin,jMax
793             DO i=iMin,iMax
794    C     diagnose work done by the wind
795              surf_flx_tke(i,j) =
796         &      halfRL*( surfaceForcingU(i,  j,bi,bj)*uVel(i  ,j,k,bi,bj)
797         &              +surfaceForcingU(i+1,j,bi,bj)*uVel(i+1,j,k,bi,bj))
798         &    + halfRL*( surfaceForcingV(i,j,  bi,bj)*vVel(i,j  ,k,bi,bj)
799         &              +surfaceForcingV(i,j+1,bi,bj)*vVel(i,j+1,k,bi,bj))
800             ENDDO
801            ENDDO
802            CALL DIAGNOSTICS_FILL( surf_flx_tke,'GGL90tau',
803         &                         0, 1, 2, bi, bj, myThid )
804    
805        ENDIF        ENDIF
806  #endif  #endif /* ALLOW_DIAGNOSTICS */
807    
808  #endif /* ALLOW_GGL90 */  #endif /* ALLOW_GGL90 */
809    
Legend:
Removed from v.1.22  
changed lines
  Added in v.1.32
  ViewVC Help
Powered by ViewVC 1.1.22