/[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.33 by heimbach, Thu Jul 2 04:42:43 2015 UTC revision 1.34 by jmc, Thu Jan 14 17:50:25 2016 UTC
# Line 91  c     _RL     SQRTTKE Line 91  c     _RL     SQRTTKE
91        _RL     IDEMIX_RiNumber        _RL     IDEMIX_RiNumber
92  #endif  #endif
93        _RL     TKEdissipation        _RL     TKEdissipation
94        _RL     tempU, tempV, prTemp        _RL     tempU, tempUp, tempV, tempVp, prTemp
95        _RL     MaxLength, tmpmlx, tmpVisc        _RL     MaxLength, tmpmlx, tmpVisc
96        _RL     TKEPrandtlNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     TKEPrandtlNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
97        _RL     GGL90mixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     GGL90mixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
# Line 126  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    #endif
130  CEOP  CEOP
131        p4=0.25 _d 0  
132        p8=0.125 _d 0        PARAMETER( iMin = 2-OLx, iMax = sNx+OLx-1 )
133        p16=0.0625 _d 0        PARAMETER( jMin = 2-OLy, jMax = sNy+OLy-1 )
134    #ifdef ALLOW_GGL90_SMOOTH
135          p4  = 0.25   _d 0
136          p8  = 0.125  _d 0
137          p16 = 0.0625 _d 0
138  #endif  #endif
       iMin = 2-OLx  
       iMax = sNx+OLx-1  
       jMin = 2-OLy  
       jMax = sNy+OLy-1  
139    
140  C     set separate time step (should be deltaTtracer)  C     set separate time step (should be deltaTtracer)
141        deltaTggl90 = dTtracerLev(1)        deltaTggl90 = dTtracerLev(1)
# Line 144  C     explicit/implicit timestepping wei Line 145  C     explicit/implicit timestepping wei
145        explDissFac = 0. _d 0        explDissFac = 0. _d 0
146        implDissFac = 1. _d 0 - explDissFac        implDissFac = 1. _d 0 - explDissFac
147    
148  C     For nonlinear free surface and especially with r*-coordinates, the  C     For nonlinear free surface and especially with r*-coordinates, the
149  C     hFacs change every timestep, so we need to update them here in the  C     hFacs change every timestep, so we need to update them here in the
150  C     case of using IDEMIX.  C     case of using IDEMIX.
151         DO K=1,Nr         DO K=1,Nr
152          Km1 = MAX(K-1,1)          Km1 = MAX(K-1,1)
153          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
154           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
155            hFac =            hFac =
156       &         MIN(.5 _d 0,_hFacC(i,j,km1,bi,bj) ) +       &         MIN(.5 _d 0,_hFacC(i,j,km1,bi,bj) ) +
157       &         MIN(.5 _d 0,_hFacC(i,j,k  ,bi,bj) )       &         MIN(.5 _d 0,_hFacC(i,j,k  ,bi,bj) )
158            recip_hFacI(I,J,K)=0. _d 0            recip_hFacI(I,J,K)=0. _d 0
159            IF ( hFac .NE. 0. _d 0 )            IF ( hFac .NE. 0. _d 0 )
160       &         recip_hFacI(I,J,K)=1. _d 0/hFac       &         recip_hFacI(I,J,K)=1. _d 0/hFac
# Line 429  C              GGL9rdiffKr etc. is robus Line 430  C              GGL9rdiffKr etc. is robus
430          ENDDO          ENDDO
431         ENDDO         ENDDO
432    
433    C     compute vertical shear (dU/dz)^2+(dV/dz)^2
434           IF ( calcMeanVertShear ) THEN
435    C     by averaging (@ grid-cell center) the 4 vertical shear compon @ U,V pos.
436            DO j=jMin,jMax
437             DO i=iMin,iMax
438              tempU  = ( uVel( i ,j,km1,bi,bj) - uVel( i ,j,k,bi,bj) )
439              tempUp = ( uVel(i+1,j,km1,bi,bj) - uVel(i+1,j,k,bi,bj) )
440              tempV  = ( vVel(i, j ,km1,bi,bj) - vVel(i, j ,k,bi,bj) )
441              tempVp = ( vVel(i,j+1,km1,bi,bj) - vVel(i,j+1,k,bi,bj) )
442              verticalShear(i,j) = (
443         &                 ( tempU*tempU + tempUp*tempUp )*halfRL
444         &               + ( tempV*tempV + tempVp*tempVp )*halfRL
445         &                         )*recip_drC(k)*recip_drC(k)
446             ENDDO
447            ENDDO
448           ELSE
449    C     from the averaged flow at grid-cell center (2 compon x 2 pos.)
450            DO j=jMin,jMax
451             DO i=iMin,iMax
452              tempU = ( ( uVel(i,j,km1,bi,bj) + uVel(i+1,j,km1,bi,bj) )
453         &             -( uVel(i,j,k  ,bi,bj) + uVel(i+1,j,k  ,bi,bj) )
454         &            )*halfRL*recip_drC(k)
455              tempV = ( ( vVel(i,j,km1,bi,bj) + vVel(i,j+1,km1,bi,bj) )
456         &             -( vVel(i,j,k  ,bi,bj) + vVel(i,j+1,k  ,bi,bj) )
457         &            )*halfRL*recip_drC(k)
458              verticalShear(i,j) = tempU*tempU + tempV*tempV
459             ENDDO
460            ENDDO
461           ENDIF
462    
463  C     compute Prandtl number (always greater than 0)  C     compute Prandtl number (always greater than 0)
464  #ifdef ALLOW_GGL90_IDEMIX  #ifdef ALLOW_GGL90_IDEMIX
465         IF ( useIDEMIX ) THEN         IF ( useIDEMIX ) THEN
466         DO j=jMin,jMax          DO j=jMin,jMax
467          DO i=iMin,iMax           DO i=iMin,iMax
468  C     vertical shear term (dU/dz)^2+(dV/dz)^2  C     account for partical cell factor in vertical shear:
469           tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)            verticalShear(i,j) = verticalShear(i,j)
470       &                 -( uVel(i,j,k  ,bi,bj)+uVel(i+1,j,k  ,bi,bj)) )       &                       * recip_hFacI(i,j,k)*recip_hFacI(i,j,k)
471       &        *recip_drC(k)*recip_hFacI(i,j,k)            RiNumber = MAX(Nsquare(i,j,k),0. _d 0)
472           tempV= .5 _d 0*( vVel(i,j,km1,bi,bj)+vVel(i,j+1,km1,bi,bj)       &         /(verticalShear(i,j)+GGL90eps)
      &                 -( vVel(i,j,k  ,bi,bj)+vVel(i,j+1,k  ,bi,bj)) )  
      &        *recip_drC(k)*recip_hFacI(i,j,k)  
          verticalShear(i,j) = tempU*tempU + tempV*tempV  
          RiNumber = MAX(Nsquare(i,j,k),0. _d 0)  
      &        /(verticalShear(i,j)+GGL90eps)  
