/[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.14 by jmc, Sat Oct 31 03:18:50 2009 UTC revision 1.15 by dfer, Sat Nov 14 14:27:56 2009 UTC
# Line 98  C-    tri-diagonal matrix Line 98  C-    tri-diagonal matrix
98        _RL     a(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     a(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
99        _RL     b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
100        _RL     c(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     c(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
101  c     INTEGER errCode        INTEGER errCode
102  #ifdef ALLOW_GGL90_HORIZDIFF  #ifdef ALLOW_GGL90_HORIZDIFF
103  C-    xA, yA   - area of lateral faces  C-    xA, yA   - area of lateral faces
104  C-    dfx, dfy - diffusive flux across lateral faces  C-    dfx, dfy - diffusive flux across lateral faces
# Line 414  C     ... across y-faces Line 414  C     ... across y-faces
414           ENDDO           ENDDO
415          ENDDO          ENDDO
416  C     Compute divergence of fluxes  C     Compute divergence of fluxes
 C- jmc: concerned about missing a deltaT since gTKE is already the future TKE.  
417          DO j=1-Oly,sNy+Oly-1          DO j=1-Oly,sNy+Oly-1
418           DO i=1-Olx,sNx+Olx-1           DO i=1-Olx,sNx+Olx-1
419            gTKE(i,j,k)=gTKE(i,j,k)            gTKE(i,j,k)=gTKE(i,j,k)
420       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj)       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj)
421       &         *( (dfx(i+1,j)-dfx(i,j))       &         *( (dfx(i+1,j)-dfx(i,j))
422       &           +(dfy(i,j+1)-dfy(i,j))       &           +(dfy(i,j+1)-dfy(i,j))
423       &           )       &           )*deltaTggl90
424           ENDDO           ENDDO
425          ENDDO          ENDDO
426  C       end of k-loop  C       end of k-loop
# Line 442  C--   Lower diagonal Line 441  C--   Lower diagonal
441         ENDDO         ENDDO
442        ENDDO        ENDDO
443        DO k=2,Nr        DO k=2,Nr
444  C- jmc: concerned that a(k=2) should always be zero         km1=MAX(2,k-1)
 C       and would be better without recip_hFacC  
        km1=max(2,k-1)  
445         DO j=jMin,jMax         DO j=jMin,jMax
446          DO i=iMin,iMax          DO i=iMin,iMax
447    C-    We keep recip_hFacC in the diffusive flux calculation,
448    C-    but no hFacC in TKE volume control
449    C-    No need for maskC(k-1) with recip_hFacC(k-1)
450           a(i,j,k) = -deltaTggl90           a(i,j,k) = -deltaTggl90
 c     &        *recip_drF(km1)*recip_hFacI(i,j,k,bi,bj)  
451       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)
452       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))
453       &        *recip_drC(k)*maskC(i,j,k,bi,bj)*maskC(i,j,k-1,bi,bj)       &        *recip_drC(k)*maskC(i,j,k,bi,bj)
454          ENDDO          ENDDO
455         ENDDO         ENDDO
456        ENDDO        ENDDO
# Line 462  C--   Upper diagonal Line 461  C--   Upper diagonal
461         ENDDO         ENDDO
462        ENDDO        ENDDO
463        DO k=2,Nr        DO k=2,Nr
 C- jmc: concerned that c(k) from k=kLow to k=Nr should always be zero  
464         DO j=jMin,jMax         DO j=jMin,jMax
465          DO i=iMin,iMax          DO i=iMin,iMax
466  C- jmc: this is dangerous since klowC could be zero if land column            kp1=MAX(1,MIN(klowC(i,j,bi,bj),k+1))
467  C       and would be better without recip_hFacC  C-    We keep recip_hFacC in the diffusive flux calculation,
468            kp1=min(klowC(i,j,bi,bj),k+1)  C-    but no hFacC in TKE volume control
469    C-    No need for maskC(k) with recip_hFacC(k)
470            c(i,j,k) = -deltaTggl90            c(i,j,k) = -deltaTggl90
 c     &        *recip_drF( k )*recip_hFacI(i,j,k,bi,bj)  
