/[MITgcm]/MITgcm_contrib/SOSE/code_ad/mom_vi_hdissip.F
ViewVC logotype

Annotation of /MITgcm_contrib/SOSE/code_ad/mom_vi_hdissip.F

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


Revision 1.1 - (hide annotations) (download)
Fri Apr 23 19:55:12 2010 UTC (15 years, 3 months ago) by mmazloff
Branch: MAIN
CVS Tags: HEAD
original files

1 mmazloff 1.1 C $Header: /u/gcmpack/MITgcm/pkg/mom_vecinv/mom_vi_hdissip.F,v 1.28 2005/09/26 15:27:11 baylor Exp $
2     C $Name: $
3    
4     #include "MOM_VECINV_OPTIONS.h"
5    
6     SUBROUTINE MOM_VI_HDISSIP(
7     I bi,bj,k,
8     I hDiv,vort3,tension,strain,KE,
9     I hFacZ,dStar,zStar,
10     I viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,
11     I harmonic,biharmonic,useVariableViscosity,
12     O uDissip,vDissip,
13     I myThid)
14    
15     cph(
16     cph The following line was commented in the argument list
17     cph TAMC cannot digest commented lines within continuing lines
18     c I viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,
19     cph)
20    
21     IMPLICIT NONE
22     C
23     C Calculate horizontal dissipation terms
24     C [del^2 - del^4] (u,v)
25     C
26    
27     C == Global variables ==
28     #include "SIZE.h"
29     #include "GRID.h"
30     #include "EEPARAMS.h"
31     #include "PARAMS.h"
32    
33     C == Routine arguments ==
34     INTEGER bi,bj,k
35     _RL hDiv(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
36     _RL vort3(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
37     _RL tension(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
38     _RL strain(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
39     _RL KE(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
40     _RS hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
41     _RL dStar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
42     _RL zStar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
43     _RL uDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
44     _RL vDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
45     _RL viscAh_Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
46     _RL viscAh_D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
47     _RL viscA4_Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
48     _RL viscA4_D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
49     LOGICAL harmonic, biharmonic, useVariableViscosity
50     INTEGER myThid
51    
52     C == Local variables ==
53     INTEGER I,J
54     _RL Zip,Zij,Zpj,Dim,Dij,Dmj,uD2,vD2,uD4,vD4
55     CMM(
56     INTEGER mm,klevbm,kplace
57     _RL SclFactD,SclFactZ
58     _RL SclFactD4,SclFactZ4
59     CMM( Saniiity check
60     CMM print*,'CMM=> rf(k),r_low',rf(k),r_low(i,j,bi,bj)
61     CMM print*,'CMM=>k,klevbm,kplace,ViscAh',k,klevbm,kplace,SclFactD
62     CMM)CMM)
63    
64     C - Laplacian terms
65     IF (harmonic) THEN
66     DO j=2-Oly,sNy+Oly-1
67     DO i=2-Olx,sNx+Olx-1
68     CMM(
69     klevbm = 0
70     do mm = 1,Nr
71     if (rf(mm).GE.r_low(i,j,bi,bj)) then
72     klevbm = mm
73     endif
74     enddo
75     kplace=klevbm-k
76     CMM kplace = 0 is first wet cell, 1 is second, etc
77     Cmm kplace GT 1 excludes bottom 2 cells
78     IF (kplace.GT.2) THEN
79     SclFactD = viscAhD
80     SclFactZ = viscAhZ
81     ELSE
82     SclFactD = 1000.D0
83     SclFactZ = 1000.D0
84     ENDIF
85     CMM)
86     Dim=hDiv( i ,j-1)
87     Dij=hDiv( i , j )
88     Dmj=hDiv(i-1, j )
89     Zip=hFacZ( i ,j+1)*vort3( i ,j+1)
90     Zij=hFacZ( i , j )*vort3( i , j )
91     Zpj=hFacZ(i+1, j )*vort3(i+1, j )
92    
93     C This bit scales the harmonic dissipation operator to be proportional
94     C to the grid-cell area over the time-step. viscAh is then non-dimensional
95     C and should be less than 1/8, for example viscAh=0.01
96     IF (useVariableViscosity) THEN
97     Dij=Dij*viscAh_D(i,j)
98     Dim=Dim*viscAh_D(i,j-1)
99     Dmj=Dmj*viscAh_D(i-1,j)
100     Zij=Zij*viscAh_Z(i,j)
101     Zip=Zip*viscAh_Z(i,j+1)
102     Zpj=Zpj*viscAh_Z(i+1,j)
103     uD2 = (
104     & cosFacU(j,bi,bj)*( Dij-Dmj )*recip_DXC(i,j,bi,bj)
105     & -recip_hFacW(i,j,k,bi,bj)*( Zip-Zij )*recip_DYG(i,j,bi,bj) )
106     vD2 = (
107     & recip_hFacS(i,j,k,bi,bj)*( Zpj-Zij )*recip_DXG(i,j,bi,bj)
108     & *cosFacV(j,bi,bj)
109     & +( Dij-Dim )*recip_DYC(i,j,bi,bj) )
110     ELSE
111     uD2 =
112     CMM( viscAhD*
113     & SclFactD*
114     CMM)
115     & cosFacU(j,bi,bj)*( Dij-Dmj )*recip_DXC(i,j,bi,bj)
116     CMM( - viscAhZ*
117     & - SclFactZ*
118     CMM)
119     & recip_hFacW(i,j,k,bi,bj)*
120     & ( Zip-Zij )*recip_DYG(i,j,bi,bj)
121     vD2 =
122     CMM( viscAhZ*
123     & SclFactZ*
124     CMM)
125     & recip_hFacS(i,j,k,bi,bj)*
126     & cosFacV(j,bi,bj)*( Zpj-Zij )*recip_DXG(i,j,bi,bj)
127     CMM( & + viscAhD*
128     & + SclFactD*
129     CMM)
130     & ( Dij-Dim )*recip_DYC(i,j,bi,bj)
131     ENDIF
132    
133     uDissip(i,j) = uD2
134     vDissip(i,j) = vD2
135    
136     ENDDO
137     ENDDO
138     ELSE
139     DO j=2-Oly,sNy+Oly-1
140     DO i=2-Olx,sNx+Olx-1
141     uDissip(i,j) = 0.
142     vDissip(i,j) = 0.
143     ENDDO
144     ENDDO
145     ENDIF
146    
147     C - Bi-harmonic terms
148     IF (biharmonic) THEN
149     DO j=2-Oly,sNy+Oly-1
150     DO i=2-Olx,sNx+Olx-1
151     CMM(
152     klevbm = 0
153     do mm = 1,Nr
154     if (rf(mm).GE.r_low(i,j,bi,bj)) then
155     klevbm = mm
156     endif
157     enddo
158     kplace=klevbm-k
159     CMM kplace = 0 is first wet cell, 1 is second, etc
160     Cmm kplace GT 1 excludes bottom 2 cells
161     IF (kplace.GT.2) THEN
162     SclFactD4 = viscA4D
163     SclFactZ4 = viscA4Z
164     ELSE
165     SclFactD4 = 1000000000.D0
166     SclFactZ4 = 1000000000.D0
167     ENDIF
168     CMM)
169     #ifdef MOM_VI_ORIGINAL_VISCA4
170     Dim=dyF( i ,j-1,bi,bj)*dStar( i ,j-1)
171     Dij=dyF( i , j ,bi,bj)*dStar( i , j )
172     Dmj=dyF(i-1, j ,bi,bj)*dStar(i-1, j )
173    
174     Zip=dxV( i ,j+1,bi,bj)*hFacZ( i ,j+1)*zStar( i ,j+1)
175     Zij=dxV( i , j ,bi,bj)*hFacZ( i , j )*zStar( i , j )
176     Zpj=dxV(i+1, j ,bi,bj)*hFacZ(i+1, j )*zStar(i+1, j )
177     #else
178     Dim=dStar( i ,j-1)
179     Dij=dStar( i , j )
180     Dmj=dStar(i-1, j )
181    
182     Zip=hFacZ( i ,j+1)*zStar( i ,j+1)
183     Zij=hFacZ( i , j )*zStar( i , j )
184     Zpj=hFacZ(i+1, j )*zStar(i+1, j )
185     #endif
186    
187     C This bit scales the harmonic dissipation operator to be proportional
188     C to the grid-cell area over the time-step. viscAh is then non-dimensional
189     C and should be less than 1/8, for example viscAh=0.01
190     IF (useVariableViscosity) THEN
191     Dij=Dij*viscA4_D(i,j)
192     Dim=Dim*viscA4_D(i,j-1)
193     Dmj=Dmj*viscA4_D(i-1,j)
194     Zij=Zij*viscA4_Z(i,j)
195     Zip=Zip*viscA4_Z(i,j+1)
196     Zpj=Zpj*viscA4_Z(i+1,j)
197    
198     #ifdef MOM_VI_ORIGINAL_VISCA4
199     uD4 = recip_rAw(i,j,bi,bj)*(
200     & ( (Dij-Dmj)*cosFacU(j,bi,bj) )
201     & -recip_hFacW(i,j,k,bi,bj)*( Zip-Zij ) )
202     vD4 = recip_rAs(i,j,bi,bj)*(
203     & recip_hFacS(i,j,k,bi,bj)*( (Zpj-Zij)*cosFacV(j,bi,bj) )
204     & + ( Dij-Dim ) )
205     ELSE
206     uD4 = recip_rAw(i,j,bi,bj)*(
207     & viscA4*( (Dij-Dmj)*cosFacU(j,bi,bj) )
208     & -recip_hFacW(i,j,k,bi,bj)*viscA4*( Zip-Zij ) )
209     vD4 = recip_rAs(i,j,bi,bj)*(
210     & recip_hFacS(i,j,k,bi,bj)*viscA4*( (Zpj-Zij)*cosFacV(j,bi,bj) )
211     & + viscA4*( Dij-Dim ) )
212     #else /* MOM_VI_ORIGINAL_VISCA4 */
213     uD4 = (
214     & cosFacU(j,bi,bj)*( Dij-Dmj )*recip_DXC(i,j,bi,bj)
215     & -recip_hFacW(i,j,k,bi,bj)*( Zip-Zij )*recip_DYG(i,j,bi,bj) )
216     vD4 = (
217     & recip_hFacS(i,j,k,bi,bj)*( Zpj-Zij )*recip_DXG(i,j,bi,bj)
218     & *cosFacV(j,bi,bj)
219     & +( Dij-Dim )*recip_DYC(i,j,bi,bj) )
220     ELSE
221     uD4 = SclFactD4*
222     CMM( viscA4D*
223     & cosFacU(j,bi,bj)*( Dij-Dmj )*recip_DXC(i,j,bi,bj)
224     & - SclFactZ4*recip_hFacW(i,j,k,bi,bj)*
225     & ( Zip-Zij )*recip_DYG(i,j,bi,bj)
226     vD4 = SclFactZ4*recip_hFacS(i,j,k,bi,bj)*
227     & cosFacV(j,bi,bj)*( Zpj-Zij )*recip_DXG(i,j,bi,bj)
228     & + SclFactD4* ( Dij-Dim )*recip_DYC(i,j,bi,bj)
229     CMM)
230     #endif /* MOM_VI_ORIGINAL_VISCA4 */
231     ENDIF
232    
233     uDissip(i,j) = uDissip(i,j) - uD4
234     vDissip(i,j) = vDissip(i,j) - vD4
235    
236     ENDDO
237     ENDDO
238     ENDIF
239    
240     RETURN
241     END

  ViewVC Help
Powered by ViewVC 1.1.22