/[MITgcm]/MITgcm/pkg/mom_common/mom_hdissip.F
ViewVC logotype

Diff of /MITgcm/pkg/mom_common/mom_hdissip.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.3 by baylor, Thu Apr 14 14:22:51 2005 UTC revision 1.4 by baylor, Fri Sep 16 19:32:20 2005 UTC
# Line 5  C $Name$ Line 5  C $Name$
5    
6        SUBROUTINE MOM_HDISSIP(        SUBROUTINE MOM_HDISSIP(
7       I        bi,bj,k,       I        bi,bj,k,
8       I        tension,strain,hFacZ,viscAt,viscAs,       I        hDiv,vort3,tension,strain,KE,
9         I        hFacZ,
10       O        uDissip,vDissip,       O        uDissip,vDissip,
11       I        myThid)       I        myThid)
12        IMPLICIT NONE        IMPLICIT NONE
# Line 23  C     == Global variables == Line 24  C     == Global variables ==
24    
25  C     == Routine arguments ==  C     == Routine arguments ==
26        INTEGER bi,bj,k        INTEGER bi,bj,k
27          _RL hDiv(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
28          _RL vort3(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
29        _RL tension(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL tension(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
30        _RL strain(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL strain(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
31          _RL KE(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
32        _RS hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL viscAt  
       _RL viscAs  
33        _RL uDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
34        _RL vDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
35        INTEGER myThid        INTEGER myThid
36    
37  C     == Local variables ==  C     == Local variables ==
38        INTEGER I,J        INTEGER I,J
39        _RL viscA_t(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscAh_t(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
40        _RL viscA_s(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscAh_s(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
41        _RL ASmag, smagfac        _RL viscA4_t(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
42        _RL vg2Min, vg2Max, AlinMax, AlinMin        _RL viscA4_s(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
43        _RL lenA2, lenAz2        LOGICAL harmonic, biharmonic, useVariableViscosity
   
       IF (deltaTmom.NE.0.) THEN  
        vg2Min=viscAhGridMin/deltaTmom  
        vg2Max=viscAhGridMax/deltaTmom  
       ELSE  
        vg2Min=0.  
        vg2Max=0.  
       ENDIF  
44    
45  C     - Calculate Smagorinsky Coefficients        CALL MOM_CALC_VISC(
46        smagfac=(viscC2smag/pi)**2       I        bi,bj,k,
47        DO j=2-Oly,sNy+Oly-1       O        viscAh_s,viscAh_t,viscA4_s,viscA4_t,
48         DO i=2-Olx,sNx+Olx-1       O        harmonic,biharmonic,useVariableViscosity,
49          IF (viscC2smag.NE.0.) THEN       I        hDiv,vort3,tension,strain,KE,
50  C    Geometric Mean is used as length scale       I        myThid)
          lenA2=(2*rA(i,j,bi,bj)/  
      &      (dxF(I,J,bi,bj)+dyF(I,J,bi,bj)))**2  
   
          Asmag=smagfac*lenA2  
      &    *sqrt(tension(i,j)**2  
      &          +0.25*(strain(i+1, j )**2+strain( i ,j+1)**2  
      &                +strain(i-1, j )**2+strain( i ,j-1)**2))  
          viscA_t(i,j)=min(viscAhMax,max(viscAt,Asmag))  
   
          IF (vg2Max.GT.0.) THEN  
             AlinMax=vg2Max*lenA2  
             viscA_t(i,j)=min(AlinMax,viscA_t(i,j))  
          ENDIF  
          AlinMin=vg2Min*lenA2  
          viscA_t(i,j)=max(AlinMin,viscA_t(i,j))  
   
 C    Geometric Mean is used as length scale  
          lenAz2=(2*rAz(i,j,bi,bj)/  
      &       (dxV(I,J,bi,bj)+dyU(I,J,bi,bj)))**2  
          Asmag=smagfac*lenAz2  
      &    *sqrt(strain(i,j)**2  
      &          +0.25*(tension( i , j )**2+tension( i ,j+1)**2  
      &                +tension(i+1, j )**2+tension(i+1,j+1)**2))  
          viscA_s(i,j)=min(viscAhMax,max(viscAs,Asmag))  
   
          IF (vg2Max.GT.0.) THEN  
             AlinMax=vg2Max*lenAz2  
             viscA_s(i,j)=min(AlinMax,viscA_s(i,j))  
          ENDIF  
          AlinMin=vg2Min*lenAz2  
          viscA_s(i,j)=max(AlinMin,viscA_s(i,j))  
   
         ELSE  
          viscA_t(i,j)=viscAt  
          viscA_s(i,j)=viscAs  
         ENDIF  
        ENDDO  
       ENDDO  
51    
52  C     - Laplacian  and bi-harmonic terms  C     - Laplacian and bi-harmonic terms
53        DO j=2-Oly,sNy+Oly-1        IF (harmonic) THEN
54         DO i=2-Olx,sNx+Olx-1         DO j=2-Oly,sNy+Oly-1
55            DO i=2-Olx,sNx+Olx-1
56          uDissip(i,j) =  
57       &   recip_dyg(i,j,bi,bj)*recip_dyg(i,j,bi,bj)           uDissip(i,j) =
58       &   *recip_dxc(i,j,bi,bj)       &    recip_dyg(i,j,bi,bj)*recip_dyg(i,j,bi,bj)
59       &   *(       &    *recip_dxc(i,j,bi,bj)
60       &      dyf( i , j ,bi,bj)*dyf( i , j ,bi,bj)       &    *(
61       &        *viscA_t( i , j )*tension( i , j )       &       dyf( i , j ,bi,bj)*dyf( i , j ,bi,bj)
62       &     -dyf(i-1, j ,bi,bj)*dyf(i-1, j ,bi,bj)       &         *viscAh_t( i , j )*tension( i , j )
63       &        *viscA_t(i-1, j )*tension(i-1, j )       &      -dyf(i-1, j ,bi,bj)*dyf(i-1, j ,bi,bj)
64       &    )       &         *viscAh_t(i-1, j )*tension(i-1, j )
65       &   +recip_dxc(i,j,bi,bj)*recip_dxc(i,j,bi,bj)       &     )
66       &   *recip_dyg(i,j,bi,bj)       &    +recip_dxc(i,j,bi,bj)*recip_dxc(i,j,bi,bj)
67       &   *(       &    *recip_dyg(i,j,bi,bj)
68       &      dxv( i ,j+1,bi,bj)*dxv( i ,j+1,bi,bj)       &    *(
69       &        *viscA_s(i,j+1)*strain( i ,j+1)       &       dxv( i ,j+1,bi,bj)*dxv( i ,j+1,bi,bj)
70       &     -dxv( i , j ,bi,bj)*dxv( i , j ,bi,bj)       &         *viscAh_s(i,j+1)*strain( i ,j+1)
71       &        *viscA_s(i, j )*strain( i , j )       &      -dxv( i , j ,bi,bj)*dxv( i , j ,bi,bj)
72       &    )       &         *viscAh_s(i, j )*strain( i , j )
73         &     )
74          vDissip(i,j) =  
75       &   recip_dyc(i,j,bi,bj)*recip_dyc(i,j,bi,bj)           vDissip(i,j) =
76       &   *recip_dxg(i,j,bi,bj)       &    recip_dyc(i,j,bi,bj)*recip_dyc(i,j,bi,bj)
77       &   *(       &    *recip_dxg(i,j,bi,bj)
78       &      dyu(i+1, j ,bi,bj)*dyu(i+1, j ,bi,bj)       &    *(
79       &        *viscA_s(i+1,j)*strain(i+1,j)       &       dyu(i+1, j ,bi,bj)*dyu(i+1, j ,bi,bj)
80       &     -dyu( i , j ,bi,bj)*dyu( i , j ,bi,bj)       &         *viscAh_s(i+1,j)*strain(i+1,j)
81       &        *viscA_s( i ,j)*strain( i ,j)       &      -dyu( i , j ,bi,bj)*dyu( i , j ,bi,bj)
82       &    )       &         *viscAh_s( i ,j)*strain( i ,j)
83       &   -recip_dxg(i,j,bi,bj)*recip_dxg(i,j,bi,bj)       &     )
84       &   *recip_dyc(i,j,bi,bj)       &    -recip_dxg(i,j,bi,bj)*recip_dxg(i,j,bi,bj)
85       &   *(       &    *recip_dyc(i,j,bi,bj)
86       &      dxf( i , j ,bi,bj)*dxf( i , j ,bi,bj)       &    *(
87       &        *viscA_t(i, j )*tension(i, j )       &       dxf( i , j ,bi,bj)*dxf( i , j ,bi,bj)
88       &     -dxf( i ,j-1,bi,bj)*dxf( i ,j-1,bi,bj)       &         *viscAh_t(i, j )*tension(i, j )
89       &        *viscA_t(i,j-1)*tension(i,j-1)       &      -dxf( i ,j-1,bi,bj)*dxf( i ,j-1,bi,bj)
90       &    )       &         *viscAh_t(i,j-1)*tension(i,j-1)
91         &     )
92    
93            ENDDO
94         ENDDO         ENDDO
95        ENDDO        ENDIF
96          IF (biharmonic) THEN
97           STOP 'MOM_HDISSIP: BIHARMONIC NOT ALLOWED WITH STRAIN-TENSION'
98          ENDIF
99    
100        RETURN        RETURN
101        END        END
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.22