/[MITgcm]/MITgcm/pkg/mom_vecinv/mom_vi_hdissip.F
ViewVC logotype

Diff of /MITgcm/pkg/mom_vecinv/mom_vi_hdissip.F

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

revision 1.11 by jmc, Tue Sep 21 12:57:50 2004 UTC revision 1.12 by jmc, Tue Nov 2 01:05:02 2004 UTC
# Line 41  C     == Local variables == Line 41  C     == Local variables ==
41        _RL Zip,Zij,Zpj,Dim,Dij,Dmj,uD2,vD2,uD4,vD4        _RL Zip,Zij,Zpj,Dim,Dij,Dmj,uD2,vD2,uD4,vD4
42        _RL Alin,Alth,grdVrt,vg2,vg4        _RL Alin,Alth,grdVrt,vg2,vg4
43        LOGICAL useVariableViscosity        LOGICAL useVariableViscosity
44        integer klev        INTEGER klev
45    
46        useVariableViscosity=        useVariableViscosity=
47       &      (viscAhGrid*deltaTmom.NE.0.)       &      (viscAhGrid*deltaTmom.NE.0.)
# Line 76  C but this approximation will work on cu Line 76  C but this approximation will work on cu
76           grdVrt=max(grdVrt,           grdVrt=max(grdVrt,
77       &        abs((vort3(i+1,j+1)-vort3(i+1,j))*recip_DYG(i+1,j,bi,bj)))       &        abs((vort3(i+1,j+1)-vort3(i+1,j))*recip_DYG(i+1,j,bi,bj)))
78           Alth=viscC2leith*grdVrt*(rA(i,j,bi,bj)**1.5)           Alth=viscC2leith*grdVrt*(rA(i,j,bi,bj)**1.5)
79           Alin=viscAh+vg2*rA ( i , j ,bi,bj)           Alin=viscAhD+vg2*rA ( i , j ,bi,bj)
80           viscAh_D(i,j)=min(viscAhMax,Alin+Alth)           viscAh_D(i,j)=min(viscAhMax,Alin+Alth)
81    
82           Alth=viscC4leith*grdVrt*0.125*(rA(i,j,bi,bj)**2.5)           Alth=viscC4leith*grdVrt*0.125*(rA(i,j,bi,bj)**2.5)
83           Alin=viscA4+vg4*(rA ( i , j ,bi,bj)**2)           Alin=viscA4D+vg4*(rA ( i , j ,bi,bj)**2)
84           viscA4_D(i,j)=min(viscA4Max,Alin+Alth)           viscA4_D(i,j)=min(viscA4Max,Alin+Alth)
85    
86  C This is the vector magnitude of the vorticity gradient  C This is the vector magnitude of the vorticity gradient
# Line 98  C but this approximation will work on cu Line 98  C but this approximation will work on cu
98           grdVrt=max(grdVrt,           grdVrt=max(grdVrt,
99       &          abs((vort3(i,j-1)-vort3(i,j))*recip_DYG(i,j-1,bi,bj)))       &          abs((vort3(i,j-1)-vort3(i,j))*recip_DYG(i,j-1,bi,bj)))
100           Alth=viscC2leith*grdVrt*(rAz(i,j,bi,bj)**1.5)           Alth=viscC2leith*grdVrt*(rAz(i,j,bi,bj)**1.5)
101           Alin=viscAh+vg2*rAz( i , j ,bi,bj)           Alin=viscAhZ+vg2*rAz( i , j ,bi,bj)
102           viscAh_Z(i,j)=min(viscAhMax,Alin+Alth)           viscAh_Z(i,j)=min(viscAhMax,Alin+Alth)
103    
104           Alth=viscC4leith*grdVrt*0.125*(rAz(i,j,bi,bj)**2.5)           Alth=viscC4leith*grdVrt*0.125*(rAz(i,j,bi,bj)**2.5)
105           Alin=viscA4+vg4*(rAz( i , j ,bi,bj)**2)           Alin=viscA4Z+vg4*(rAz( i , j ,bi,bj)**2)
106           viscA4_Z(i,j)=min(viscA4Max,Alin+Alth)           viscA4_Z(i,j)=min(viscA4Max,Alin+Alth)
107          ENDDO          ENDDO
108         ENDDO         ENDDO
109        ELSE        ELSE
110         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
111          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
112           viscAh_D(i,j)=viscAh           viscAh_D(i,j)=viscAhD
113           viscAh_Z(i,j)=viscAh           viscAh_Z(i,j)=viscAhZ
114           viscA4_D(i,j)=viscA4           viscA4_D(i,j)=viscA4D
115           viscA4_Z(i,j)=viscA4           viscA4_Z(i,j)=viscA4Z
116          ENDDO          ENDDO
117         ENDDO         ENDDO
118        ENDIF        ENDIF
119    
120  C     - Laplacian  and bi-harmonic terms  C     - Laplacian  terms
121        DO j=2-Oly,sNy+Oly-1        IF ( viscC2leith.NE.0. .OR. viscAhGrid.NE.0.
122         DO i=2-Olx,sNx+Olx-1       &    .OR. viscAhD.NE.0. .OR. viscAhZ.NE.0. ) THEN
123           DO j=2-Oly,sNy+Oly-1
124          Dim=hDiv( i ,j-1)          DO i=2-Olx,sNx+Olx-1
125          Dij=hDiv( i , j )  
126          Dmj=hDiv(i-1, j )           Dim=hDiv( i ,j-1)
127          Zip=hFacZ( i ,j+1)*vort3( i ,j+1)           Dij=hDiv( i , j )
128          Zij=hFacZ( i , j )*vort3( i , j )           Dmj=hDiv(i-1, j )
129          Zpj=hFacZ(i+1, j )*vort3(i+1, j )           Zip=hFacZ( i ,j+1)*vort3( i ,j+1)
130             Zij=hFacZ( i , j )*vort3( i , j )
131             Zpj=hFacZ(i+1, j )*vort3(i+1, j )
132    
133  C This bit scales the harmonic dissipation operator to be proportional  C This bit scales the harmonic dissipation operator to be proportional
134  C to the grid-cell area over the time-step. viscAh is then non-dimensional  C to the grid-cell area over the time-step. viscAh is then non-dimensional
135  C and should be less than 1/8, for example viscAh=0.01  C and should be less than 1/8, for example viscAh=0.01
136          if (useVariableViscosity) then           IF (useVariableViscosity) THEN
137            Dij=Dij*viscAh_D(i,j)            Dij=Dij*viscAh_D(i,j)
138            Dim=Dim*viscAh_D(i,j-1)            Dim=Dim*viscAh_D(i,j-1)
139            Dmj=Dmj*viscAh_D(i-1,j)            Dmj=Dmj*viscAh_D(i-1,j)
# Line 145  C and should be less than 1/8, for examp Line 147  C and should be less than 1/8, for examp
147       &       recip_hFacS(i,j,k,bi,bj)*( Zpj-Zij )*recip_DXG(i,j,bi,bj)       &       recip_hFacS(i,j,k,bi,bj)*( Zpj-Zij )*recip_DXG(i,j,bi,bj)
148       &                                           *cosFacV(j,bi,bj)       &                                           *cosFacV(j,bi,bj)
149       &                               +( Dij-Dim )*recip_DYC(i,j,bi,bj) )       &                               +( Dij-Dim )*recip_DYC(i,j,bi,bj) )
150          else           ELSE
151            uD2 = viscAh*(            uD2 = viscAhD*
152       &               cosFacU(j,bi,bj)*( Dij-Dmj )*recip_DXC(i,j,bi,bj)       &               cosFacU(j,bi,bj)*( Dij-Dmj )*recip_DXC(i,j,bi,bj)
153       &      -recip_hFacW(i,j,k,bi,bj)*( Zip-Zij )*recip_DYG(i,j,bi,bj) )       &        - viscAhZ*recip_hFacW(i,j,k,bi,bj)*
154            vD2 = viscAh*(       &                                ( Zip-Zij )*recip_DYG(i,j,bi,bj)
155       &       recip_hFacS(i,j,k,bi,bj)*( Zpj-Zij )*recip_DXG(i,j,bi,bj)            vD2 = viscAhZ*recip_hFacS(i,j,k,bi,bj)*
156       &                                           *cosFacV(j,bi,bj)       &               cosFacV(j,bi,bj)*( Zpj-Zij )*recip_DXG(i,j,bi,bj)
157       &                               +( Dij-Dim )*recip_DYC(i,j,bi,bj) )       &        + viscAhD*              ( Dij-Dim )*recip_DYC(i,j,bi,bj)
158          endif           ENDIF
159    
160             uDissip(i,j) = uD2
161             vDissip(i,j) = vD2
162    
163            ENDDO
164           ENDDO
165          ELSE
166           DO j=2-Oly,sNy+Oly-1
167            DO i=2-Olx,sNx+Olx-1
168             uDissip(i,j) = 0.
169             vDissip(i,j) = 0.
170            ENDDO
171           ENDDO
172          ENDIF
173    
174    C     - Bi-harmonic terms
175          IF ( viscC4leith.NE.0. .OR. viscA4Grid.NE.0.
176         &    .OR. viscA4D.NE.0. .OR. viscA4Z.NE.0. ) THEN
177           DO j=2-Oly,sNy+Oly-1
178            DO i=2-Olx,sNx+Olx-1
179    
180  #ifdef MOM_VI_ORIGINAL_VISCA4  #ifdef MOM_VI_ORIGINAL_VISCA4
181          Dim=dyF( i ,j-1,bi,bj)*dStar( i ,j-1)           Dim=dyF( i ,j-1,bi,bj)*dStar( i ,j-1)
182          Dij=dyF( i , j ,bi,bj)*dStar( i , j )           Dij=dyF( i , j ,bi,bj)*dStar( i , j )
183          Dmj=dyF(i-1, j ,bi,bj)*dStar(i-1, j )           Dmj=dyF(i-1, j ,bi,bj)*dStar(i-1, j )
184    
185          Zip=dxV( i ,j+1,bi,bj)*hFacZ( i ,j+1)*zStar( i ,j+1)           Zip=dxV( i ,j+1,bi,bj)*hFacZ( i ,j+1)*zStar( i ,j+1)
186          Zij=dxV( i , j ,bi,bj)*hFacZ( i , j )*zStar( i , j )           Zij=dxV( i , j ,bi,bj)*hFacZ( i , j )*zStar( i , j )
187          Zpj=dxV(i+1, j ,bi,bj)*hFacZ(i+1, j )*zStar(i+1, j )           Zpj=dxV(i+1, j ,bi,bj)*hFacZ(i+1, j )*zStar(i+1, j )
188  #else  #else
189          Dim=dStar( i ,j-1)           Dim=dStar( i ,j-1)
190          Dij=dStar( i , j )           Dij=dStar( i , j )
191          Dmj=dStar(i-1, j )           Dmj=dStar(i-1, j )
192    
193          Zip=hFacZ( i ,j+1)*zStar( i ,j+1)           Zip=hFacZ( i ,j+1)*zStar( i ,j+1)
194          Zij=hFacZ( i , j )*zStar( i , j )           Zij=hFacZ( i , j )*zStar( i , j )
195          Zpj=hFacZ(i+1, j )*zStar(i+1, j )           Zpj=hFacZ(i+1, j )*zStar(i+1, j )
196  #endif  #endif
197    
   
