/[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.2.2.1 by adcroft, Tue Feb 26 16:04:49 2002 UTC revision 1.36 by jmc, Thu Sep 10 18:08:51 2015 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "MOM_VECINV_OPTIONS.h"
5    
6        SUBROUTINE MOM_VI_HDISSIP(        SUBROUTINE MOM_VI_HDISSIP(
7       I        bi,bj,k,       I        bi, bj, k,
8       I        hDiv,vort3,hFacZ,dStar,zStar,       I        hDiv, vort3, dStar, zStar, hFacZ,
9       O        uDissip,vDissip,       I        viscAh_Z, viscAh_D, viscA4_Z, viscA4_D,
10       I        myThid)       I        harmonic, biharmonic, useVariableViscosity,
11         O        uDissip, vDissip,
12         I        myThid )
13    
14        IMPLICIT NONE        IMPLICIT NONE
15  C  
16  C     Calculate horizontal dissipation terms  C     Calculate horizontal dissipation terms
17  C     [del^2 - del^4] (u,v)  C     [del^2 - del^4] (u,v)
 C  
18    
19  C     == Global variables ==  C     == Global variables ==
20  #include "SIZE.h"  #include "SIZE.h"
 #include "GRID.h"  
21  #include "EEPARAMS.h"  #include "EEPARAMS.h"
22  #include "PARAMS.h"  #include "PARAMS.h"
23    #include "GRID.h"
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)        _RL hDiv (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
28        _RL vort3(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vort3(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RS hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
29        _RL dStar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL dStar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
30        _RL zStar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL zStar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
31          _RS hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
32          _RL viscAh_Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
33          _RL viscAh_D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
34          _RL viscA4_Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
35          _RL viscA4_D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
36          LOGICAL harmonic, biharmonic, useVariableViscosity
37        _RL uDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
38        _RL vDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
39        INTEGER myThid        INTEGER myThid
40    
41  C     == Local variables ==  C     == Local variables ==
42        INTEGER I,J        INTEGER i, j
43        _RL Zip,Zij,Zpj,Dim,Dij,Dmj,uD2,vD2,uD4,vD4        _RL Zip, Zij, Zpj, Dim, Dij, Dmj, uD2, vD2, uD4, vD4
44          _RL Zip1, Zij1, Zpj1
45    
46    C     - Laplacian  terms
47          IF (harmonic) THEN
48    C This bit scales the harmonic dissipation operator to be proportional
49    C to the grid-cell area over the time-step. viscAh is then non-dimensional
50    C and should be less than 1/8, for example viscAh=0.01
51           IF (useVariableViscosity) THEN
52            DO j=2-OLy,sNy+OLy-1
53             DO i=2-OLx,sNx+OLx-1
54    
55              Dij=hDiv( i , j )*viscAh_D(i,j)
56              Dim=hDiv( i ,j-1)*viscAh_D(i,j-1)
57              Dmj=hDiv(i-1, j )*viscAh_D(i-1,j)
58              Zij=hFacZ( i , j )*vort3( i , j )*viscAh_Z(i,j)
59              Zip=hFacZ( i ,j+1)*vort3( i ,j+1)*viscAh_Z(i,j+1)
60              Zpj=hFacZ(i+1, j )*vort3(i+1, j )*viscAh_Z(i+1,j)
61    
62  C     - Laplacian  and bi-harmonic terms            uD2 = (
       DO j=2-Oly,sNy+Oly-1  
        DO i=2-Olx,sNx+Olx-1  
   
 c       Dim=dyF( i ,j-1,bi,bj)*hFacC( i ,j-1,k,bi,bj)*hDiv( i ,j-1)  
 c       Dij=dyF( i , j ,bi,bj)*hFacC( i , j ,k,bi,bj)*hDiv( i , j )  
 c       Dmj=dyF(i-1, j ,bi,bj)*hFacC(i-1, j ,k,bi,bj)*hDiv(i-1, j )  
 c       Dim=dyF( i ,j-1,bi,bj)*                       hDiv( i ,j-1)  
 c       Dij=dyF( i , j ,bi,bj)*                       hDiv( i , j )  
 c       Dmj=dyF(i-1, j ,bi,bj)*                       hDiv(i-1, j )  
         Dim=                                          hDiv( i ,j-1)  
         Dij=                                          hDiv( i , j )  
         Dmj=                                          hDiv(i-1, j )  
   
 c       Zip=dxV( i ,j+1,bi,bj)*hFacZ( i ,j+1)*vort3( i ,j+1)  
 c       Zij=dxV( i , j ,bi,bj)*hFacZ( i , j )*vort3( i , j )  
 c       Zpj=dxV(i+1, j ,bi,bj)*hFacZ(i+1, j )*vort3(i+1, j )  
         Zip=                   hFacZ( i ,j+1)*vort3( i ,j+1)  
         Zij=                   hFacZ( i , j )*vort3( i , j )  
         Zpj=                   hFacZ(i+1, j )*vort3(i+1, j )  
   
 c       uD2 = recip_rAw(i,j,bi,bj)*(  
 c    &    recip_hFacW(i,j,k,bi,bj)*viscAh*( (Dij-Dmj)*cosFacU(j,bi,bj) )  
 c    &   -recip_hFacW(i,j,k,bi,bj)*viscAh*( Zip-Zij ) )  
 c       uD2 = recip_rAw(i,j,bi,bj)*(  
 c    &                             viscAh*( (Dij-Dmj)*cosFacU(j,bi,bj) )  
 c    &   -recip_hFacW(i,j,k,bi,bj)*viscAh*( Zip-Zij ) )  
         uD2 = viscAh*(  
63       &               cosFacU(j,bi,bj)*( Dij-Dmj )*recip_DXC(i,j,bi,bj)       &               cosFacU(j,bi,bj)*( Dij-Dmj )*recip_DXC(i,j,bi,bj)
64       &      -recip_hFacW(i,j,k,bi,bj)*( Zip-Zij )*recip_DYG(i,j,bi,bj) )       &     -_recip_hFacW(i,j,k,bi,bj)*( Zip-Zij )*recip_DYG(i,j,bi,bj) )
65    #ifdef ISOTROPIC_COS_SCALING
66         &                                           *cosFacU(j,bi,bj)
67    #endif /* ISOTROPIC_COS_SCALING */
68              vD2 = (
69         &      _recip_hFacS(i,j,k,bi,bj)*( Zpj-Zij )*recip_DXG(i,j,bi,bj)
70         &                                           *cosFacV(j,bi,bj)
71         &                               +( Dij-Dim )*recip_DYC(i,j,bi,bj) )
72    #ifdef ISOTROPIC_COS_SCALING
73         &                                           *cosFacV(j,bi,bj)
74    #endif /* ISOTROPIC_COS_SCALING */
75    
76  c       vD2 = recip_rAs(i,j,bi,bj)*(            uDissip(i,j) = uD2
77  c    &    recip_hFacS(i,j,k,bi,bj)*viscAh*( (Zpj-Zij)*cosFacV(j,bi,bj) )            vDissip(i,j) = vD2
78  c    &   +recip_hFacS(i,j,k,bi,bj)*viscAh*( Dij-Dim ) )  
79  c       vD2 = recip_rAs(i,j,bi,bj)*(           ENDDO
80  c    &    recip_hFacS(i,j,k,bi,bj)*viscAh*( (Zpj-Zij)*cosFacV(j,bi,bj) )          ENDDO
81  c    &   +                         viscAh*( Dij-Dim ) )         ELSE
82          vD2 = viscAh*(          DO j=2-OLy,sNy+OLy-1
83       &       recip_hFacS(i,j,k,bi,bj)*( Zpj-Zij )*recip_DXG(i,j,bi,bj)           DO i=2-OLx,sNx+OLx-1
84    
85              Dim=hDiv( i ,j-1)
86              Dij=hDiv( i , j )
87              Dmj=hDiv(i-1, j )
88              Zip=hFacZ( i ,j+1)*vort3( i ,j+1)
89              Zij=hFacZ( i , j )*vort3( i , j )
90              Zpj=hFacZ(i+1, j )*vort3(i+1, j )
91    
92              uD2 = viscAhD*
93         &               cosFacU(j,bi,bj)*( Dij-Dmj )*recip_DXC(i,j,bi,bj)
94         &        - viscAhZ*_recip_hFacW(i,j,k,bi,bj)*
95         &                                ( Zip-Zij )*recip_DYG(i,j,bi,bj)
96    #ifdef ISOTROPIC_COS_SCALING
97         &                                           *cosFacU(j,bi,bj)
98    #endif /* ISOTROPIC_COS_SCALING */
99              vD2 = viscAhZ*_recip_hFacS(i,j,k,bi,bj)*
100         &               cosFacV(j,bi,bj)*( Zpj-Zij )*recip_DXG(i,j,bi,bj)
101         &        + viscAhD*              ( Dij-Dim )*recip_DYC(i,j,bi,bj)
102    #ifdef ISOTROPIC_COS_SCALING
103         &                                           *cosFacV(j,bi,bj)
104    #endif /* ISOTROPIC_COS_SCALING */
105    
106              uDissip(i,j) = uD2
107              vDissip(i,j) = vD2
108    
109             ENDDO
110            ENDDO
111           ENDIF
112          ELSE
113           DO j=2-OLy,sNy+OLy-1
114            DO i=2-OLx,sNx+OLx-1
115             uDissip(i,j) = 0.
116             vDissip(i,j) = 0.
117            ENDDO
118           ENDDO
119          ENDIF
120    
121    C     - Bi-harmonic terms
122          IF (biharmonic) THEN
123    
124    C This bit scales the harmonic dissipation operator to be proportional
125    C to the grid-cell area over the time-step. viscAh is then non-dimensional
126    C and should be less than 1/8, for example viscAh=0.01
127           IF (useVariableViscosity) THEN
128            DO j=2-OLy,sNy+OLy-1
129             DO i=2-OLx,sNx+OLx-1
130    
131    #ifdef MOM_VI_ORIGINAL_VISCA4
132              Dim=dyF( i ,j-1,bi,bj)*dStar( i ,j-1)
133              Dij=dyF( i , j ,bi,bj)*dStar( i , j )
134              Dmj=dyF(i-1, j ,bi,bj)*dStar(i-1, j )
135    
136              Zip1=dxV( i ,j+1,bi,bj)*hFacZ( i ,j+1)*zStar( i ,j+1)
137              Zij1=dxV( i , j ,bi,bj)*hFacZ( i , j )*zStar( i , j )
138              Zpj1=dxV(i+1, j ,bi,bj)*hFacZ(i+1, j )*zStar(i+1, j )
139    #else
140              Dim=dStar( i ,j-1)
141              Dij=dStar( i , j )
142              Dmj=dStar(i-1, j )
143    
144              Zip1=hFacZ( i ,j+1)*zStar( i ,j+1)
145              Zij1=hFacZ( i , j )*zStar( i , j )
146              Zpj1=hFacZ(i+1, j )*zStar(i+1, j )
147    #endif
148              Dij=Dij*viscA4_D(i,j)
149              Dim=Dim*viscA4_D(i,j-1)
150              Dmj=Dmj*viscA4_D(i-1,j)
151              Zij=Zij1*viscA4_Z(i,j)
152              Zip=Zip1*viscA4_Z(i,j+1)
153              Zpj=Zpj1*viscA4_Z(i+1,j)
154    
155    #ifdef MOM_VI_ORIGINAL_VISCA4
156              uD4 = recip_rAw(i,j,bi,bj)*(
157         &                             ( (Dij-Dmj)*cosFacU(j,bi,bj) )
158         &  -_recip_hFacW(i,j,k,bi,bj)*( Zip-Zij )
159    # ifdef ISOTROPIC_COS_SCALING
160         &                                        *cosFacU(j,bi,bj)
161    # endif /* ISOTROPIC_COS_SCALING */
162         &         )
163              vD4 = recip_rAs(i,j,bi,bj)*(
164         &   _recip_hFacS(i,j,k,bi,bj)*( (Zpj-Zij)*cosFacV(j,bi,bj) )
165         &   +                         ( Dij-Dim )
166    # ifdef ISOTROPIC_COS_SCALING
167         &                                        *cosFacV(j,bi,bj)
168    # endif /* ISOTROPIC_COS_SCALING */
169         &         )
170    #else /* MOM_VI_ORIGINAL_VISCA4 */
171              uD4 = (
172         &               cosFacU(j,bi,bj)*( Dij-Dmj )*recip_DXC(i,j,bi,bj)
173         &     -_recip_hFacW(i,j,k,bi,bj)*( Zip-Zij )*recip_DYG(i,j,bi,bj) )
174    # ifdef ISOTROPIC_COS_SCALING
175         &                                           *cosFacU(j,bi,bj)
176    # endif /* ISOTROPIC_COS_SCALING */
177              vD4 = (
178         &      _recip_hFacS(i,j,k,bi,bj)*( Zpj-Zij )*recip_DXG(i,j,bi,bj)
179       &                                           *cosFacV(j,bi,bj)       &                                           *cosFacV(j,bi,bj)
180       &                               +( Dij-Dim )*recip_DYC(i,j,bi,bj) )       &                               +( Dij-Dim )*recip_DYC(i,j,bi,bj) )
181    # ifdef ISOTROPIC_COS_SCALING
182         &                                           *cosFacV(j,bi,bj)
183    # endif /* ISOTROPIC_COS_SCALING */
184    #endif  /* MOM_VI_ORIGINAL_VISCA4 */
185    
186              uDissip(i,j) = uDissip(i,j) - uD4
187              vDissip(i,j) = vDissip(i,j) - vD4
188    
189  c       Dim=dyF( i ,j-1,bi,bj)*hFacC( i ,j-1,k,bi,bj)*dStar( i ,j-1)           ENDDO
190  c       Dij=dyF( i , j ,bi,bj)*hFacC( i , j ,k,bi,bj)*dStar( i , j )          ENDDO
191  c       Dmj=dyF(i-1, j ,bi,bj)*hFacC(i-1, j ,k,bi,bj)*dStar(i-1, j )         ELSE
192          Dim=dyF( i ,j-1,bi,bj)*                       dStar( i ,j-1)          DO j=2-OLy,sNy+OLy-1
193          Dij=dyF( i , j ,bi,bj)*                       dStar( i , j )           DO i=2-OLx,sNx+OLx-1
194          Dmj=dyF(i-1, j ,bi,bj)*                       dStar(i-1, j )  
195    #ifdef MOM_VI_ORIGINAL_VISCA4
196          Zip=dxV( i ,j+1,bi,bj)*hFacZ( i ,j+1)*zStar( i ,j+1)            Dim=dyF( i ,j-1,bi,bj)*dStar( i ,j-1)
197          Zij=dxV( i , j ,bi,bj)*hFacZ( i , j )*zStar( i , j )            Dij=dyF( i , j ,bi,bj)*dStar( i , j )
198          Zpj=dxV(i+1, j ,bi,bj)*hFacZ(i+1, j )*zStar(i+1, j )            Dmj=dyF(i-1, j ,bi,bj)*dStar(i-1, j )
199    
200  c       uD4 = recip_rAw(i,j,bi,bj)*(            Zip1=dxV( i ,j+1,bi,bj)*hFacZ( i ,j+1)*zStar( i ,j+1)
201  c    &    recip_hFacW(i,j,k,bi,bj)*viscA4*( (Dij-Dmj)*cosFacU(j,bi,bj) )            Zij1=dxV( i , j ,bi,bj)*hFacZ( i , j )*zStar( i , j )
202  c    &   -recip_hFacW(i,j,k,bi,bj)*viscA4*( Zip-Zij ) )            Zpj1=dxV(i+1, j ,bi,bj)*hFacZ(i+1, j )*zStar(i+1, j )
203          uD4 = recip_rAw(i,j,bi,bj)*(  #else
204       &                             viscA4*( (Dij-Dmj)*cosFacU(j,bi,bj) )            Dim=dStar( i ,j-1)
205       &   -recip_hFacW(i,j,k,bi,bj)*viscA4*( Zip-Zij ) )            Dij=dStar( i , j )
206              Dmj=dStar(i-1, j )
207  c       vD4 = recip_rAs(i,j,bi,bj)*(  
208  c    &    recip_hFacS(i,j,k,bi,bj)*viscA4*( (Zpj-Zij)*cosFacV(j,bi,bj) )            Zip1=hFacZ( i ,j+1)*zStar( i ,j+1)
209  c    &   +recip_hFacS(i,j,k,bi,bj)*viscA4*( Dij-Dim ) )            Zij1=hFacZ( i , j )*zStar( i , j )
210          vD4 = recip_rAs(i,j,bi,bj)*(            Zpj1=hFacZ(i+1, j )*zStar(i+1, j )
211       &    recip_hFacS(i,j,k,bi,bj)*viscA4*( (Zpj-Zij)*cosFacV(j,bi,bj) )  #endif
212       &   +                         viscA4*( Dij-Dim ) )            Zij=Zij1
213              Zip=Zip1
214              Zpj=Zpj1
215    
216    #ifdef MOM_VI_ORIGINAL_VISCA4
217              uD4 = recip_rAw(i,j,bi,bj)*(
218         &                             viscA4D*( Dij-Dmj )*cosFacU(j,bi,bj)
219         &  -_recip_hFacW(i,j,k,bi,bj)*viscA4Z*( Zip-Zij )
220    # ifdef ISOTROPIC_COS_SCALING
221         &                                                *cosFacU(j,bi,bj)
222    # endif /* ISOTROPIC_COS_SCALING */
223         &                               )
224              vD4 = recip_rAs(i,j,bi,bj)*(
225         &   _recip_hFacS(i,j,k,bi,bj)*viscA4Z*( Zpj-Zij )*cosFacV(j,bi,bj)
226         &   +                         viscA4D*( Dij-Dim )
227    # ifdef ISOTROPIC_COS_SCALING
228         &                                                *cosFacV(j,bi,bj)
229    # endif /* ISOTROPIC_COS_SCALING */
230         &                               )
231    #else /* MOM_VI_ORIGINAL_VISCA4 */
232              uD4 = viscA4D*
233         &               cosFacU(j,bi,bj)*( Dij-Dmj )*recip_DXC(i,j,bi,bj)
234         &        - viscA4Z*_recip_hFacW(i,j,k,bi,bj)*
235         &                                ( Zip-Zij )*recip_DYG(i,j,bi,bj)
236    # ifdef ISOTROPIC_COS_SCALING
237         &                                           *cosFacU(j,bi,bj)
238    # endif /* ISOTROPIC_COS_SCALING */
239              vD4 = viscA4Z*_recip_hFacS(i,j,k,bi,bj)*
240         &               cosFacV(j,bi,bj)*( Zpj-Zij )*recip_DXG(i,j,bi,bj)
241         &        + viscA4D*              ( Dij-Dim )*recip_DYC(i,j,bi,bj)
242    # ifdef ISOTROPIC_COS_SCALING
243         &                                           *cosFacV(j,bi,bj)
244    # endif /* ISOTROPIC_COS_SCALING */
245    #endif /* MOM_VI_ORIGINAL_VISCA4 */
246    
247          uDissip(i,j) = uD2 - uD4            uDissip(i,j) = uDissip(i,j) - uD4
248          vDissip(i,j) = vD2 - vD4            vDissip(i,j) = vDissip(i,j) - vD4
249    
250             ENDDO
251            ENDDO
252           ENDIF
253          ENDIF
254    
255          IF ( harmonic .OR. biharmonic ) THEN
256           DO j=1-OLy,sNy+OLy-1
257            DO i=1-OLx,sNx+OLx-1
258             uDissip(i,j) = uDissip(i,j)*maskW(i,j,k,bi,bj)
259         &                              *recip_deepFacC(k)
260             vDissip(i,j) = vDissip(i,j)*maskS(i,j,k,bi,bj)
261         &                              *recip_deepFacC(k)
262            ENDDO
263         ENDDO         ENDDO
264        ENDDO        ENDIF
265    
266        RETURN        RETURN
267        END        END
Legend:
Removed from v.1.2.2.1  
changed lines
  Added in v.1.36
  ViewVC Help
Powered by ViewVC 1.1.22