473  CML         IDEMIX_RiNumber = 1./GGL90eps  CML         IDEMIX_RiNumber = 1./GGL90eps
474           IDEMIX_RiNumber = MAX( KappaM(i,j)*Nsquare(i,j,k), 0. _d 0)/            IDEMIX_RiNumber = MAX( KappaM(i,j)*Nsquare(i,j,k), 0. _d 0)/
475       &    (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)
476           prTemp         = MIN(5.*RiNumber, 6.6 _d 0*IDEMIX_RiNumber)            prTemp         = MIN(5.*RiNumber, 6.6 _d 0*IDEMIX_RiNumber)
477           TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)            TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
478           TKEPrandtlNumber(i,j,k) = MAX( 1. _d 0,TKEPrandtlNumber(i,j,k))            TKEPrandtlNumber(i,j,k) = MAX( 1. _d 0,TKEPrandtlNumber(i,j,k))
479             ENDDO
480          ENDDO          ENDDO
        ENDDO  
481         ELSE         ELSE
482  #else /* ndef ALLOW_GGL90_IDEMIX */  #else /* ndef ALLOW_GGL90_IDEMIX */
483         IF (.TRUE.) THEN         IF (.TRUE.) THEN
484  #endif /* ALLOW_GGL90_IDEMIX */  #endif /* ALLOW_GGL90_IDEMIX */
485         DO j=jMin,jMax          DO j=jMin,jMax
486          DO i=iMin,iMax           DO i=iMin,iMax
487           tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)            RiNumber = MAX(Nsquare(i,j,k),0. _d 0)
488       &                 -( uVel(i,j,k  ,bi,bj)+uVel(i+1,j,k  ,bi,bj)) )       &         /(verticalShear(i,j)+GGL90eps)
489       &        *recip_drC(k)            prTemp = 1. _d 0
490           tempV= .5 _d 0*( vVel(i,j,km1,bi,bj)+vVel(i,j+1,km1,bi,bj)            IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
491       &                 -( vVel(i,j,k  ,bi,bj)+vVel(i,j+1,k  ,bi,bj)) )            TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
492       &        *recip_drC(k)           ENDDO
          verticalShear(i,j) = tempU*tempU + tempV*tempV  
          RiNumber = MAX(Nsquare(i,j,k),0. _d 0)  
      &        /(verticalShear(i,j)+GGL90eps)  
          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)  
