/[MITgcm]/MITgcm/pkg/mom_common/mom_calc_visc.F
ViewVC logotype

Diff of /MITgcm/pkg/mom_common/mom_calc_visc.F

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

revision 1.15 by jmc, Mon Oct 3 21:43:03 2005 UTC revision 1.16 by jmc, Tue Oct 4 02:49:13 2005 UTC
# Line 95  C     == Local variables == Line 95  C     == Local variables ==
95        _RL Alin,grdVrt,grdDiv, keZpt        _RL Alin,grdVrt,grdDiv, keZpt
96        _RL recip_dt,L2,L3,L4,L5,L2rdt,L4rdt        _RL recip_dt,L2,L3,L4,L5,L2rdt,L4rdt
97        _RL Uscl,U4scl        _RL Uscl,U4scl
98          _RL divDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
99          _RL divDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
100        _RL viscAh_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscAh_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
101        _RL viscAh_DMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscAh_DMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102        _RL viscA4_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscA4_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 183  C     == Local variables == Line 185  C     == Local variables ==
185    
186  C     - Viscosity  C     - Viscosity
187        IF (useVariableViscosity) THEN        IF (useVariableViscosity) THEN
188    
189    C      horizontal gradient of horizontal divergence:
190           DO j=1-Oly,sNy+Oly
191            DO i=1-Olx,sNx+Olx
192              divDx(i,j) = 0.
193              divDy(i,j) = 0.
194            ENDDO
195           ENDDO
196           IF (calcleith) THEN
197    C-       gradient in x direction:
198    #ifndef ALLOW_AUTODIFF_TAMC
199             IF (useCubedSphereExchange) THEN
200    C        to compute d/dx(hDiv), fill corners with appropriate values:
201               CALL FILL_CS_CORNER_TR_RL( .TRUE., hDiv, bi,bj, myThid )
202             ENDIF
203    #endif
204             DO j=2-Oly,sNy+Oly-1
205              DO i=2-Olx,sNx+Olx-1
206                divDx(i,j) = (hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj)
207              ENDDO
208             ENDDO
209    
210    C-       gradient in y direction:
211    #ifndef ALLOW_AUTODIFF_TAMC
212             IF (useCubedSphereExchange) THEN
213    C        to compute d/dy(hDiv), fill corners with appropriate values:
214               CALL FILL_CS_CORNER_TR_RL(.FALSE., hDiv, bi,bj, myThid )
215             ENDIF
216    #endif
217             DO j=2-Oly,sNy+Oly-1
218              DO i=2-Olx,sNx+Olx-1
219                divDy(i,j) = (hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj)
220              ENDDO
221             ENDDO
222           ENDIF
223    
224         DO j=2-Oly,sNy+Oly-1         DO j=2-Oly,sNy+Oly-1
225          DO i=2-Olx,sNx+Olx-1          DO i=2-Olx,sNx+Olx-1
226  CCCCCCCCCCCCCCC Divergence Point CalculationsCCCCCCCCCCCCCCCCCCCC  CCCCCCCCCCCCCCC Divergence Point CalculationsCCCCCCCCCCCCCCCCCCCC
# Line 232  C This is the vector magnitude of the vo Line 270  C This is the vector magnitude of the vo
270    
271  C This is the vector magnitude of grad (div.v) squared  C This is the vector magnitude of grad (div.v) squared
272  C Using it in Leith serves to damp instabilities in w.  C Using it in Leith serves to damp instabilities in w.
273            grdDiv=0.25 _d 0*(            grdDiv=0.25 _d 0*( (divDx(i+1,j)*divDx(i+1,j)
274       &     ((hDiv(i+1,j)-hDiv(i,j))*recip_DXC(i+1,j,bi,bj))**2       &                        + divDx(i,j)*divDx(i,j) )
275       &     +((hDiv(i,j+1)-hDiv(i,j))*recip_DYC(i,j+1,bi,bj))**2       &                     + (divDy(i,j+1)*divDy(i,j+1)
276       &     +((hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj))**2       &                        + divDy(i,j)*divDy(i,j) )  )
      &     +((hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj))**2)  