198  C This bit scales the harmonic dissipation operator to be proportional  C This bit scales the harmonic dissipation operator to be proportional
199  C to the grid-cell area over the time-step. viscAh is then non-dimensional  C to the grid-cell area over the time-step. viscAh is then non-dimensional
200  C and should be less than 1/8, for example viscAh=0.01  C and should be less than 1/8, for example viscAh=0.01
201          IF (useVariableViscosity) THEN           IF (useVariableViscosity) THEN
202            Dij=Dij*viscA4_D(i,j)            Dij=Dij*viscA4_D(i,j)
203            Dim=Dim*viscA4_D(i,j-1)            Dim=Dim*viscA4_D(i,j-1)
204            Dmj=Dmj*viscA4_D(i-1,j)            Dmj=Dmj*viscA4_D(i-1,j)
# Line 192  C and should be less than 1/8, for examp Line 213  C and should be less than 1/8, for examp
213            vD4 = recip_rAs(i,j,bi,bj)*(            vD4 = recip_rAs(i,j,bi,bj)*(
214       &    recip_hFacS(i,j,k,bi,bj)*( (Zpj-Zij)*cosFacV(j,bi,bj) )       &    recip_hFacS(i,j,k,bi,bj)*( (Zpj-Zij)*cosFacV(j,bi,bj) )
215       &   +                         ( Dij-Dim ) )       &   +                         ( Dij-Dim ) )
216          ELSE           ELSE
217            uD4 = recip_rAw(i,j,bi,bj)*(            uD4 = recip_rAw(i,j,bi,bj)*(
218       &                             viscA4*( (Dij-Dmj)*cosFacU(j,bi,bj) )       &                             viscA4*( (Dij-Dmj)*cosFacU(j,bi,bj) )
219       &   -recip_hFacW(i,j,k,bi,bj)*viscA4*( Zip-Zij ) )       &   -recip_hFacW(i,j,k,bi,bj)*viscA4*( Zip-Zij ) )
# Line 207  C and should be less than 1/8, for examp Line 228  C and should be less than 1/8, for examp
228       &       recip_hFacS(i,j,k,bi,bj)*( Zpj-Zij )*recip_DXG(i,j,bi,bj)       &       recip_hFacS(i,j,k,bi,bj)*( Zpj-Zij )*recip_DXG(i,j,bi,bj)
229       &                                           *cosFacV(j,bi,bj)       &                                           *cosFacV(j,bi,bj)
230       &                               +( Dij-Dim )*recip_DYC(i,j,bi,bj) )       &                               +( Dij-Dim )*recip_DYC(i,j,bi,bj) )
231          ELSE           ELSE
232            uD4 = viscA4*(            uD4 = viscA4D*
233       &               cosFacU(j,bi,bj)*( Dij-Dmj )*recip_DXC(i,j,bi,bj)       &               cosFacU(j,bi,bj)*( Dij-Dmj )*recip_DXC(i,j,bi,bj)
234       &      -recip_hFacW(i,j,k,bi,bj)*( Zip-Zij )*recip_DYG(i,j,bi,bj) )       &        - viscA4Z*recip_hFacW(i,j,k,bi,bj)*
235            vD4 = viscA4*(       &                                ( Zip-Zij )*recip_DYG(i,j,bi,bj)
236       &       recip_hFacS(i,j,k,bi,bj)*( Zpj-Zij )*recip_DXG(i,j,bi,bj)            vD4 = viscA4Z*recip_hFacS(i,j,k,bi,bj)*
237       &                                           *cosFacV(j,bi,bj)       &               cosFacV(j,bi,bj)*( Zpj-Zij )*recip_DXG(i,j,bi,bj)
238       &                               +( Dij-Dim )*recip_DYC(i,j,bi,bj) )       &        + viscA4D*              ( Dij-Dim )*recip_DYC(i,j,bi,bj)
239  #endif /* MOM_VI_ORIGINAL_VISCA4 */  #endif /* MOM_VI_ORIGINAL_VISCA4 */
240          ENDIF           ENDIF
241    
242          uDissip(i,j) = uD2 - uD4           uDissip(i,j) = uDissip(i,j) - uD4
243          vDissip(i,j) = vD2 - vD4           vDissip(i,j) = vDissip(i,j) - vD4
244    
245            ENDDO
246         ENDDO         ENDDO
247        ENDDO        ENDIF
248    
249  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
250        if (useDiagnostics) then        if (useDiagnostics) then
Legend:
Removed from v.1.11  
changed lines
  Added in v.1.12
  ViewVC Help
Powered by ViewVC 1.1.22