/[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.15 by dfer, Sat Nov 14 14:27:56 2009 UTC revision 1.16 by gforget, Fri Aug 6 18:37:05 2010 UTC
# Line 84  c     _RL     SQRTTKE Line 84  c     _RL     SQRTTKE
84        _RL     RiNumber        _RL     RiNumber
85        _RL     TKEdissipation        _RL     TKEdissipation
86        _RL     tempU, tempV, prTemp        _RL     tempU, tempV, prTemp
87        _RL     MaxLength, tmpmlx        _RL     MaxLength, tmpmlx, tmpVisc
88        _RL     TKEPrandtlNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     TKEPrandtlNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
89        _RL     GGL90mixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     GGL90mixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
90        _RL     rMixingLength    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     rMixingLength    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
# Line 94  c     _RL     SQRTTKE Line 94  c     _RL     SQRTTKE
94        _RL     rhoKm1           (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     rhoKm1           (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
95        _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
96        _RL     gTKE             (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     gTKE             (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
97          _RL     GGL90visctmp     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
98  C-    tri-diagonal matrix  C-    tri-diagonal matrix
99        _RL     a(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     a(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
100        _RL     b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
# Line 108  C-    dfx, dfy - diffusive flux across l Line 109  C-    dfx, dfy - diffusive flux across l
109        _RL     dfy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     dfy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
110  #endif /* ALLOW_GGL90_HORIZDIFF */  #endif /* ALLOW_GGL90_HORIZDIFF */
111  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
112        _RL p4, p8, p16, tmpdiffKrS        _RL p4, p8, p16
113        p4=0.25 _d 0        p4=0.25 _d 0
114        p8=0.125 _d 0        p8=0.125 _d 0
115        p16=0.0625 _d 0        p16=0.0625 _d 0
# Line 240  c         MaxLength=MAX(MaxLength,20. _d Line 241  c         MaxLength=MAX(MaxLength,20. _d
241    
242        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN
243  C-  C-
244    cgf ensure mixing between first and second level
245    c      DO J=jMin,jMax
246    c        DO I=iMin,iMax
247    c         GGL90mixingLength(I,J,2)=drF(1)
248    c        ENDDO
249    c      ENDDO
250    cgf
251         DO k=2,Nr         DO k=2,Nr
252          DO J=jMin,jMax          DO J=jMin,jMax
253           DO I=iMin,iMax           DO I=iMin,iMax
# Line 325  c       ENDDO Line 333  c       ENDDO
333  c      ENDDO  c      ENDDO
334  c     ENDDO  c     ENDDO
335    
336    
337        DO k=2,Nr        DO k=2,Nr
338         Km1 = K-1         Km1 = K-1
339         DO J=jMin,jMax         DO J=jMin,jMax
# Line 346  c         TKEPrandtlNumber(I,J,K) = 1. _ Line 355  c         TKEPrandtlNumber(I,J,K) = 1. _
355    
356  C     viscosity and diffusivity  C     viscosity and diffusivity
357           KappaM = GGL90ck*GGL90mixingLength(I,J,K)*SQRTTKE(i,j,k)           KappaM = GGL90ck*GGL90mixingLength(I,J,K)*SQRTTKE(i,j,k)
358             GGL90visctmp(I,J,K) = MAX(KappaM,diffKrNrT(k))
359         &                            * maskC(I,J,K,bi,bj)
360    c        note: storing GGL90visctmp like this, and using it later to compute
361    c              GGL9rdiffKr etc. is robust in case of smoothing (e.g. see OPA)
362             KappaM = MAX(KappaM,viscArNr(k)) * maskC(I,J,K,bi,bj)
363           KappaH = KappaM/TKEPrandtlNumber(I,J,K)           KappaH = KappaM/TKEPrandtlNumber(I,J,K)
364             KappaE(I,J,K) = GGL90alpha * KappaM * maskC(I,J,K,bi,bj)
 C     Set a minium (= background) and maximum value  
          KappaM = MAX(KappaM,viscArNr(k))  
          KappaH = MAX(KappaH,diffKrNrT(k))  
          KappaM = MIN(KappaM,GGL90viscMax)  
          KappaH = MIN(KappaH,GGL90diffMax)  
   
 C     Mask land points and storage  
          GGL90viscAr(I,J,K,bi,bj) = KappaM * maskC(I,J,K,bi,bj)  
          GGL90diffKr(I,J,K,bi,bj) = KappaH * maskC(I,J,K,bi,bj)  
          KappaE(I,J,K) = GGL90alpha * GGL90viscAr(I,J,K,bi,bj)  
365    
366  C     dissipation term  C     dissipation term
367           TKEdissipation = ab05*GGL90ceps           TKEdissipation = ab05*GGL90ceps
# Line 519  C     solve tri-diagonal system, and sto Line 523  C     solve tri-diagonal system, and sto
523       U                        gTKE,       U                        gTKE,
524       O                        errCode,       O                        errCode,
525       I                        bi, bj, myThid )       I                        bi, bj, myThid )
 c     CALL GGL90_SOLVE( bi, bj, iMin, iMax, jMin, jMax,  
 c    I     a, b, c,  
 c    U     gTKE,  
 c    I     myThid )  
526    
527  C     now update TKE  C     now update TKE
528        DO K=1,Nr        DO K=1,Nr
# Line 538  C     impose minimum TKE to avoid numeri Line 538  C     impose minimum TKE to avoid numeri
538  C     end of time step  C     end of time step
539  C     ===============================  C     ===============================
540    
 #ifdef ALLOW_GGL90_SMOOTH  
541        DO K=1,Nr        DO K=1,Nr
542         DO J=jMin,jMax         DO J=jMin,jMax
543          DO I=iMin,iMax          DO I=iMin,iMax
544           tmpdiffKrS=  #ifdef ALLOW_GGL90_SMOOTH
545             tmpVisc=
546       &  (       &  (
547       &   p4 *  GGL90viscAr(i  ,j  ,k,bi,bj) * mskCor(i  ,j  ,bi,bj)       &   p4 *  GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
548       &  +p8 *( GGL90viscAr(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)       &  +p8 *( GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj)
549       &       + GGL90viscAr(i  ,j-1,k,bi,bj) * mskCor(i  ,j-1,bi,bj)       &       + GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj)
550       &       + GGL90viscAr(i+1,j  ,k,bi,bj) * mskCor(i+1,j  ,bi,bj)       &       + GGL90visctmp(i+1,j  ,k) * mskCor(i+1,j  ,bi,bj)
551       &       + GGL90viscAr(i  ,j+1,k,bi,bj) * mskCor(i  ,j+1,bi,bj))       &       + GGL90visctmp(i  ,j+1,k) * mskCor(i  ,j+1,bi,bj))
552       &  +p16*( GGL90viscAr(i+1,j+1,k,bi,bj) * mskCor(i+1,j+1,bi,bj)       &  +p16*( GGL90visctmp(i+1,j+1,k) * mskCor(i+1,j+1,bi,bj)
553       &       + GGL90viscAr(i+1,j-1,k,bi,bj) * mskCor(i+1,j-1,bi,bj)       &       + GGL90visctmp(i+1,j-1,k) * mskCor(i+1,j-1,bi,bj)
554       &       + GGL90viscAr(i-1,j+1,k,bi,bj) * mskCor(i-1,j+1,bi,bj)       &       + GGL90visctmp(i-1,j+1,k) * mskCor(i-1,j+1,bi,bj)
555       &       + GGL90viscAr(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj))       &       + GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj))
556       &  )       &  )
557       & /(p4       & /(p4
558       &  +p8 *(       maskC(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)       &  +p8 *(       maskC(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)
# Line 564  C     =============================== Line 564  C     ===============================
564       &       +       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)
565       &       +       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))
566       &  )*maskC(i,j,k,bi,bj)*mskCor(i,j,bi,bj)       &  )*maskC(i,j,k,bi,bj)*mskCor(i,j,bi,bj)
567       &   /TKEPrandtlNumber(i,j,k)  #else
568           GGL90diffKrS(I,J,K,bi,bj)= MAX( tmpdiffKrS , diffKrNrT(k) )           tmpVisc = GGL90visctmp(I,J,K)
569    #endif
570             tmpVisc = MIN(tmpVisc/TKEPrandtlNumber(i,j,k),GGL90diffMax)
571             GGL90diffKr(I,J,K,bi,bj)= MAX( tmpVisc , diffKrNrT(k) )
572          ENDDO          ENDDO
573         ENDDO         ENDDO
574        ENDDO        ENDDO
575    
576    
577    
578          DO K=1,Nr
579           DO J=jMin,jMax
580            DO I=iMin,iMax
581    #ifdef ALLOW_GGL90_SMOOTH
582            tmpVisc =
583         & (
584         &   p4 *(GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
585         &       +GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj))
586         &  +p8 *(GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj)
587         &       +GGL90visctmp(i-1,j+1,k) * mskCor(i-1,j+1,bi,bj)
588         &       +GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj)
589         &       +GGL90visctmp(i  ,j+1,k) * mskCor(i  ,j+1,bi,bj))
590         &  )
591         & /(p4 * 2. _d 0
592         &  +p8 *(      maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj)
593         &       +      maskC(i-1,j+1,k,bi,bj) * mskCor(i-1,j+1,bi,bj)
594         &       +      maskC(i  ,j-1,k,bi,bj) * mskCor(i  ,j-1,bi,bj)
595         &       +      maskC(i  ,j+1,k,bi,bj) * mskCor(i  ,j+1,bi,bj))
596         &  )
597         &  *maskC(i  ,j,k,bi,bj)*mskCor(i  ,j,bi,bj)
598         &  *maskC(i-1,j,k,bi,bj)*mskCor(i-1,j,bi,bj)
599    #else
600            tmpVisc = _maskW(i,j,k,bi,bj) *
601         &                   (.5 _d 0*(GGL90visctmp(i,j,k)
602         &                            +GGL90visctmp(i-1,j,k))
603         &                   )
604  #endif  #endif
605            tmpVisc = MIN( tmpVisc , GGL90viscMax )
606            GGL90viscArU(i,j,k,bi,bj) = MAX( tmpVisc , viscArNr(k)  )
607            ENDDO
608           ENDDO
609          ENDDO
610    
611    
612          DO K=1,Nr
613           DO J=jMin,jMax
614            DO I=iMin,iMax
615    #ifdef ALLOW_GGL90_SMOOTH
616            tmpVisc =
617         & (
618         &   p4 *(GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
619         &       +GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj))
620         &  +p8 *(GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj)
621         &       +GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj)
622         &       +GGL90visctmp(i+1,j  ,k) * mskCor(i+1,j  ,bi,bj)
623         &       +GGL90visctmp(i+1,j-1,k) * mskCor(i+1,j-1,bi,bj))
624         &  )
625         & /(p4 * 2. _d 0
626         &  +p8 *(      maskC(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)
627         &       +      maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj)
628         &       +      maskC(i+1,j  ,k,bi,bj) * mskCor(i+1,j  ,bi,bj)
629         &       +      maskC(i+1,j-1,k,bi,bj) * mskCor(i+1,j-1,bi,bj))
630         &  )
631         &   *maskC(i,j  ,k,bi,bj)*mskCor(i,j  ,bi,bj)
632         &   *maskC(i,j-1,k,bi,bj)*mskCor(i,j-1,bi,bj)
633    #else
634            tmpVisc = _maskS(i,j,k,bi,bj) *
635         &                   (.5 _d 0*(GGL90visctmp(i,j,k)
636         &                            +GGL90visctmp(i,j-1,k))
637         &                   )
638    
639    #endif
640            tmpVisc = MIN( tmpVisc , GGL90viscMax )
641            GGL90viscArV(i,j,k,bi,bj) = MAX( tmpVisc , viscArNr(k)  )
642            ENDDO
643           ENDDO
644          ENDDO
645    
646  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
647        IF ( useDiagnostics ) THEN        IF ( useDiagnostics ) THEN
648           CALL DIAGNOSTICS_FILL( GGL90TKE   ,'GGL90TKE',           CALL DIAGNOSTICS_FILL( GGL90TKE   ,'GGL90TKE',
649       &                          0,Nr, 1, bi, bj, myThid )       &                          0,Nr, 1, bi, bj, myThid )
650           CALL DIAGNOSTICS_FILL( GGL90viscAr,'GGL90Ar ',           CALL DIAGNOSTICS_FILL( GGL90viscArU,'GGL90ArU',
651         &                          0,Nr, 1, bi, bj, myThid )
652             CALL DIAGNOSTICS_FILL( GGL90viscArV,'GGL90ArV',
653       &                          0,Nr, 1, bi, bj, myThid )       &                          0,Nr, 1, bi, bj, myThid )
654           CALL DIAGNOSTICS_FILL( GGL90diffKr,'GGL90Kr ',           CALL DIAGNOSTICS_FILL( GGL90diffKr,'GGL90Kr ',
655       &                          0,Nr, 1, bi, bj, myThid )       &                          0,Nr, 1, bi, bj, myThid )
Legend:
Removed from v.1.15  
changed lines
  Added in v.1.16
  ViewVC Help
Powered by ViewVC 1.1.22