/[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.31 by jmc, Mon Feb 23 21:20:15 2015 UTC revision 1.32 by mlosch, Thu Feb 26 16:45:24 2015 UTC
# Line 59  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 77  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
# Line 101  c     _RL     SQRTTKE Line 103  c     _RL     SQRTTKE
103  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
104        _RL     surf_flx_tke     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     surf_flx_tke     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105  #endif /* ALLOW_DIAGNOSTICS */  #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)
115        _RL     c3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     c3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
116        INTEGER errCode        INTEGER errCode
117  #ifdef ALLOW_GGL90_HORIZDIFF  #ifdef ALLOW_GGL90_HORIZDIFF
 C     hFac     :: fractional thickness of W-cell  
118  C     xA, yA   :: area of lateral faces  C     xA, yA   :: area of lateral faces
119  C     dfx, dfy :: diffusive flux across lateral faces  C     dfx, dfy :: diffusive flux across lateral faces
120  C     gTKE     :: right hand side of diffusion equation  C     gTKE     :: right hand side of diffusion equation
       _RL     hFac  
121        _RL     xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
122        _RL     yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
123        _RL     dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 138  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
# Line 161  C     Initialize local fields Line 187  C     Initialize local fields
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
# Line 177  C     Initialize local fields Line 205  C     Initialize local fields
205    
206  #ifdef ALLOW_GGL90_IDEMIX  #ifdef ALLOW_GGL90_IDEMIX
207        IF ( useIDEMIX) CALL GGL90_IDEMIX(        IF ( useIDEMIX) CALL GGL90_IDEMIX(
208       &   bi, bj, sigmaR, myTime, myIter, myThid )       &   bi, bj, hFacI, recip_hFacI, sigmaR, myTime, myIter, myThid )
209  #endif /* ALLOW_GGL90_IDEMIX */  #endif /* ALLOW_GGL90_IDEMIX */
210    
 C     start k-loop  
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 191  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 220  C- ensure mixing between first and secon Line 234  C- ensure mixing between first and secon
234         ENDDO         ENDDO
235        ENDIF        ENDIF
236    
237  C- Impose upper and lower bound for mixing length  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 342  c         MaxLength=MAX(MaxLength,20. _d Line 357  c         MaxLength=MAX(MaxLength,20. _d
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    
# Line 401  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
 #ifdef ALLOW_GGL90_IDEMIX  
410            gTKE(i,j) = -recip_drC(k)*recip_rA(i,j,bi,bj)            gTKE(i,j) = -recip_drC(k)*recip_rA(i,j,bi,bj)
411       &         *recip_hFacI(i,j,k,bi,bj)       &         *recip_hFacI(i,j,k)
 #else  
           hFac = MIN(.5 _d 0,_hFacC(i,j,k-1,bi,bj) ) +  
      &           MIN(.5 _d 0,_hFacC(i,j,k  ,bi,bj) )  
           gTKE(i,j) = 0.0  
           IF ( hFac .ne. 0.0 )  
      &      gTKE(i,j) = -recip_drC(k)*recip_rA(i,j,bi,bj)/hFac  
 #endif  
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         DO j=jMin,jMax
422          DO i=iMin,iMax          DO i=iMin,iMax
423  C     vertical shear term (dU/dz)^2+(dV/dz)^2           KappaM(i,j) = GGL90ck*GGL90mixingLength(i,j,k)*SQRTTKE(i,j,k)
424           tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)           GGL90visctmp(i,j,k) = MAX(KappaM(i,j),diffKrNrS(k))
      &                 -( uVel(i,j,k  ,bi,bj)+uVel(i+1,j,k  ,bi,bj)) )  
      &        *recip_drC(k)  
 #ifdef ALLOW_GGL90_IDEMIX  
      &        *recip_hFacI(i,j,k,bi,bj)  
 #endif  
          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)  
 #ifdef ALLOW_GGL90_IDEMIX  
      &        *recip_hFacI(i,j,k,bi,bj)  
 #endif  
          verticalShear = tempU*tempU + tempV*tempV  
   
 C     viscosity and diffusivity  
          KappaM = GGL90ck*GGL90mixingLength(i,j,k)*SQRTTKE(i,j,k)  
          GGL90visctmp(i,j,k) = MAX(KappaM,diffKrNrS(k))  