493          ENDDO          ENDDO
        ENDDO  
494         ENDIF         ENDIF
495    
496         DO j=jMin,jMax         DO j=jMin,jMax
# Line 662  C     =============================== Line 681  C     ===============================
681         DO j=1,sNy         DO j=1,sNy
682          DO i=1,sNx          DO i=1,sNx
683  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
684           tmpVisc=           tmpVisc = (
685       &  (       &     p4 *    GGL90visctmp(i  ,j  ,k)*mskCor(i  ,j  ,bi,bj)
686       &   p4 *  GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)       &    +p8 *( ( GGL90visctmp(i-1,j  ,k)*mskCor(i-1,j  ,bi,bj)
687       &  +p8 *( GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj)       &           + GGL90visctmp(i+1,j  ,k)*mskCor(i+1,j  ,bi,bj) )
688       &       + GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj)       &         + ( GGL90visctmp(i  ,j-1,k)*mskCor(i  ,j-1,bi,bj)
689       &       + GGL90visctmp(i+1,j  ,k) * mskCor(i+1,j  ,bi,bj)       &           + GGL90visctmp(i  ,j+1,k)*mskCor(i  ,j+1,bi,bj) ) )
690       &       + GGL90visctmp(i  ,j+1,k) * mskCor(i  ,j+1,bi,bj))       &    +p16*( ( GGL90visctmp(i+1,j+1,k)*mskCor(i+1,j+1,bi,bj)
691       &  +p16*( GGL90visctmp(i+1,j+1,k) * mskCor(i+1,j+1,bi,bj)       &           + GGL90visctmp(i-1,j-1,k)*mskCor(i-1,j-1,bi,bj) )
692       &       + GGL90visctmp(i+1,j-1,k) * mskCor(i+1,j-1,bi,bj)       &         + ( GGL90visctmp(i+1,j-1,k)*mskCor(i+1,j-1,bi,bj)
693       &       + GGL90visctmp(i-1,j+1,k) * mskCor(i-1,j+1,bi,bj)       &           + GGL90visctmp(i-1,j+1,k)*mskCor(i-1,j+1,bi,bj) ) )
694       &       + GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj))       &             )/(
695       &  )       &     p4
696       & /(p4       &    +p8 *( (  maskC(i-1,j  ,k,bi,bj)*mskCor(i-1,j  ,bi,bj)
697       &  +p8 *(       maskC(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)       &           +  maskC(i+1,j  ,k,bi,bj)*mskCor(i+1,j  ,bi,bj) )
698       &       +       maskC(i  ,j-1,k,bi,bj) * mskCor(i  ,j-1,bi,bj)       &         + (  maskC(i  ,j-1,k,bi,bj)*mskCor(i  ,j-1,bi,bj)
699       &       +       maskC(i+1,j  ,k,bi,bj) * mskCor(i+1,j  ,bi,bj)       &           +  maskC(i  ,j+1,k,bi,bj)*mskCor(i  ,j+1,bi,bj) ) )
700       &       +       maskC(i  ,j+1,k,bi,bj) * mskCor(i  ,j+1,bi,bj))       &    +p16*( (  maskC(i+1,j+1,k,bi,bj)* mskCor(i+1,j+1,bi,bj)
701       &  +p16*(       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) )
702       &       +       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)
703       &       +       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) ) )
704       &       +       maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj))       &               )*maskC(i,j,k,bi,bj)*mskCor(i,j,bi,bj)
      &  )*maskC(i,j,k,bi,bj)*mskCor(i,j,bi,bj)  