471       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)
472       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))
473       &        *recip_drC(k)*maskC(i,j,k,bi,bj)*maskC(i,j,k-1,bi,bj)       &        *recip_drC(k)*maskC(i,j,k-1,bi,bj)
474          ENDDO          ENDDO
475         ENDDO         ENDDO
476        ENDDO        ENDDO
477  C--   Center diagonal  C--   Center diagonal
478        DO k=1,Nr        DO k=1,Nr
479           km1 = MAX(k-1,1)
480         DO j=jMin,jMax         DO j=jMin,jMax
481          DO i=iMin,iMax          DO i=iMin,iMax
482            b(i,j,k) = 1. _d 0 - c(i,j,k) - a(i,j,k)            b(i,j,k) = 1. _d 0 - c(i,j,k) - a(i,j,k)
483       &        + ab15*deltaTggl90*GGL90ceps*SQRTTKE(I,J,K)       &        + ab15*deltaTggl90*GGL90ceps*SQRTTKE(I,J,K)
484       &        *rMixingLength(I,J,K)*maskC(i,j,k,bi,bj)       &        * rMixingLength(I,J,K)
485         &        * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)
486           ENDDO           ENDDO
487         ENDDO         ENDDO
488        ENDDO        ENDDO
489  C     end set up matrix  C     end set up matrix
490    
491  C     Apply boundary condition  C     Apply boundary condition
492  C- jmc: concerned about conservation when a or c are changed after computing b        kp1 = MIN(Nr,kSurf+1)
493        DO J=jMin,jMax        DO J=jMin,jMax
494         DO I=iMin,iMax         DO I=iMin,iMax
495  C     estimate friction velocity uStar from surface forcing  C     estimate friction velocity uStar from surface forcing
# Line 500  C     estimate friction velocity uStar f Line 500  C     estimate friction velocity uStar f
500       &              + surfaceForcingV(I,  J+1,bi,bj) ) )**2       &              + surfaceForcingV(I,  J+1,bi,bj) ) )**2
501       &                     )       &                     )
502  C     Dirichlet surface boundary condition for TKE  C     Dirichlet surface boundary condition for TKE
 C- jmc: would be much better to update the provisional TKE (i.e. gTKE) at k=2  
503          gTKE(I,J,kSurf) = MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)          gTKE(I,J,kSurf) = MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)
504       &                     *maskC(I,J,kSurf,bi,bj)       &                     *maskC(I,J,kSurf,bi,bj)
505            gTKE(i,j,kp1) = gTKE(i,j,kp1)
506         &                - a(i,j,kp1)*gTKE(i,j,kSurf)
507            a(i,j,kp1) = 0. _d 0
508  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom
509          kBottom   = MAX(kLowC(I,J,bi,bj),1)          kBottom   = MAX(kLowC(I,J,bi,bj),1)
510          gTKE(I,J,kBottom) = gTKE(I,J,kBottom)          gTKE(I,J,kBottom) = gTKE(I,J,kBottom)
# Line 512  C     Dirichlet bottom boundary conditio Line 514  C     Dirichlet bottom boundary conditio
514        ENDDO        ENDDO
515    
516  C     solve tri-diagonal system, and store solution on gTKE (previously rhs)  C     solve tri-diagonal system, and store solution on gTKE (previously rhs)
517  c     CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
518  c    I                        a, b, c,       I                        a, b, c,
519  c    U                        gTKE,       U                        gTKE,
520  c    O                        errCode,       O                        errCode,
521  c    I                        bi, bj, myThid )       I                        bi, bj, myThid )
522        CALL GGL90_SOLVE( bi, bj, iMin, iMax, jMin, jMax,  c     CALL GGL90_SOLVE( bi, bj, iMin, iMax, jMin, jMax,
523       I     a, b, c,  c    I     a, b, c,
524       U     gTKE,  c    U     gTKE,
525       I     myThid )  c    I     myThid )
526    
527  C     now update TKE  C     now update TKE
528        DO K=1,Nr        DO K=1,Nr
Legend:
Removed from v.1.14  
changed lines
  Added in v.1.15
  ViewVC Help
Powered by ViewVC 1.1.22