/[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.20 by jmc, Thu Mar 15 15:23:22 2012 UTC
# Line 98  c     _RL     SQRTTKE Line 98  c     _RL     SQRTTKE
98        _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
99        _RL     GGL90visctmp     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     GGL90visctmp     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
100  C-    tri-diagonal matrix  C-    tri-diagonal matrix
101        _RL     a(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     a3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
102        _RL     b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     b3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
103        _RL     c(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     c3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
104        INTEGER errCode        INTEGER errCode
105  #ifdef ALLOW_GGL90_HORIZDIFF  #ifdef ALLOW_GGL90_HORIZDIFF
106  C     xA, yA   :: area of lateral faces  C     xA, yA   :: area of lateral faces
# Line 139  C     Initialize local fields Line 139  C     Initialize local fields
139           TKEPrandtlNumber(i,j,k)  = 1. _d 0           TKEPrandtlNumber(i,j,k)  = 1. _d 0
140           GGL90mixingLength(i,j,k) = GGL90mixingLengthMin           GGL90mixingLength(i,j,k) = GGL90mixingLengthMin
141           GGL90visctmp(i,j,k)      = 0. _d 0           GGL90visctmp(i,j,k)      = 0. _d 0
142    #ifndef SOLVE_DIAGONAL_LOWMEMORY
143             a3d(i,j,k) = 0. _d 0
144             b3d(i,j,k) = 1. _d 0
145             c3d(i,j,k) = 0. _d 0
146    #endif
147          ENDDO          ENDDO
148         ENDDO         ENDDO
149        ENDDO        ENDDO
# Line 454  C     set up matrix Line 459  C     set up matrix
459  C--   Lower diagonal  C--   Lower diagonal
460        DO j=jMin,jMax        DO j=jMin,jMax
461         DO i=iMin,iMax         DO i=iMin,iMax
462           a(i,j,1) = 0. _d 0           a3d(i,j,1) = 0. _d 0
463         ENDDO         ENDDO
464        ENDDO        ENDDO
465        DO k=2,Nr        DO k=2,Nr
# Line 464  C--   Lower diagonal Line 469  C--   Lower diagonal
469  C-    We keep recip_hFacC in the diffusive flux calculation,  C-    We keep recip_hFacC in the diffusive flux calculation,
470  C-    but no hFacC in TKE volume control  C-    but no hFacC in TKE volume control
471  C-    No need for maskC(k-1) with recip_hFacC(k-1)  C-    No need for maskC(k-1) with recip_hFacC(k-1)
472           a(i,j,k) = -deltaTggl90           a3d(i,j,k) = -deltaTggl90
473       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)
474       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))
475       &        *recip_drC(k)*maskC(i,j,k,bi,bj)       &        *recip_drC(k)*maskC(i,j,k,bi,bj)
# Line 474  C-    No need for maskC(k-1) with recip_ Line 479  C-    No need for maskC(k-1) with recip_
479  C--   Upper diagonal  C--   Upper diagonal
480        DO j=jMin,jMax        DO j=jMin,jMax
481         DO i=iMin,iMax         DO i=iMin,iMax
482           c(i,j,1)  = 0. _d 0           c3d(i,j,1)  = 0. _d 0
483         ENDDO         ENDDO
484        ENDDO        ENDDO
485        DO k=2,Nr        DO k=2,Nr
# Line 484  C--   Upper diagonal Line 489  C--   Upper diagonal
489  C-    We keep recip_hFacC in the diffusive flux calculation,  C-    We keep recip_hFacC in the diffusive flux calculation,
490  C-    but no hFacC in TKE volume control  C-    but no hFacC in TKE volume control
491  C-    No need for maskC(k) with recip_hFacC(k)  C-    No need for maskC(k) with recip_hFacC(k)
492            c(i,j,k) = -deltaTggl90            c3d(i,j,k) = -deltaTggl90
493       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)
494       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))
495       &        *recip_drC(k)*maskC(i,j,k-1,bi,bj)       &        *recip_drC(k)*maskC(i,j,k-1,bi,bj)
# Line 496  C--   Center diagonal Line 501  C--   Center diagonal
501         km1 = MAX(k-1,1)         km1 = MAX(k-1,1)
502         DO j=jMin,jMax         DO j=jMin,jMax
503          DO i=iMin,iMax          DO i=iMin,iMax
504            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)
505       &        + implDissFac*deltaTggl90*GGL90ceps*SQRTTKE(i,j,k)       &        + implDissFac*deltaTggl90*GGL90ceps*SQRTTKE(i,j,k)
506       &        * rMixingLength(i,j,k)       &        * rMixingLength(i,j,k)
507       &        * 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 525  C     Dirichlet surface boundary conditi
525          GGL90TKE(i,j,kSurf,bi,bj) = maskC(i,j,kSurf,bi,bj)          GGL90TKE(i,j,kSurf,bi,bj) = maskC(i,j,kSurf,bi,bj)
526       &           *MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)       &           *MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)
527          GGL90TKE(i,j,kp1,bi,bj) = GGL90TKE(i,j,kp1,bi,bj)          GGL90TKE(i,j,kp1,bi,bj) = GGL90TKE(i,j,kp1,bi,bj)
528       &               - a(i,j,kp1)*GGL90TKE(i,j,kSurf,bi,bj)       &               - a3d(i,j,kp1)*GGL90TKE(i,j,kSurf,bi,bj)
529          a(i,j,kp1) = 0. _d 0          a3d(i,j,kp1) = 0. _d 0
530  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom
531          kBottom   = MAX(kLowC(i,j,bi,bj),1)          kBottom   = MAX(kLowC(i,j,bi,bj),1)
532          GGL90TKE(i,j,kBottom,bi,bj) = GGL90TKE(i,j,kBottom,bi,bj)          GGL90TKE(i,j,kBottom,bi,bj) = GGL90TKE(i,j,kBottom,bi,bj)
533       &                              - GGL90TKEbottom*c(i,j,kBottom)       &                              - GGL90TKEbottom*c3d(i,j,kBottom)
534          c(i,j,kBottom) = 0. _d 0          c3d(i,j,kBottom) = 0. _d 0
535         ENDDO         ENDDO
536        ENDDO        ENDDO
537    
538  C     solve tri-diagonal system  C     solve tri-diagonal system
539        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
540       I                        a, b, c,       I                        a3d, b3d, c3d,
541       U                        GGL90TKE,       U                        GGL90TKE,
542       O                        errCode,       O                        errCode,
543       I                        bi, bj, myThid )       I                        bi, bj, myThid )
Legend:
Removed from v.1.19  
changed lines
  Added in v.1.20
  ViewVC Help
Powered by ViewVC 1.1.22