705  #else  #else
706           tmpVisc = GGL90visctmp(i,j,k)           tmpVisc = GGL90visctmp(i,j,k)
707  #endif  #endif
# Line 697  C     =============================== Line 715  C     ===============================
715         DO j=1,sNy         DO j=1,sNy
716          DO i=1,sNx+1          DO i=1,sNx+1
717  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
718          tmpVisc =           tmpVisc = (
719       & (       &     p4 *(   GGL90visctmp(i-1,j  ,k)*mskCor(i-1,j  ,bi,bj)
720       &   p4 *(GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)       &           + GGL90visctmp(i  ,j  ,k)*mskCor(i  ,j  ,bi,bj) )
721       &       +GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj))       &    +p8 *( ( GGL90visctmp(i-1,j-1,k)*mskCor(i-1,j-1,bi,bj)
722       &  +p8 *(GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj)       &           + GGL90visctmp(i  ,j-1,k)*mskCor(i  ,j-1,bi,bj) )
723       &       +GGL90visctmp(i-1,j+1,k) * mskCor(i-1,j+1,bi,bj)       &         + ( GGL90visctmp(i-1,j+1,k)*mskCor(i-1,j+1,bi,bj)
724       &       +GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj)       &           + GGL90visctmp(i  ,j+1,k)*mskCor(i  ,j+1,bi,bj) ) )
725       &       +GGL90visctmp(i  ,j+1,k) * mskCor(i  ,j+1,bi,bj))       &             )/(
726       &  )       &     p4 * 2. _d 0
727       & /(p4 * 2. _d 0       &    +p8 *( (  maskC(i-1,j-1,k,bi,bj)*mskCor(i-1,j-1,bi,bj)
728       &  +p8 *(      maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj)       &           +  maskC(i  ,j-1,k,bi,bj)*mskCor(i  ,j-1,bi,bj) )
729       &       +      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)
730       &       +      maskC(i  ,j-1,k,bi,bj) * mskCor(i  ,j-1,bi,bj)       &           +  maskC(i  ,j+1,k,bi,bj)*mskCor(i  ,j+1,bi,bj) ) )
731       &       +      maskC(i  ,j+1,k,bi,bj) * mskCor(i  ,j+1,bi,bj))       &               )*maskC(i-1,j,k,bi,bj)*mskCor(i-1,j,bi,bj)
732       &  )       &                *maskC(i  ,j,k,bi,bj)*mskCor(i  ,j,bi,bj)
      &  *maskC(i  ,j,k,bi,bj)*mskCor(i  ,j,bi,bj)  
      &  *maskC(i-1,j,k,bi,bj)*mskCor(i-1,j,bi,bj)  
733  #else  #else
734          tmpVisc = _maskW(i,j,k,bi,bj) *           tmpVisc = _maskW(i,j,k,bi,bj) * halfRL
735       &                   (.5 _d 0*(GGL90visctmp(i,j,k)       &          *( GGL90visctmp(i-1,j,k)
736       &                            +GGL90visctmp(i-1,j,k))       &           + GGL90visctmp(i,j,k) )
      &                   )  