277    
278            viscAh_DLth(i,j)=            viscAh_DLth(i,j)=
279       &     sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3       &     sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3
# Line 257  c (and differs by as much as 4X) Line 294  c (and differs by as much as 4X)
294            grdVrt=max(grdVrt,            grdVrt=max(grdVrt,
295       &     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)))
296    
297            grdDiv=abs((hDiv(i+1,j)-hDiv(i,j))*recip_DXC(i+1,j,bi,bj))            grdDiv=max( abs(divDx(i+1,j)), abs(divDx(i,j)) )
298            grdDiv=max(grdDiv,            grdDiv=max( grdDiv, abs(divDy(i,j+1)) )
299       &     abs((hDiv(i,j+1)-hDiv(i,j))*recip_DYC(i,j+1,bi,bj)))            grdDiv=max( grdDiv, abs(divDy(i,j))   )
           grdDiv=max(grdDiv,  
      &     abs((hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj)))  
           grdDiv=max(grdDiv,  
      &     abs((hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj)))  
300    
301  c This approximation is good to the same order as above...  c This approximation is good to the same order as above...
302            viscAh_Dlth(i,j)=            viscAh_Dlth(i,j)=
# Line 355  C This is the vector magnitude of the vo Line 388  C This is the vector magnitude of the vo
388       &     +((vort3(i,j-1)-vort3(i,j))*recip_DYG(i,j-1,bi,bj))**2)       &     +((vort3(i,j-1)-vort3(i,j))*recip_DYG(i,j-1,bi,bj))**2)
389    
390  C This is the vector magnitude of grad(div.v) squared  C This is the vector magnitude of grad(div.v) squared
391            grdDiv=0.25 _d 0*(            grdDiv=0.25 _d 0*( (divDx(i,j-1)*divDx(i,j-1)
392       &     ((hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj))**2       &                        + divDx(i,j)*divDx(i,j) )
393       &     +((hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj))**2       &                     + (divDy(i-1,j)*divDy(i-1,j)
394       &     +((hDiv(i,j-1)-hDiv(i-1,j-1))*recip_DXC(i,j-1,bi,bj))**2       &                        + divDy(i,j)*divDy(i,j) )  )
      &     +((hDiv(i-1,j)-hDiv(i-1,j-1))*recip_DYC(i-1,j,bi,bj))**2)  
395    
396            viscAh_ZLth(i,j)=            viscAh_ZLth(i,j)=
397       &     sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3       &     sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3
# Line 380  C but this approximation will work on cu Line 412  C but this approximation will work on cu
412            grdVrt=max(grdVrt,            grdVrt=max(grdVrt,
413       &     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)))
414    
415            grdDiv=abs((hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj))            grdDiv=max( abs(divDx(i,j)), abs(divDx(i,j-1)) )
416            grdDiv=max(grdDiv,            grdDiv=max( grdDiv, abs(divDy(i,j))   )
417       &     abs((hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj)))            grdDiv=max( grdDiv, abs(divDy(i-1,j)) )
           grdDiv=max(grdDiv,  
      &     abs((hDiv(i,j-1)-hDiv(i-1,j-1))*recip_DXC(i,j-1,bi,bj)))  
           grdDiv=max(grdDiv,  
      &     abs((hDiv(i-1,j)-hDiv(i-1,j-1))*recip_DYC(i-1,j,bi,bj)))  
418    
419            viscAh_ZLth(i,j)=(viscC2leith*grdVrt            viscAh_ZLth(i,j)=(viscC2leith*grdVrt
420       &                     +(viscC2leithD*grdDiv))*L3       &                     +(viscC2leithD*grdDiv))*L3
Legend:
Removed from v.1.15  
changed lines
  Added in v.1.16
  ViewVC Help
Powered by ViewVC 1.1.22