/[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.16 by jmc, Mon Dec 20 19:08:21 2004 UTC revision 1.17 by baylor, Thu Mar 10 02:39:56 2005 UTC
# Line 45  C     == Routine arguments == Line 45  C     == Routine arguments ==
45  C     == Local variables ==  C     == Local variables ==
46        INTEGER I,J        INTEGER I,J
47        _RL Zip,Zij,Zpj,Dim,Dij,Dmj,uD2,vD2,uD4,vD4        _RL Zip,Zij,Zpj,Dim,Dij,Dmj,uD2,vD2,uD4,vD4
48        _RL Alin,AlinMin,AlinMax,Alth,grdVrt,vg2,vg4,vg4Min,vg4Max        _RL Alin,AlinMin,AlinMax,Alth,grdVrt,grdDiv
49          _RL vg2,vg2Min,vg2Max,vg4,vg4Min,vg4Max
50        LOGICAL useVariableViscosity        LOGICAL useVariableViscosity
51    
52        useVariableViscosity=        useVariableViscosity=
53       &      (viscAhGrid*deltaTmom.NE.0.)       &      (viscAhGrid*deltaTmom.NE.0.)
54       &  .OR.(viscA4Grid*deltaTmom.NE.0.)       &  .OR.(viscA4Grid*deltaTmom.NE.0.)
55       &  .OR.(viscC2leith.NE.0.)       &  .OR.(viscC2leith.NE.0.)
56         &  .OR.(viscC2leithD.NE.0.)
57       &  .OR.(viscC4leith.NE.0.)       &  .OR.(viscC4leith.NE.0.)
58         &  .OR.(viscC4leithD.NE.0.)
59    
60        IF (deltaTmom.NE.0.) THEN        IF (deltaTmom.NE.0.) THEN
61         vg2=viscAhGrid/deltaTmom         vg2=viscAhGrid/deltaTmom
62           vg2Min=viscAhGridMin/deltaTmom
63           vg2Max=viscAhGridMax/deltaTmom
64         vg4=viscA4Grid/deltaTmom         vg4=viscA4Grid/deltaTmom
65         vg4Min=viscA4GridMin/deltaTmom         vg4Min=viscA4GridMin/deltaTmom
66         vg4Max=viscA4GridMax/deltaTmom         vg4Max=viscA4GridMax/deltaTmom
67        ELSE        ELSE
68         vg2=0.         vg2=0.
69           vg2Min=0.
70           vg2Max=0.
71         vg4=0.         vg4=0.
72         vg4Min=0.         vg4Min=0.
73         vg4Max=0.         vg4Max=0.
# Line 70  C     - Viscosity Line 77  C     - Viscosity
77        IF (useVariableViscosity) THEN        IF (useVariableViscosity) THEN
78         DO j=2-Oly,sNy+Oly-1         DO j=2-Oly,sNy+Oly-1
79          DO i=2-Olx,sNx+Olx-1          DO i=2-Olx,sNx+Olx-1
80  C This is the vector magnitude of the vorticity gradient           IF (useFullLeith) THEN
81  c        grdVrt=sqrt(0.25*(  C This is the vector magnitude of the vorticity gradient squared
82  c    &      ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2            grdVrt=0.25*(
83  c    &     +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2       &     ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2
84  c    &     +((vort3(i+1,j+1)-vort3(i,j+1))*recip_DXG(i,j+1,bi,bj))**2       &     +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2
85  c    &     +((vort3(i+1,j+1)-vort3(i+1,j))*recip_DYG(i+1,j,bi,bj))**2 ))       &     +((vort3(i+1,j+1)-vort3(i,j+1))*recip_DXG(i,j+1,bi,bj))**2
86  C but this approximation will work on cube (and differs by as much as 4X)       &     +((vort3(i+1,j+1)-vort3(i+1,j))*recip_DYG(i+1,j,bi,bj))**2)
87           grdVrt=abs((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))  
88           grdVrt=max(grdVrt,  C This is the vector magnitude of grad (div.v) squared
89       &        abs((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj)))  C Using it in Leith serves to damp instabilities in w.
90           grdVrt=max(grdVrt,             grdDiv=0.25*(
91       &        abs((vort3(i+1,j+1)-vort3(i,j+1))*recip_DXG(i,j+1,bi,bj)))       &      ((hDiv(i+1,j)-hDiv(i,j))*recip_DXG(i,j,bi,bj))**2
92           grdVrt=max(grdVrt,       &      +((hDiv(i,j+1)-hDiv(i,j))*recip_DYG(i,j,bi,bj))**2
93       &        abs((vort3(i+1,j+1)-vort3(i+1,j))*recip_DYG(i+1,j,bi,bj)))       &      +((hDiv(i-1,j)-hDiv(i,j))*recip_DXG(i-1,j,bi,bj))**2
94           Alth=viscC2leith*grdVrt*(rA(i,j,bi,bj)**1.5)       &      +((hDiv(i,j-1)-hDiv(i,j))*recip_DYG(i,j-1,bi,bj))**2)
95             ELSE
96    C but this approximation will work on cube
97    c (and differs by as much as 4X)
98               grdVrt=abs((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))
99               grdVrt=max(grdVrt,
100         &      abs((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj)))
101               grdVrt=max(grdVrt,
102         &      abs((vort3(i+1,j+1)-vort3(i,j+1))*recip_DXG(i,j+1,bi,bj)))
103               grdVrt=max(grdVrt,
104         &      abs((vort3(i+1,j+1)-vort3(i+1,j))*recip_DYG(i+1,j,bi,bj)))
105               grdVrt=grdVrt**2
106               grdDiv=0.
107             ENDIF
108    
109    C  Harmonic Modified Leith on Div.u points
110             Alth=sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)
111         &         *(rA(i,j,bi,bj)**1.5)
112           Alin=viscAhD+vg2*rA ( i , j ,bi,bj)           Alin=viscAhD+vg2*rA ( i , j ,bi,bj)
113           viscAh_D(i,j)=min(viscAhMax,Alin+Alth)           viscAh_D(i,j)=min(viscAhMax,Alin+Alth)
114           Alth=viscC4leith*grdVrt*0.125*(rA(i,j,bi,bj)**2.5)           IF (vg2Max.GT.0.) THEN
115                AlinMax=vg2Max*(rA ( i , j ,bi,bj))
116                viscAh_D(i,j)=min(AlinMax,viscAh_D(i,j))
117             ENDIF
118             AlinMin=vg2Min*(rA ( i , j ,bi,bj))
119             viscAh_D(i,j)=max(AlinMin,viscAh_D(i,j))
120    
121    C  BiHarmonic Modified Leith on Div.u points
122             Alth=sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)
123         &          *0.125*(rAz(i,j,bi,bj)**2.5)
124           Alin=viscA4D+vg4*(rA ( i , j ,bi,bj)**2)           Alin=viscA4D+vg4*(rA ( i , j ,bi,bj)**2)
125           viscA4_D(i,j)=min(viscA4Max,Alin+Alth)           viscA4_D(i,j)=min(viscA4Max,Alin+Alth)
126           IF (vg4Max.GT.0.) THEN           IF (vg4Max.GT.0.) THEN
# Line 97  C but this approximation will work on cu Line 130  C but this approximation will work on cu
130           AlinMin=vg4Min*(rA ( i , j ,bi,bj)**2)           AlinMin=vg4Min*(rA ( i , j ,bi,bj)**2)
131           viscA4_D(i,j)=max(AlinMin,viscA4_D(i,j))           viscA4_D(i,j)=max(AlinMin,viscA4_D(i,j))
132    
133  C This is the vector magnitude of the vorticity gradient  C This is the vector magnitude of the vorticity gradient squared
134  c        grdVrt=sqrt(0.25*(           IF (useFullLeith) THEN
135  c    &      ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2             grdVrt=0.25*(
136  c    &     +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2       &      ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2
137  c    &     +((vort3(i-1,j)-vort3(i,j))*recip_DXG(i-1,j,bi,bj))**2       &      +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2
138  c    &     +((vort3(i,j-1)-vort3(i,j))*recip_DYG(i,j-1,bi,bj))**2 ))       &      +((vort3(i-1,j)-vort3(i,j))*recip_DXG(i-1,j,bi,bj))**2
139         &      +((vort3(i,j-1)-vort3(i,j))*recip_DYG(i,j-1,bi,bj))**2)
140    
141    C This is the vector magnitude of grad(div.v) squared
142               grdDiv=0.25*(
143         &       ((hDiv(i+1,j)-hDiv(i,j))*recip_DXG(i,j,bi,bj))**2
144         &      +((hDiv(i,j+1)-hDiv(i,j))*recip_DYG(i,j,bi,bj))**2
145         &      +((hDiv(i+1,j+1)-hDiv(i,j+1))*recip_DXG(i,j+1,bi,bj))**2
146         &      +((hDiv(i+1,j+1)-hDiv(i+1,j))*recip_DYG(i+1,j,bi,bj))**2)
147             ELSE
148  C but this approximation will work on cube (and differs by as much as 4X)  C but this approximation will work on cube (and differs by as much as 4X)
149           grdVrt=abs((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))             grdVrt=abs((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))
150           grdVrt=max(grdVrt,             grdVrt=max(grdVrt,
151       &          abs((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj)))       &       abs((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj)))
152           grdVrt=max(grdVrt,             grdVrt=max(grdVrt,
153       &          abs((vort3(i-1,j)-vort3(i,j))*recip_DXG(i-1,j,bi,bj)))       &       abs((vort3(i-1,j)-vort3(i,j))*recip_DXG(i-1,j,bi,bj)))
154           grdVrt=max(grdVrt,             grdVrt=max(grdVrt,
155       &          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)))
156           Alth=viscC2leith*grdVrt*(rAz(i,j,bi,bj)**1.5)             grdVrt=grdVrt**2
157               grdDiv=0.
158             ENDIF
159    
160    C  Harmonic Modified Leith on Zeta points
161             Alth=sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)
162         &         *(rA(i,j,bi,bj)**1.5)
163           Alin=viscAhZ+vg2*rAz( i , j ,bi,bj)           Alin=viscAhZ+vg2*rAz( i , j ,bi,bj)
164           viscAh_Z(i,j)=min(viscAhMax,Alin+Alth)           viscAh_Z(i,j)=min(viscAhMax,Alin+Alth)
165             IF (vg2Max.GT.0.) THEN
166                AlinMax=vg2Max*(rAz( i , j ,bi,bj))
167                viscAh_Z(i,j)=min(AlinMax,viscAh_Z(i,j))
168             ENDIF
169             AlinMin=vg2Min*(rA ( i , j ,bi,bj))
170             viscAh_Z(i,j)=max(AlinMin,viscAh_Z(i,j))
171    
172           Alth=viscC4leith*grdVrt*0.125*(rAz(i,j,bi,bj)**2.5)  C  BiHarmonic Modified Leith on Zeta Points
173             Alth=sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)
174         &          *0.125*(rAz(i,j,bi,bj)**2.5)
175             grdVrt=sqrt(grdVrt)
176           Alin=viscA4Z+vg4*(rAz( i , j ,bi,bj)**2)           Alin=viscA4Z+vg4*(rAz( i , j ,bi,bj)**2)
177           viscA4_Z(i,j)=min(viscA4Max,Alin+Alth)           viscA4_Z(i,j)=min(viscA4Max,Alin+Alth)
178           IF (vg4Max.GT.0.) THEN           IF (vg4Max.GT.0.) THEN
Legend:
Removed from v.1.16  
changed lines
  Added in v.1.17
  ViewVC Help
Powered by ViewVC 1.1.22