737  #endif  #endif
738          tmpVisc = MIN( tmpVisc , GGL90viscMax )           tmpVisc = MIN( tmpVisc , GGL90viscMax )
739          GGL90viscArU(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )           GGL90viscArU(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
740          ENDDO          ENDDO
741         ENDDO         ENDDO
742        ENDDO        ENDDO
# Line 730  C     =============================== Line 745  C     ===============================
745         DO j=1,sNy+1         DO j=1,sNy+1
746          DO i=1,sNx          DO i=1,sNx
747  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
748          tmpVisc =           tmpVisc = (
749       & (       &     p4 *(   GGL90visctmp(i  ,j-1,k)*mskCor(i  ,j-1,bi,bj)
750       &   p4 *(GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)       &           + GGL90visctmp(i  ,j  ,k)*mskCor(i  ,j  ,bi,bj) )
751       &       +GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj))       &    +p8 *( ( GGL90visctmp(i-1,j-1,k)*mskCor(i-1,j-1,bi,bj)
752       &  +p8 *(GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj)       &           + GGL90visctmp(i-1,j  ,k)*mskCor(i-1,j  ,bi,bj) )
753       &       +GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj)       &         + ( GGL90visctmp(i+1,j-1,k)*mskCor(i+1,j-1,bi,bj)
754       &       +GGL90visctmp(i+1,j  ,k) * mskCor(i+1,j  ,bi,bj)       &           + GGL90visctmp(i+1,j  ,k)*mskCor(i+1,j  ,bi,bj) ) )
755       &       +GGL90visctmp(i+1,j-1,k) * mskCor(i+1,j-1,bi,bj))       &             )/(
756       &  )       &     p4 * 2. _d 0
757       & /(p4 * 2. _d 0       &    +p8 *( (  maskC(i-1,j-1,k,bi,bj)*mskCor(i-1,j-1,bi,bj)
758       &  +p8 *(      maskC(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)       &           +  maskC(i-1,j  ,k,bi,bj)*mskCor(i-1,j  ,bi,bj) )
759       &       +      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)
760       &       +      maskC(i+1,j  ,k,bi,bj) * mskCor(i+1,j  ,bi,bj)       &           +  maskC(i+1,j  ,k,bi,bj)*mskCor(i+1,j  ,bi,bj) ) )
761       &       +      maskC(i+1,j-1,k,bi,bj) * mskCor(i+1,j-1,bi,bj))       &               )*maskC(i,j-1,k,bi,bj)*mskCor(i,j-1,bi,bj)
762       &  )       &                *maskC(i,j  ,k,bi,bj)*mskCor(i,j  ,bi,bj)
      &   *maskC(i,j  ,k,bi,bj)*mskCor(i,j  ,bi,bj)  
      &   *maskC(i,j-1,k,bi,bj)*mskCor(i,j-1,bi,bj)  
763  #else  #else
764          tmpVisc = _maskS(i,j,k,bi,bj) *           tmpVisc = _maskS(i,j,k,bi,bj) * halfRL
765       &                   (.5 _d 0*(GGL90visctmp(i,j,k)       &          *( GGL90visctmp(i,j-1,k)
766       &                            +GGL90visctmp(i,j-1,k))       &           + GGL90visctmp(i,j,k) )
      &                   )  
   
767  #endif  #endif
768          tmpVisc = MIN( tmpVisc , GGL90viscMax )           tmpVisc = MIN( tmpVisc , GGL90viscMax )
769          GGL90viscArV(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )           GGL90viscArV(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
770          ENDDO          ENDDO
771         ENDDO         ENDDO
772        ENDDO        ENDDO
Legend:
Removed from v.1.33  
changed lines
  Added in v.1.34
  ViewVC Help
Powered by ViewVC 1.1.22