425       &                            * maskC(i,j,k,bi,bj)       &                            * maskC(i,j,k,bi,bj)
426  C        note: storing GGL90visctmp like this, and using it later to compute  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)  C              GGL9rdiffKr etc. is robust in case of smoothing (e.g. see OPA)
428           KappaM = MAX(KappaM,viscArNr(k)) * maskC(i,j,k,bi,bj)           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)  C     compute Prandtl number (always greater than 0)
          RiNumber = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)  
433  #ifdef ALLOW_GGL90_IDEMIX  #ifdef ALLOW_GGL90_IDEMIX
434           IF ( useIDEMIX ) THEN
435           DO j=jMin,jMax
436            DO i=iMin,iMax
437    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)
439         &                 -( 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  CML         IDEMIX_RiNumber = 1./GGL90eps
448           IDEMIX_RiNumber = MAX( KappaM*Nsquare(i,j,k), 0. _d 0)/           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)       &    (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)           prTemp         = MIN(5.*RiNumber, 6.6 _d 0*IDEMIX_RiNumber)
451  #else           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)
464             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)) )
466         &        *recip_drC(k)
467             verticalShear(i,j) = tempU*tempU + tempV*tempV
468             RiNumber = MAX(Nsquare(i,j,k),0. _d 0)
469         &        /(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
 #endif /* ALLOW_GGL90_IDEMIX */  
472           TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)           TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
473           TKEPrandtlNumber(i,j,k) = MAX( 1. _d 0,TKEPrandtlNumber(i,j,k))          ENDDO
474           ENDDO
475           ENDIF
476    
477           DO j=jMin,jMax
478            DO i=iMin,iMax
479  C        diffusivity  C        diffusivity
480           KappaH = KappaM/TKEPrandtlNumber(i,j,k)           KappaH = KappaM(i,j)/TKEPrandtlNumber(i,j,k)
481           KappaE(i,j,k) = GGL90alpha * KappaM * maskC(i,j,k,bi,bj)           KappaE(i,j,k) = GGL90alpha * KappaM(i,j) * maskC(i,j,k,bi,bj)
482    
483  C     dissipation term  C     dissipation term
484           TKEdissipation = explDissFac*GGL90ceps           TKEdissipation = explDissFac*GGL90ceps
# Line 469  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
 #ifdef ALLOW_GGL90_IDEMIX  
      &        + IDEMIX_tau_d(i,j,k,bi,bj)*IDEMIX_E(i,j,k,bi,bj)**2  
 #endif  
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 516  C-    No need for maskC(k-1) with recip_ Line 545  C-    No need for maskC(k-1) with recip_
545       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)
546       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))
547       &        *recip_drC(k)*maskC(i,j,k,bi,bj)       &        *recip_drC(k)*maskC(i,j,k,bi,bj)
 #ifdef ALLOW_GGL90_IDEMIX  
      &        *recip_hFacI(i,j,k,bi,bj)  
 #endif  
548          ENDDO          ENDDO
549         ENDDO         ENDDO
550        ENDDO        ENDDO
# Line 531  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)
 #ifdef ALLOW_GGL90_IDEMIX  
      &        *recip_hFacI(i,j,k,bi,bj)  
 #endif  
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        IF (.NOT.GGL90_dirichlet) THEN
586  C      Neumann bottom boundary condition for TKE: no flux from bottom  C      Neumann bottom boundary condition for TKE: no flux from bottom
587         DO j=jMin,jMax         DO j=jMin,jMax
# Line 561  C--   Center diagonal Line 597  C--   Center diagonal
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
Legend:
Removed from v.1.31  
changed lines
  Added in v.1.32
  ViewVC Help
Powered by ViewVC 1.1.22