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

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

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


Revision 1.18 - (hide annotations) (download)
Fri Mar 11 00:53:17 2005 UTC (19 years, 2 months ago) by jmc
Branch: MAIN
Changes since 1.17: +33 -30 lines
fix few rA,rAz confusion (from the previous check-in).

1 jmc 1.18 C $Header: /u/gcmpack/MITgcm/pkg/mom_vecinv/mom_vi_hdissip.F,v 1.17 2005/03/10 02:39:56 baylor Exp $
2 jmc 1.3 C $Name: $
3 adcroft 1.2
4 adcroft 1.8 #include "MOM_VECINV_OPTIONS.h"
5 adcroft 1.2
6     SUBROUTINE MOM_VI_HDISSIP(
7     I bi,bj,k,
8     I hDiv,vort3,hFacZ,dStar,zStar,
9     O uDissip,vDissip,
10     I myThid)
11 heimbach 1.14
12     cph(
13 jmc 1.15 cph The following line was commented in the argument list
14 heimbach 1.14 cph TAMC cannot digest commented lines within continuing lines
15     c I viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,
16     cph)
17    
18 adcroft 1.2 IMPLICIT NONE
19     C
20     C Calculate horizontal dissipation terms
21     C [del^2 - del^4] (u,v)
22     C
23    
24     C == Global variables ==
25     #include "SIZE.h"
26     #include "GRID.h"
27     #include "EEPARAMS.h"
28     #include "PARAMS.h"
29    
30     C == Routine arguments ==
31     INTEGER bi,bj,k
32     _RL hDiv(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
33     _RL vort3(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
34     _RS hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
35     _RL dStar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
36     _RL zStar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
37     _RL uDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
38     _RL vDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
39     INTEGER myThid
40    
41     C == Local variables ==
42 jmc 1.18 _RL viscAh_Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
43     _RL viscAh_D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
44     _RL viscA4_Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
45     _RL viscA4_D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
46 adcroft 1.2 INTEGER I,J
47     _RL Zip,Zij,Zpj,Dim,Dij,Dmj,uD2,vD2,uD4,vD4
48 jmc 1.18 _RL Alin,AlinMin,AlinMax,Alth2,Alth4,grdVrt,grdDiv
49 baylor 1.17 _RL vg2,vg2Min,vg2Max,vg4,vg4Min,vg4Max
50 adcroft 1.6 LOGICAL useVariableViscosity
51    
52     useVariableViscosity=
53 jmc 1.18 & (viscAhGrid.NE.0.)
54     & .OR.(viscA4Grid.NE.0.)
55 adcroft 1.6 & .OR.(viscC2leith.NE.0.)
56 baylor 1.17 & .OR.(viscC2leithD.NE.0.)
57 adcroft 1.6 & .OR.(viscC4leith.NE.0.)
58 baylor 1.17 & .OR.(viscC4leithD.NE.0.)
59 adcroft 1.6
60     IF (deltaTmom.NE.0.) THEN
61     vg2=viscAhGrid/deltaTmom
62 baylor 1.17 vg2Min=viscAhGridMin/deltaTmom
63     vg2Max=viscAhGridMax/deltaTmom
64 adcroft 1.6 vg4=viscA4Grid/deltaTmom
65 dimitri 1.13 vg4Min=viscA4GridMin/deltaTmom
66     vg4Max=viscA4GridMax/deltaTmom
67 adcroft 1.6 ELSE
68     vg2=0.
69 baylor 1.17 vg2Min=0.
70     vg2Max=0.
71 adcroft 1.6 vg4=0.
72 dimitri 1.13 vg4Min=0.
73     vg4Max=0.
74 adcroft 1.6 ENDIF
75    
76     C - Viscosity
77     IF (useVariableViscosity) THEN
78     DO j=2-Oly,sNy+Oly-1
79     DO i=2-Olx,sNx+Olx-1
80 jmc 1.18 IF (useFullLeith) THEN
81 baylor 1.17 C This is the vector magnitude of the vorticity gradient squared
82     grdVrt=0.25*(
83     & ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2
84     & +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2
85     & +((vort3(i+1,j+1)-vort3(i,j+1))*recip_DXG(i,j+1,bi,bj))**2
86     & +((vort3(i+1,j+1)-vort3(i+1,j))*recip_DYG(i+1,j,bi,bj))**2)
87    
88     C This is the vector magnitude of grad (div.v) squared
89     C Using it in Leith serves to damp instabilities in w.
90     grdDiv=0.25*(
91     & ((hDiv(i+1,j)-hDiv(i,j))*recip_DXG(i,j,bi,bj))**2
92     & +((hDiv(i,j+1)-hDiv(i,j))*recip_DYG(i,j,bi,bj))**2
93     & +((hDiv(i-1,j)-hDiv(i,j))*recip_DXG(i-1,j,bi,bj))**2
94     & +((hDiv(i,j-1)-hDiv(i,j))*recip_DYG(i,j-1,bi,bj))**2)
95 jmc 1.18
96     Alth2=sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)
97     & *(rA(i,j,bi,bj)**1.5)
98     Alth4=sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)
99     & *0.125*(rA(i,j,bi,bj)**2.5)
100 baylor 1.17 ELSE
101 jmc 1.18 C but this approximation will work on cube
102 baylor 1.17 c (and differs by as much as 4X)
103     grdVrt=abs((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))
104     grdVrt=max(grdVrt,
105     & abs((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj)))
106     grdVrt=max(grdVrt,
107     & abs((vort3(i+1,j+1)-vort3(i,j+1))*recip_DXG(i,j+1,bi,bj)))
108     grdVrt=max(grdVrt,
109     & abs((vort3(i+1,j+1)-vort3(i+1,j))*recip_DYG(i+1,j,bi,bj)))
110 jmc 1.18
111     Alth2=viscC2leith*grdVrt*(rA(i,j,bi,bj)**1.5)
112     Alth4=viscC4leith*grdVrt*0.125*(rA(i,j,bi,bj)**2.5)
113 baylor 1.17 ENDIF
114    
115     C Harmonic Modified Leith on Div.u points
116 jmc 1.12 Alin=viscAhD+vg2*rA ( i , j ,bi,bj)
117 jmc 1.18 viscAh_D(i,j)=min(viscAhMax,Alin+Alth2)
118 baylor 1.17 IF (vg2Max.GT.0.) THEN
119     AlinMax=vg2Max*(rA ( i , j ,bi,bj))
120     viscAh_D(i,j)=min(AlinMax,viscAh_D(i,j))
121     ENDIF
122     AlinMin=vg2Min*(rA ( i , j ,bi,bj))
123     viscAh_D(i,j)=max(AlinMin,viscAh_D(i,j))
124    
125     C BiHarmonic Modified Leith on Div.u points
126 jmc 1.12 Alin=viscA4D+vg4*(rA ( i , j ,bi,bj)**2)
127 jmc 1.18 viscA4_D(i,j)=min(viscA4Max,Alin+Alth4)
128 dimitri 1.13 IF (vg4Max.GT.0.) THEN
129     AlinMax=vg4Max*(rA ( i , j ,bi,bj)**2)
130     viscA4_D(i,j)=min(AlinMax,viscA4_D(i,j))
131     ENDIF
132     AlinMin=vg4Min*(rA ( i , j ,bi,bj)**2)
133     viscA4_D(i,j)=max(AlinMin,viscA4_D(i,j))
134 adcroft 1.6
135 baylor 1.17 C This is the vector magnitude of the vorticity gradient squared
136 jmc 1.18 IF (useFullLeith) THEN
137 baylor 1.17 grdVrt=0.25*(
138     & ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2
139     & +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2
140     & +((vort3(i-1,j)-vort3(i,j))*recip_DXG(i-1,j,bi,bj))**2
141     & +((vort3(i,j-1)-vort3(i,j))*recip_DYG(i,j-1,bi,bj))**2)
142    
143     C This is the vector magnitude of grad(div.v) squared
144     grdDiv=0.25*(
145     & ((hDiv(i+1,j)-hDiv(i,j))*recip_DXG(i,j,bi,bj))**2
146     & +((hDiv(i,j+1)-hDiv(i,j))*recip_DYG(i,j,bi,bj))**2
147     & +((hDiv(i+1,j+1)-hDiv(i,j+1))*recip_DXG(i,j+1,bi,bj))**2
148     & +((hDiv(i+1,j+1)-hDiv(i+1,j))*recip_DYG(i+1,j,bi,bj))**2)
149 jmc 1.18
150     Alth2=sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)
151     & *(rAz(i,j,bi,bj)**1.5)
152     Alth4=sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)
153     & *0.125*(rAz(i,j,bi,bj)**2.5)
154 baylor 1.17 ELSE
155 adcroft 1.6 C but this approximation will work on cube (and differs by as much as 4X)
156 baylor 1.17 grdVrt=abs((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))
157     grdVrt=max(grdVrt,
158     & abs((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj)))
159     grdVrt=max(grdVrt,
160     & abs((vort3(i-1,j)-vort3(i,j))*recip_DXG(i-1,j,bi,bj)))
161     grdVrt=max(grdVrt,
162     & abs((vort3(i,j-1)-vort3(i,j))*recip_DYG(i,j-1,bi,bj)))
163 jmc 1.18
164     Alth2=viscC2leith*grdVrt*(rAz(i,j,bi,bj)**1.5)
165     Alth4=viscC4leith*grdVrt*0.125*(rAz(i,j,bi,bj)**2.5)
166 baylor 1.17 ENDIF
167    
168     C Harmonic Modified Leith on Zeta points
169 jmc 1.12 Alin=viscAhZ+vg2*rAz( i , j ,bi,bj)
170 jmc 1.18 viscAh_Z(i,j)=min(viscAhMax,Alin+Alth2)
171 baylor 1.17 IF (vg2Max.GT.0.) THEN
172     AlinMax=vg2Max*(rAz( i , j ,bi,bj))
173     viscAh_Z(i,j)=min(AlinMax,viscAh_Z(i,j))
174     ENDIF
175 jmc 1.18 AlinMin=vg2Min*(rAz( i , j ,bi,bj))
176 baylor 1.17 viscAh_Z(i,j)=max(AlinMin,viscAh_Z(i,j))
177 adcroft 1.6
178 baylor 1.17 C BiHarmonic Modified Leith on Zeta Points
179 jmc 1.12 Alin=viscA4Z+vg4*(rAz( i , j ,bi,bj)**2)
180 jmc 1.18 viscA4_Z(i,j)=min(viscA4Max,Alin+Alth4)
181 dimitri 1.13 IF (vg4Max.GT.0.) THEN
182     AlinMax=vg4Max*(rAz( i , j ,bi,bj)**2)
183     viscA4_Z(i,j)=min(AlinMax,viscA4_Z(i,j))
184     ENDIF
185     AlinMin=vg4Min*(rAz( i , j ,bi,bj)**2)
186     viscA4_Z(i,j)=max(AlinMin,viscA4_Z(i,j))
187 adcroft 1.6 ENDDO
188     ENDDO
189     ELSE
190     DO j=1-Oly,sNy+Oly
191     DO i=1-Olx,sNx+Olx
192 jmc 1.12 viscAh_D(i,j)=viscAhD
193     viscAh_Z(i,j)=viscAhZ
194     viscA4_D(i,j)=viscA4D
195     viscA4_Z(i,j)=viscA4Z
196 adcroft 1.6 ENDDO
197     ENDDO
198     ENDIF
199 adcroft 1.2
200 jmc 1.12 C - Laplacian terms
201     IF ( viscC2leith.NE.0. .OR. viscAhGrid.NE.0.
202     & .OR. viscAhD.NE.0. .OR. viscAhZ.NE.0. ) THEN
203     DO j=2-Oly,sNy+Oly-1
204     DO i=2-Olx,sNx+Olx-1
205    
206     Dim=hDiv( i ,j-1)
207     Dij=hDiv( i , j )
208     Dmj=hDiv(i-1, j )
209     Zip=hFacZ( i ,j+1)*vort3( i ,j+1)
210     Zij=hFacZ( i , j )*vort3( i , j )
211     Zpj=hFacZ(i+1, j )*vort3(i+1, j )
212 adcroft 1.2
213 adcroft 1.4 C This bit scales the harmonic dissipation operator to be proportional
214     C to the grid-cell area over the time-step. viscAh is then non-dimensional
215     C and should be less than 1/8, for example viscAh=0.01
216 jmc 1.12 IF (useVariableViscosity) THEN
217 adcroft 1.6 Dij=Dij*viscAh_D(i,j)
218     Dim=Dim*viscAh_D(i,j-1)
219     Dmj=Dmj*viscAh_D(i-1,j)
220     Zij=Zij*viscAh_Z(i,j)
221     Zip=Zip*viscAh_Z(i,j+1)
222     Zpj=Zpj*viscAh_Z(i+1,j)
223 adcroft 1.4 uD2 = (
224     & cosFacU(j,bi,bj)*( Dij-Dmj )*recip_DXC(i,j,bi,bj)
225     & -recip_hFacW(i,j,k,bi,bj)*( Zip-Zij )*recip_DYG(i,j,bi,bj) )
226     vD2 = (
227     & recip_hFacS(i,j,k,bi,bj)*( Zpj-Zij )*recip_DXG(i,j,bi,bj)
228     & *cosFacV(j,bi,bj)
229     & +( Dij-Dim )*recip_DYC(i,j,bi,bj) )
230 jmc 1.12 ELSE
231     uD2 = viscAhD*
232 adcroft 1.2 & cosFacU(j,bi,bj)*( Dij-Dmj )*recip_DXC(i,j,bi,bj)
233 jmc 1.12 & - viscAhZ*recip_hFacW(i,j,k,bi,bj)*
234     & ( Zip-Zij )*recip_DYG(i,j,bi,bj)
235     vD2 = viscAhZ*recip_hFacS(i,j,k,bi,bj)*
236     & cosFacV(j,bi,bj)*( Zpj-Zij )*recip_DXG(i,j,bi,bj)
237     & + viscAhD* ( Dij-Dim )*recip_DYC(i,j,bi,bj)
238     ENDIF
239    
240     uDissip(i,j) = uD2
241     vDissip(i,j) = vD2
242    
243     ENDDO
244     ENDDO
245     ELSE
246     DO j=2-Oly,sNy+Oly-1
247     DO i=2-Olx,sNx+Olx-1
248     uDissip(i,j) = 0.
249     vDissip(i,j) = 0.
250     ENDDO
251     ENDDO
252     ENDIF
253    
254     C - Bi-harmonic terms
255     IF ( viscC4leith.NE.0. .OR. viscA4Grid.NE.0.
256     & .OR. viscA4D.NE.0. .OR. viscA4Z.NE.0. ) THEN
257     DO j=2-Oly,sNy+Oly-1
258     DO i=2-Olx,sNx+Olx-1
259 adcroft 1.2
260 jmc 1.11 #ifdef MOM_VI_ORIGINAL_VISCA4
261 jmc 1.12 Dim=dyF( i ,j-1,bi,bj)*dStar( i ,j-1)
262     Dij=dyF( i , j ,bi,bj)*dStar( i , j )
263     Dmj=dyF(i-1, j ,bi,bj)*dStar(i-1, j )
264    
265     Zip=dxV( i ,j+1,bi,bj)*hFacZ( i ,j+1)*zStar( i ,j+1)
266     Zij=dxV( i , j ,bi,bj)*hFacZ( i , j )*zStar( i , j )
267     Zpj=dxV(i+1, j ,bi,bj)*hFacZ(i+1, j )*zStar(i+1, j )
268 jmc 1.11 #else
269 jmc 1.12 Dim=dStar( i ,j-1)
270     Dij=dStar( i , j )
271     Dmj=dStar(i-1, j )
272    
273     Zip=hFacZ( i ,j+1)*zStar( i ,j+1)
274     Zij=hFacZ( i , j )*zStar( i , j )
275     Zpj=hFacZ(i+1, j )*zStar(i+1, j )
276 jmc 1.11 #endif
277    
278 adcroft 1.4 C This bit scales the harmonic dissipation operator to be proportional
279     C to the grid-cell area over the time-step. viscAh is then non-dimensional
280     C and should be less than 1/8, for example viscAh=0.01
281 jmc 1.12 IF (useVariableViscosity) THEN
282 adcroft 1.6 Dij=Dij*viscA4_D(i,j)
283     Dim=Dim*viscA4_D(i,j-1)
284     Dmj=Dmj*viscA4_D(i-1,j)
285     Zij=Zij*viscA4_Z(i,j)
286     Zip=Zip*viscA4_Z(i,j+1)
287     Zpj=Zpj*viscA4_Z(i+1,j)
288 jmc 1.11
289     #ifdef MOM_VI_ORIGINAL_VISCA4
290 adcroft 1.4 uD4 = recip_rAw(i,j,bi,bj)*(
291     & ( (Dij-Dmj)*cosFacU(j,bi,bj) )
292     & -recip_hFacW(i,j,k,bi,bj)*( Zip-Zij ) )
293     vD4 = recip_rAs(i,j,bi,bj)*(
294     & recip_hFacS(i,j,k,bi,bj)*( (Zpj-Zij)*cosFacV(j,bi,bj) )
295     & + ( Dij-Dim ) )
296 jmc 1.12 ELSE
297 adcroft 1.4 uD4 = recip_rAw(i,j,bi,bj)*(
298 adcroft 1.2 & viscA4*( (Dij-Dmj)*cosFacU(j,bi,bj) )
299     & -recip_hFacW(i,j,k,bi,bj)*viscA4*( Zip-Zij ) )
300 adcroft 1.4 vD4 = recip_rAs(i,j,bi,bj)*(
301 adcroft 1.2 & recip_hFacS(i,j,k,bi,bj)*viscA4*( (Zpj-Zij)*cosFacV(j,bi,bj) )
302     & + viscA4*( Dij-Dim ) )
303 jmc 1.11 #else /* MOM_VI_ORIGINAL_VISCA4 */
304     uD4 = (
305     & cosFacU(j,bi,bj)*( Dij-Dmj )*recip_DXC(i,j,bi,bj)
306     & -recip_hFacW(i,j,k,bi,bj)*( Zip-Zij )*recip_DYG(i,j,bi,bj) )
307     vD4 = (
308     & recip_hFacS(i,j,k,bi,bj)*( Zpj-Zij )*recip_DXG(i,j,bi,bj)
309     & *cosFacV(j,bi,bj)
310     & +( Dij-Dim )*recip_DYC(i,j,bi,bj) )
311 jmc 1.12 ELSE
312 jmc 1.18 uD4 = viscA4D*
313 jmc 1.11 & cosFacU(j,bi,bj)*( Dij-Dmj )*recip_DXC(i,j,bi,bj)
314 jmc 1.12 & - viscA4Z*recip_hFacW(i,j,k,bi,bj)*
315     & ( Zip-Zij )*recip_DYG(i,j,bi,bj)
316     vD4 = viscA4Z*recip_hFacS(i,j,k,bi,bj)*
317     & cosFacV(j,bi,bj)*( Zpj-Zij )*recip_DXG(i,j,bi,bj)
318     & + viscA4D* ( Dij-Dim )*recip_DYC(i,j,bi,bj)
319 jmc 1.11 #endif /* MOM_VI_ORIGINAL_VISCA4 */
320 jmc 1.12 ENDIF
321 adcroft 1.2
322 jmc 1.12 uDissip(i,j) = uDissip(i,j) - uD4
323     vDissip(i,j) = vDissip(i,j) - vD4
324 adcroft 1.2
325 jmc 1.12 ENDDO
326 adcroft 1.2 ENDDO
327 jmc 1.12 ENDIF
328 molod 1.7
329     #ifdef ALLOW_DIAGNOSTICS
330 jmc 1.15 IF (useDiagnostics) THEN
331 jmc 1.16 CALL DIAGNOSTICS_FILL(viscAh_D,'VISCAH ',k,1,2,bi,bj,myThid)
332     CALL DIAGNOSTICS_FILL(viscA4_D,'VISCA4 ',k,1,2,bi,bj,myThid)
333 jmc 1.15 ENDIF
334 molod 1.7 #endif
335 adcroft 1.2
336     RETURN
337     END

  ViewVC Help
Powered by ViewVC 1.1.22