/[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.18 by baylor, Mon Oct 10 18:45:30 2005 UTC revision 1.20 by jmc, Wed Oct 12 20:24:22 2005 UTC
# Line 99  C     == Local variables == Line 99  C     == Local variables ==
99        _RL Uscl,U4scl        _RL Uscl,U4scl
100        _RL divDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL divDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
101        _RL divDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL divDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102          _RL vrtDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
103          _RL vrtDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
104        _RL viscAh_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscAh_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105        _RL viscAh_DMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscAh_DMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
106        _RL viscA4_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscA4_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 187  C     == Local variables == Line 189  C     == Local variables ==
189    
190        IF (calcleith) THEN        IF (calcleith) THEN
191          IF (useFullLeith) THEN          IF (useFullLeith) THEN
192           leith2fac=(viscC2leith/pi)**6           leith2fac =(viscC2leith /pi)**6
          leith4fac=0.015625 _d 0*(viscC4leith/pi)**6  
193           leithD2fac=(viscC2leithD/pi)**6           leithD2fac=(viscC2leithD/pi)**6
194             leith4fac =0.015625 _d 0*(viscC4leith /pi)**6
195           leithD4fac=0.015625 _d 0*(viscC4leithD/pi)**6           leithD4fac=0.015625 _d 0*(viscC4leithD/pi)**6
196          ELSE          ELSE
197           leith2fac=(viscC2leith/pi)**3           leith2fac =(viscC2leith /pi)**3
198           leithD2fac=(viscC2leithD/pi)**3           leithD2fac=(viscC2leithD/pi)**3
199           leith4fac=0.0125 _d 0*(viscC4leith/pi)**3           leith4fac =0.125 _d 0*(viscC4leith /pi)**3
200           leithD4fac=0.0125 _d 0*(viscC4leithD/pi)**3           leithD4fac=0.125 _d 0*(viscC4leithD/pi)**3
201          ENDIF          ENDIF
202        ELSE        ELSE
203          leith2fac=0. _d 0          leith2fac=0. _d 0
# Line 207  C     == Local variables == Line 209  C     == Local variables ==
209  C     - Viscosity  C     - Viscosity
210        IF (useVariableViscosity) THEN        IF (useVariableViscosity) THEN
211    
212  C      horizontal gradient of horizontal divergence:  C-     Initialise to zero gradient of vorticity & divergence:
213         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
214          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
215            divDx(i,j) = 0.            divDx(i,j) = 0.
216            divDy(i,j) = 0.            divDy(i,j) = 0.
217              vrtDx(i,j) = 0.
218              vrtDy(i,j) = 0.
219          ENDDO          ENDDO
220         ENDDO         ENDDO
221    
222         IF (calcleith) THEN         IF (calcleith) THEN
223    C      horizontal gradient of horizontal divergence:
224    
225  C-       gradient in x direction:  C-       gradient in x direction:
226  #ifndef ALLOW_AUTODIFF_TAMC  #ifndef ALLOW_AUTODIFF_TAMC
227           IF (useCubedSphereExchange) THEN           IF (useCubedSphereExchange) THEN
# Line 240  C        to compute d/dy(hDiv), fill cor Line 247  C        to compute d/dy(hDiv), fill cor
247              divDy(i,j) = (hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj)              divDy(i,j) = (hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj)
248            ENDDO            ENDDO
249           ENDDO           ENDDO
250    
251    C      horizontal gradient of vertical vorticity:
252    C-       gradient in x direction:
253             DO j=2-Oly,sNy+Oly
254              DO i=2-Olx,sNx+Olx-1
255                vrtDx(i,j) = (vort3(i+1,j)-vort3(i,j))
256         &                  *recip_DXG(i,j,bi,bj)
257         &                  *maskS(i,j,k,bi,bj)
258              ENDDO
259             ENDDO
260    C-       gradient in y direction:
261             DO j=2-Oly,sNy+Oly-1
262              DO i=2-Olx,sNx+Olx
263                vrtDy(i,j) = (vort3(i,j+1)-vort3(i,j))
264         &                  *recip_DYG(i,j,bi,bj)
265         &                  *maskW(i,j,k,bi,bj)
266              ENDDO
267             ENDDO
268    
269         ENDIF         ENDIF
270    
271         DO j=2-Oly,sNy+Oly-1         DO j=2-Oly,sNy+Oly-1
# Line 281  C Velocity Reynolds Scale Line 307  C Velocity Reynolds Scale
307    
308           IF (useFullLeith.and.calcleith) THEN           IF (useFullLeith.and.calcleith) THEN
309  C This is the vector magnitude of the vorticity gradient squared  C This is the vector magnitude of the vorticity gradient squared
310            grdVrt=0.25 _d 0*(            grdVrt=0.25 _d 0*( (vrtDx(i,j+1)*vrtDx(i,j+1)
311       &     ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2       &                        + vrtDx(i,j)*vrtDx(i,j) )
312       &     +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2       &                     + (vrtDy(i+1,j)*vrtDy(i+1,j)
313       &     +((vort3(i+1,j+1)-vort3(i,j+1))       &                        + vrtDy(i,j)*vrtDy(i,j) )  )
      &               *recip_DXG(i,j+1,bi,bj))**2  
      &     +((vort3(i+1,j+1)-vort3(i+1,j))  
      &               *recip_DYG(i+1,j,bi,bj))**2)  
314    
315  C This is the vector magnitude of grad (div.v) squared  C This is the vector magnitude of grad (div.v) squared
316  C Using it in Leith serves to damp instabilities in w.  C Using it in Leith serves to damp instabilities in w.
# Line 307  C Using it in Leith serves to damp insta Line 330  C Using it in Leith serves to damp insta
330           ELSEIF (calcleith) THEN           ELSEIF (calcleith) THEN
331  C but this approximation will work on cube  C but this approximation will work on cube
332  c (and differs by as much as 4X)  c (and differs by as much as 4X)
333            grdVrt=abs((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))            grdVrt=max( abs(vrtDx(i,j+1)), abs(vrtDx(i,j)) )
334            grdVrt=max(grdVrt,            grdVrt=max( grdVrt, abs(vrtDy(i+1,j)) )
335       &     abs((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj)))            grdVrt=max( grdVrt, abs(vrtDy(i,j))   )
           grdVrt=max(grdVrt,  
      &     abs((vort3(i+1,j+1)-vort3(i,j+1))*recip_DXG(i,j+1,bi,bj)))  
           grdVrt=max(grdVrt,  
      &     abs((vort3(i+1,j+1)-vort3(i+1,j))*recip_DYG(i+1,j,bi,bj)))  
336    
337    c This approximation is good to the same order as above...
338            grdDiv=max( abs(divDx(i+1,j)), abs(divDx(i,j)) )            grdDiv=max( abs(divDx(i+1,j)), abs(divDx(i,j)) )
339            grdDiv=max( grdDiv, abs(divDy(i,j+1)) )            grdDiv=max( grdDiv, abs(divDy(i,j+1)) )
340            grdDiv=max( grdDiv, abs(divDy(i,j))   )            grdDiv=max( grdDiv, abs(divDy(i,j))   )
341    
 c This approximation is good to the same order as above...  
342            viscAh_Dlth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3            viscAh_Dlth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3
343            viscA4_Dlth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5            viscA4_Dlth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5
344            viscAh_DlthD(i,j)=((leithD2fac*grdDiv))*L3            viscAh_DlthD(i,j)=((leithD2fac*grdDiv))*L3
# Line 398  C Velocity Reynolds Scale (Pb here at CS Line 417  C Velocity Reynolds Scale (Pb here at CS
417    
418  C This is the vector magnitude of the vorticity gradient squared  C This is the vector magnitude of the vorticity gradient squared
419           IF (useFullLeith.and.calcleith) THEN           IF (useFullLeith.and.calcleith) THEN
420            grdVrt=0.25 _d 0*(            grdVrt=0.25 _d 0*( (vrtDx(i-1,j)*vrtDx(i-1,j)
421       &     ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2       &                        + vrtDx(i,j)*vrtDx(i,j) )
422       &     +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2       &                     + (vrtDy(i,j-1)*vrtDy(i,j-1)
423       &     +((vort3(i-1,j)-vort3(i,j))*recip_DXG(i-1,j,bi,bj))**2       &                        + vrtDy(i,j)*vrtDy(i,j) )  )
      &     +((vort3(i,j-1)-vort3(i,j))*recip_DYG(i,j-1,bi,bj))**2)  
424    
425  C This is the vector magnitude of grad(div.v) squared  C This is the vector magnitude of grad(div.v) squared
426            grdDiv=0.25 _d 0*( (divDx(i,j-1)*divDx(i,j-1)            grdDiv=0.25 _d 0*( (divDx(i,j-1)*divDx(i,j-1)
# Line 421  C This is the vector magnitude of grad(d Line 439  C This is the vector magnitude of grad(d
439    
440           ELSEIF (calcleith) THEN           ELSEIF (calcleith) THEN
441  C but this approximation will work on cube (and differs by 4X)  C but this approximation will work on cube (and differs by 4X)
442            grdVrt=abs((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))            grdVrt=max( abs(vrtDx(i-1,j)), abs(vrtDx(i,j)) )
443            grdVrt=max(grdVrt,            grdVrt=max( grdVrt, abs(vrtDy(i,j-1)) )
444       &     abs((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj)))            grdVrt=max( grdVrt, abs(vrtDy(i,j))   )
           grdVrt=max(grdVrt,  
      &     abs((vort3(i-1,j)-vort3(i,j))*recip_DXG(i-1,j,bi,bj)))  
           grdVrt=max(grdVrt,  
      &     abs((vort3(i,j-1)-vort3(i,j))*recip_DYG(i,j-1,bi,bj)))  
445    
446            grdDiv=max( abs(divDx(i,j)), abs(divDx(i,j-1)) )            grdDiv=max( abs(divDx(i,j)), abs(divDx(i,j-1)) )
447            grdDiv=max( grdDiv, abs(divDy(i,j))   )            grdDiv=max( grdDiv, abs(divDy(i,j))   )
Legend:
Removed from v.1.18  
changed lines
  Added in v.1.20
  ViewVC Help
Powered by ViewVC 1.1.22