/[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.13 by baylor, Mon Sep 26 15:27:11 2005 UTC revision 1.23 by baylor, Fri Jul 7 18:52:10 2006 UTC
# Line 1  Line 1 
1    C $Header$
2    C $Name$
3    
4  #include "MOM_COMMON_OPTIONS.h"  #include "MOM_COMMON_OPTIONS.h"
5    
# Line 15  C     Calculate horizontal viscosities ( Line 17  C     Calculate horizontal viscosities (
17  C     harmonic viscosity=  C     harmonic viscosity=
18  C       viscAh (or viscAhD on div pts and viscAhZ on zeta pts)  C       viscAh (or viscAhD on div pts and viscAhZ on zeta pts)
19  C       +0.25*L**2*viscAhGrid/deltaT  C       +0.25*L**2*viscAhGrid/deltaT
20  C       +sqrt(viscC2leith**2*grad(Vort3)**2  C       +sqrt((viscC2leith/pi)**6*grad(Vort3)**2
21  C             +viscC2leithD**2*grad(hDiv)**2)*L**3  C             +(viscC2leithD/pi)**6*grad(hDiv)**2)*L**3
22  C       +(viscC2smag/pi)**2*L**2*sqrt(Tension**2+Strain**2)  C       +(viscC2smag/pi)**2*L**2*sqrt(Tension**2+Strain**2)
23  C  C
24  C     biharmonic viscosity=  C     biharmonic viscosity=
25  C       viscA4 (or viscA4D on div pts and viscA4Z on zeta pts)  C       viscA4 (or viscA4D on div pts and viscA4Z on zeta pts)
26  C       +0.25*0.125*L**4*viscA4Grid/deltaT (approx)  C       +0.25*0.125*L**4*viscA4Grid/deltaT (approx)
27  C       +0.125*L**5*sqrt(viscC4leith**2*grad(Vort3)**2  C       +0.125*L**5*sqrt((viscC4leith/pi)**6*grad(Vort3)**2
28  C                        +viscC4leithD**2*grad(hDiv)**2)  C                        +(viscC4leithD/pi)**6*grad(hDiv)**2)
29  C       +0.125*L**4*(viscC4smag/pi)**2*sqrt(Tension**2+Strain**2)  C       +0.125*L**4*(viscC4smag/pi)**2*sqrt(Tension**2+Strain**2)
30  C  C
31  C     Note that often 0.125*L**2 is the scale between harmonic and  C     Note that often 0.125*L**2 is the scale between harmonic and
# Line 49  C      harmonic viscosity>0.25*viscAhgri Line 51  C      harmonic viscosity>0.25*viscAhgri
51  C      biharmonic viscosity>viscA4gridmax*L**4/32/deltaT (approx)  C      biharmonic viscosity>viscA4gridmax*L**4/32/deltaT (approx)
52  C  C
53  C     RECOMMENDED VALUES  C     RECOMMENDED VALUES
54  C     viscC2Leith=?  C     viscC2Leith=1-3
55  C     viscC2LeithD=?  C     viscC2LeithD=1-3
56  C     viscC4Leith=?  C     viscC4Leith=1-3
57  C     viscC4LeithD=?  C     viscC4LeithD=1.5-3
58  C     viscC2smag=2.2-4 (Griffies and Hallberg,2000)  C     viscC2smag=2.2-4 (Griffies and Hallberg,2000)
59  C               0.2-0.9 (Smagorinsky,1993)  C               0.2-0.9 (Smagorinsky,1993)
60  C     viscC4smag=2.2-4 (Griffies and Hallberg,2000)  C     viscC4smag=2.2-4 (Griffies and Hallberg,2000)
# Line 70  C     == Global variables == Line 72  C     == Global variables ==
72  #include "GRID.h"  #include "GRID.h"
73  #include "EEPARAMS.h"  #include "EEPARAMS.h"
74  #include "PARAMS.h"  #include "PARAMS.h"
75  #ifdef ALLOW_EXCH2  #ifdef ALLOW_NONHYDROSTATIC
76  #include "W2_EXCH2_TOPOLOGY.h"  #include "NH_VARS.h"
77  #include "W2_EXCH2_PARAMS.h"  #endif
 #endif /* ALLOW_EXCH2 */  
78    
79  C     == Routine arguments ==  C     == Routine arguments ==
80        INTEGER bi,bj,k        INTEGER bi,bj,k
# Line 92  C     == Routine arguments == Line 93  C     == Routine arguments ==
93    
94  C     == Local variables ==  C     == Local variables ==
95        INTEGER I,J        INTEGER I,J
96          INTEGER kp1
97        _RL smag2fac, smag4fac        _RL smag2fac, smag4fac
98          _RL leith2fac, leith4fac
99          _RL leithD2fac, leithD4fac
100        _RL viscAhRe_max, viscA4Re_max        _RL viscAhRe_max, viscA4Re_max
101        _RL Alin,grdVrt,grdDiv        _RL Alin,grdVrt,grdDiv, keZpt
102        _RL recip_dt,L2,L3,L4,L5,L2rdt,L4rdt        _RL recip_dt,L2,L3,L4,L5,L2rdt,L4rdt
103        _RL Uscl,U4scl        _RL Uscl,U4scl
104          _RL divDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105          _RL divDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
106          _RL vrtDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
107          _RL vrtDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
108        _RL viscAh_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscAh_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
109        _RL viscAh_DMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscAh_DMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
110        _RL viscA4_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscA4_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 118  C     == Local variables == Line 126  C     == Local variables ==
126        _RL viscA4_ZSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscA4_ZSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
127        _RL viscA4_DSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscA4_DSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
128        LOGICAL calcLeith,calcSmag        LOGICAL calcLeith,calcSmag
       LOGICAL northWestCorner, northEastCorner,  
      &        southWestCorner, southEastCorner  
 #ifdef ALLOW_EXCH2  
       INTEGER myTile  
 #endif /* ALLOW_EXCH2 */  
   
 C     Special stuff for Cubed Sphere  
       southWestCorner = .FALSE.  
       southEastCorner = .FALSE.  
       northWestCorner = .FALSE.  
       northEastCorner = .FALSE.  
       IF (useCubedSphereExchange) THEN  
 #ifdef ALLOW_EXCH2  
        myTile = W2_myTileList(bi)  
        IF ( exch2_isWedge(myTile) .EQ. 1 .AND.  
      &      exch2_isSedge(myTile) .EQ. 1 ) THEN  
         southWestCorner = .TRUE.  
        ENDIF  
        IF ( exch2_isEedge(myTile) .EQ. 1 .AND.  
      &      exch2_isSedge(myTile) .EQ. 1 ) THEN  
         southEastCorner = .TRUE.  
        ENDIF  
        IF ( exch2_isEedge(myTile) .EQ. 1 .AND.  
      &      exch2_isNedge(myTile) .EQ. 1 ) THEN  
         northEastCorner = .TRUE.  
        ENDIF  
        IF ( exch2_isWedge(myTile) .EQ. 1 .AND.  
      &      exch2_isNedge(myTile) .EQ. 1 ) THEN  
         northWestCorner = .TRUE.  
        ENDIF  
 #else  
        southWestCorner = .TRUE.  
        southEastCorner = .TRUE.  
        northWestCorner = .TRUE.  
        northEastCorner = .TRUE.  
 #endif /* ALLOW_EXCH2 */  
       ENDIF  
129    
130        useVariableViscosity=        useVariableViscosity=
131       &      (viscAhGrid.NE.0.)       &      (viscAhGrid.NE.0.)
# Line 220  C     Special stuff for Cubed Sphere Line 191  C     Special stuff for Cubed Sphere
191          smag4fac=0. _d 0          smag4fac=0. _d 0
192        ENDIF        ENDIF
193    
194          IF (calcleith) THEN
195            IF (useFullLeith) THEN
196             leith2fac =(viscC2leith /pi)**6
197             leithD2fac=(viscC2leithD/pi)**6
198             leith4fac =0.015625 _d 0*(viscC4leith /pi)**6
199             leithD4fac=0.015625 _d 0*(viscC4leithD/pi)**6
200            ELSE
201             leith2fac =(viscC2leith /pi)**3
202             leithD2fac=(viscC2leithD/pi)**3
203             leith4fac =0.125 _d 0*(viscC4leith /pi)**3
204             leithD4fac=0.125 _d 0*(viscC4leithD/pi)**3
205            ENDIF
206          ELSE
207            leith2fac=0. _d 0
208            leith4fac=0. _d 0
209            leithD2fac=0. _d 0
210            leithD4fac=0. _d 0
211          ENDIF
212    
213    #ifdef ALLOW_AUTODIFF_TAMC
214          IF ( calcLeith .OR. calcSmag ) THEN
215           STOP 'calcLeith or calcSmag not implemented for ADJOINT'
216          ENDIF
217          DO j=1-Oly,sNy+Oly
218            DO i=1-Olx,sNx+Olx
219              viscAh_D(i,j)=viscAhD
220              viscAh_Z(i,j)=viscAhZ
221              viscA4_D(i,j)=viscA4D
222              viscA4_Z(i,j)=viscA4Z
223    c
224              visca4_zsmg(i,j) = 0. _d 0
225              viscah_zsmg(i,j) = 0. _d 0
226    c
227              viscAh_Dlth(i,j) = 0. _d 0
228              viscA4_Dlth(i,j) = 0. _d 0
229              viscAh_DlthD(i,j)= 0. _d 0
230              viscA4_DlthD(i,j)= 0. _d 0
231    c
232              viscAh_DSmg(i,j) = 0. _d 0
233              viscA4_DSmg(i,j) = 0. _d 0
234    c
235              viscAh_ZLth(i,j) = 0. _d 0
236              viscA4_ZLth(i,j) = 0. _d 0
237              viscAh_ZLthD(i,j)= 0. _d 0
238              viscA4_ZLthD(i,j)= 0. _d 0
239            ENDDO
240          ENDDO
241    #endif
242    
243    
244    
245  C     - Viscosity  C     - Viscosity
246        IF (useVariableViscosity) THEN        IF (useVariableViscosity) THEN
247    
248    C-     Initialise to zero gradient of vorticity & divergence:
249           DO j=1-Oly,sNy+Oly
250            DO i=1-Olx,sNx+Olx
251              divDx(i,j) = 0.
252              divDy(i,j) = 0.
253              vrtDx(i,j) = 0.
254              vrtDy(i,j) = 0.
255            ENDDO
256           ENDDO
257    
258           IF (calcleith) THEN
259    C      horizontal gradient of horizontal divergence:
260    
261    C-       gradient in x direction:
262    #ifndef ALLOW_AUTODIFF_TAMC
263             IF (useCubedSphereExchange) THEN
264    C        to compute d/dx(hDiv), fill corners with appropriate values:
265               CALL FILL_CS_CORNER_TR_RL( .TRUE., hDiv, bi,bj, myThid )
266             ENDIF
267    #endif
268             DO j=2-Oly,sNy+Oly-1
269              DO i=2-Olx,sNx+Olx-1
270                divDx(i,j) = (hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj)
271              ENDDO
272             ENDDO
273    
274    C-       gradient in y direction:
275    #ifndef ALLOW_AUTODIFF_TAMC
276             IF (useCubedSphereExchange) THEN
277    C        to compute d/dy(hDiv), fill corners with appropriate values:
278               CALL FILL_CS_CORNER_TR_RL(.FALSE., hDiv, bi,bj, myThid )
279             ENDIF
280    #endif
281             DO j=2-Oly,sNy+Oly-1
282              DO i=2-Olx,sNx+Olx-1
283                divDy(i,j) = (hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj)
284              ENDDO
285             ENDDO
286    
287    C      horizontal gradient of vertical vorticity:
288    C-       gradient in x direction:
289             DO j=2-Oly,sNy+Oly
290              DO i=2-Olx,sNx+Olx-1
291                vrtDx(i,j) = (vort3(i+1,j)-vort3(i,j))
292         &                  *recip_DXG(i,j,bi,bj)
293         &                  *maskS(i,j,k,bi,bj)
294              ENDDO
295             ENDDO
296    C-       gradient in y direction:
297             DO j=2-Oly,sNy+Oly-1
298              DO i=2-Olx,sNx+Olx
299                vrtDy(i,j) = (vort3(i,j+1)-vort3(i,j))
300         &                  *recip_DYG(i,j,bi,bj)
301         &                  *maskW(i,j,k,bi,bj)
302              ENDDO
303             ENDDO
304    
305           ENDIF
306    
307         DO j=2-Oly,sNy+Oly-1         DO j=2-Oly,sNy+Oly-1
308          DO i=2-Olx,sNx+Olx-1          DO i=2-Olx,sNx+Olx-1
309  CCCCCCCCCCCCCCC Divergence Point CalculationsCCCCCCCCCCCCCCCCCCCC  CCCCCCCCCCCCCCC Divergence Point CalculationsCCCCCCCCCCCCCCCCCCCC
# Line 248  C These are (powers of) length scales Line 330  C These are (powers of) length scales
330           ENDIF           ENDIF
331    
332  C Velocity Reynolds Scale  C Velocity Reynolds Scale
333           Uscl=sqrt(KE(i,j)*L2)*viscAhRe_max           IF ( viscAhRe_max.GT.0. .AND. KE(i,j).GT.0. ) THEN
334           U4scl=sqrt(KE(i,j))*L3*viscA4Re_max             Uscl=sqrt(KE(i,j)*L2)*viscAhRe_max
335             ELSE
336               Uscl=0.
337             ENDIF
338             IF ( viscA4Re_max.GT.0. .AND. KE(i,j).GT.0. ) THEN
339               U4scl=sqrt(KE(i,j))*L3*viscA4Re_max
340             ELSE
341               U4scl=0.
342             ENDIF
343    
344    #ifndef ALLOW_AUTODIFF_TAMC
345           IF (useFullLeith.and.calcleith) THEN           IF (useFullLeith.and.calcleith) THEN
346  C This is the vector magnitude of the vorticity gradient squared  C This is the vector magnitude of the vorticity gradient squared
347            grdVrt=0.25 _d 0*(            grdVrt=0.25 _d 0*( (vrtDx(i,j+1)*vrtDx(i,j+1)
348       &     ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2       &                        + vrtDx(i,j)*vrtDx(i,j) )
349       &     +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2       &                     + (vrtDy(i+1,j)*vrtDy(i+1,j)
350       &     +((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)  
351    
352  C This is the vector magnitude of grad (div.v) squared  C This is the vector magnitude of grad (div.v) squared
353  C Using it in Leith serves to damp instabilities in w.  C Using it in Leith serves to damp instabilities in w.
354            grdDiv=0.25 _d 0*(            grdDiv=0.25 _d 0*( (divDx(i+1,j)*divDx(i+1,j)
355       &     ((hDiv(i+1,j)-hDiv(i,j))*recip_DXC(i+1,j,bi,bj))**2       &                        + divDx(i,j)*divDx(i,j) )
356       &     +((hDiv(i,j+1)-hDiv(i,j))*recip_DYC(i,j+1,bi,bj))**2       &                     + (divDy(i,j+1)*divDy(i,j+1)
357       &     +((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)  
358    
359            viscAh_DLth(i,j)=            viscAh_DLth(i,j)=
360       &     sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3       &     sqrt(leith2fac*grdVrt+leithD2fac*grdDiv)*L3
361            viscA4_DLth(i,j)=0.125 _d 0*            viscA4_DLth(i,j)=
362       &     sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5       &     sqrt(leith4fac*grdVrt+leithD4fac*grdDiv)*L5
363            viscAh_DLthd(i,j)=            viscAh_DLthd(i,j)=
364       &     sqrt(viscC2leithD**2*grdDiv)*L3       &     sqrt(leithD2fac*grdDiv)*L3
365            viscA4_DLthd(i,j)=0.125 _d 0*            viscA4_DLthd(i,j)=
366       &     sqrt(viscC4leithD**2*grdDiv)*L5       &     sqrt(leithD4fac*grdDiv)*L5
367           ELSEIF (calcleith) THEN           ELSEIF (calcleith) THEN
368  C but this approximation will work on cube  C but this approximation will work on cube
369  c (and differs by as much as 4X)  c (and differs by as much as 4X)
370            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)) )
371            grdVrt=max(grdVrt,            grdVrt=max( grdVrt, abs(vrtDy(i+1,j)) )
372       &     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)))  
   
           grdDiv=abs((hDiv(i+1,j)-hDiv(i,j))*recip_DXC(i+1,j,bi,bj))  
           grdDiv=max(grdDiv,  
      &     abs((hDiv(i,j+1)-hDiv(i,j))*recip_DYC(i,j+1,bi,bj)))  
           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)))  
373    
374  c This approximation is good to the same order as above...  c This approximation is good to the same order as above...
375            viscAh_Dlth(i,j)=            grdDiv=max( abs(divDx(i+1,j)), abs(divDx(i,j)) )
376       &      (viscC2leith*grdVrt+(viscC2leithD*grdDiv))*L3            grdDiv=max( grdDiv, abs(divDy(i,j+1)) )
377            viscA4_Dlth(i,j)=0.125 _d 0*            grdDiv=max( grdDiv, abs(divDy(i,j))   )
378       &      (viscC4leith*grdVrt+(viscC4leithD*grdDiv))*L5  
379            viscAh_DlthD(i,j)=            viscAh_Dlth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3
380       &      ((viscC2leithD*grdDiv))*L3            viscA4_Dlth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5
381            viscA4_DlthD(i,j)=0.125 _d 0*            viscAh_DlthD(i,j)=((leithD2fac*grdDiv))*L3
382       &      ((viscC4leithD*grdDiv))*L5            viscA4_DlthD(i,j)=((leithD4fac*grdDiv))*L5
383           ELSE           ELSE
384            viscAh_Dlth(i,j)=0. _d 0            viscAh_Dlth(i,j)=0. _d 0
385            viscA4_Dlth(i,j)=0. _d 0            viscA4_Dlth(i,j)=0. _d 0
# Line 323  c This approximation is good to the same Line 398  c This approximation is good to the same
398            viscAh_DSmg(i,j)=0. _d 0            viscAh_DSmg(i,j)=0. _d 0
399            viscA4_DSmg(i,j)=0. _d 0            viscA4_DSmg(i,j)=0. _d 0
400           ENDIF           ENDIF
401    #endif /* ALLOW_AUTODIFF_TAMC */
402    
403  C  Harmonic on Div.u points  C  Harmonic on Div.u points
404           Alin=viscAhD+viscAhGrid*L2rdt           Alin=viscAhD+viscAhGrid*L2rdt
# Line 340  C  BiHarmonic on Div.u points Line 416  C  BiHarmonic on Div.u points
416           viscA4_DMax(i,j)=min(viscA4GridMax*L4rdt,viscA4Max)           viscA4_DMax(i,j)=min(viscA4GridMax*L4rdt,viscA4Max)
417           viscA4_D(i,j)=min(viscA4_DMax(i,j),viscA4_D(i,j))           viscA4_D(i,j)=min(viscA4_DMax(i,j),viscA4_D(i,j))
418    
419    #ifdef ALLOW_NONHYDROSTATIC
420    C     /* Pass Viscosities to calc_gw, if constant, not necessary */
421    
422             kp1  = MIN(k+1,Nr)
423    
424             if (k .eq. 1) then
425              viscAh_W(i,j,kp1,bi,bj)=0.5*viscAh_D(i,j)
426              viscA4_W(i,j,kp1,bi,bj)=0.5*viscA4_D(i,j)
427    
428              viscAh_W(i,j,k,bi,bj)=viscAh_D(i,j) /* These values dont get used */
429              viscA4_W(i,j,k,bi,bj)=viscA4_D(i,j)
430             else
431    C Note that previous call of this function has already added half.
432              viscAh_W(i,j,kp1,bi,bj)=0.5*viscAh_D(i,j)
433              viscA4_W(i,j,kp1,bi,bj)=0.5*viscA4_D(i,j)
434    
435              viscAh_W(i,j,k,bi,bj)=viscAh_W(i,j,k,bi,bj)+0.5*viscAh_D(i,j)
436              viscA4_W(i,j,k,bi,bj)=viscA4_W(i,j,k,bi,bj)+0.5*viscA4_D(i,j)
437             endif
438    #endif /* ALLOW_NONHYDROSTATIC */
439    
440  CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC  CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC
441  C These are (powers of) length scales  C These are (powers of) length scales
442           IF (useAreaViscLength) THEN           IF (useAreaViscLength) THEN
# Line 354  C These are (powers of) length scales Line 451  C These are (powers of) length scales
451    
452           L2rdt=0.25 _d 0*recip_dt*L2           L2rdt=0.25 _d 0*recip_dt*L2
453           IF (useAreaViscLength) THEN           IF (useAreaViscLength) THEN
454            L4rdt=0.125 _d 0*recip_dt*RaZ(i,j,bi,bj)**2            L4rdt=0.125 _d 0*recip_dt*rAz(i,j,bi,bj)**2
455           ELSE           ELSE
456            L4rdt=recip_dt/            L4rdt=recip_dt/
457       &     ( 6. _d 0*(recip_DXV(I,J,bi,bj)**4+recip_DYU(I,J,bi,bj)**4)       &     ( 6. _d 0*(recip_DXV(I,J,bi,bj)**4+recip_DYU(I,J,bi,bj)**4)
458       &      +8. _d 0*((recip_DXV(I,J,bi,bj)*recip_DYU(I,J,bi,bj))**2))       &      +8. _d 0*((recip_DXV(I,J,bi,bj)*recip_DYU(I,J,bi,bj))**2))
459           ENDIF           ENDIF
460    
461  C Velocity Reynolds Scale  C Velocity Reynolds Scale (Pb here at CS-grid corners !)
462           Uscl=sqrt(0.25 _d 0*(KE(i,j)+KE(i,j+1)+KE(i+1,j)+KE(i+1,j+1))           IF ( viscAhRe_max.GT.0. .OR. viscA4Re_max.GT.0. ) THEN
463       &         *L2)*viscAhRe_max             keZpt=0.25 _d 0*( (KE(i,j)+KE(i-1,j-1))
464           U4scl=sqrt(0.25 _d 0*(KE(i,j)+KE(i,j+1)+KE(i+1,j)+KE(i+1,j+1)))       &                      +(KE(i-1,j)+KE(i,j-1)) )
465       &         *L3*viscA4Re_max             IF ( keZpt.GT.0. ) THEN
466                 Uscl = sqrt(keZpt*L2)*viscAhRe_max
467                 U4scl= sqrt(keZpt)*L3*viscA4Re_max
468               ELSE
469                 Uscl =0.
470                 U4scl=0.
471               ENDIF
472             ELSE
473               Uscl =0.
474               U4scl=0.
475             ENDIF
476    
477    #ifndef ALLOW_AUTODIFF_TAMC
478  C This is the vector magnitude of the vorticity gradient squared  C This is the vector magnitude of the vorticity gradient squared
479           IF (useFullLeith.and.calcleith) THEN           IF (useFullLeith.and.calcleith) THEN
480            grdVrt=0.25 _d 0*(            grdVrt=0.25 _d 0*( (vrtDx(i-1,j)*vrtDx(i-1,j)
481       &     ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2       &                        + vrtDx(i,j)*vrtDx(i,j) )
482       &     +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2       &                     + (vrtDy(i,j-1)*vrtDy(i,j-1)
483       &     +((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)  
484    
485  C This is the vector magnitude of grad(div.v) squared  C This is the vector magnitude of grad(div.v) squared
486            grdDiv=0.25 _d 0*(            grdDiv=0.25 _d 0*( (divDx(i,j-1)*divDx(i,j-1)
487       &     ((hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj))**2       &                        + divDx(i,j)*divDx(i,j) )
488       &     +((hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj))**2       &                     + (divDy(i-1,j)*divDy(i-1,j)
489       &     +((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)  
490    
491            viscAh_ZLth(i,j)=            viscAh_ZLth(i,j)=
492       &     sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3       &     sqrt(leith2fac*grdVrt+leithD2fac*grdDiv)*L3
493            viscA4_ZLth(i,j)=0.125 _d 0*            viscA4_ZLth(i,j)=
494       &     sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5       &     sqrt(leith4fac*grdVrt+leithD4fac*grdDiv)*L5
495            viscAh_ZLthD(i,j)=            viscAh_ZLthD(i,j)=
496       &     sqrt(viscC2leithD**2*grdDiv)*L3       &     sqrt(leithD2fac*grdDiv)*L3
497            viscA4_ZLthD(i,j)=0.125 _d 0*            viscA4_ZLthD(i,j)=
498       &     sqrt(viscC4leithD**2*grdDiv)*L5       &     sqrt(leithD4fac*grdDiv)*L5
499    
500           ELSEIF (calcleith) THEN           ELSEIF (calcleith) THEN
501  C but this approximation will work on cube (and differs by 4X)  C but this approximation will work on cube (and differs by 4X)
502            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)) )
503            grdVrt=max(grdVrt,            grdVrt=max( grdVrt, abs(vrtDy(i,j-1)) )
504       &     abs((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj)))            grdVrt=max( grdVrt, abs(vrtDy(i,j))   )
505            grdVrt=max(grdVrt,  
506       &     abs((vort3(i-1,j)-vort3(i,j))*recip_DXG(i-1,j,bi,bj)))            grdDiv=max( abs(divDx(i,j)), abs(divDx(i,j-1)) )
507            grdVrt=max(grdVrt,            grdDiv=max( grdDiv, abs(divDy(i,j))   )
508       &     abs((vort3(i,j-1)-vort3(i,j))*recip_DYG(i,j-1,bi,bj)))            grdDiv=max( grdDiv, abs(divDy(i-1,j)) )
509    
510            grdDiv=abs((hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj))            viscAh_ZLth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3
511            grdDiv=max(grdDiv,            viscA4_ZLth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5
512       &     abs((hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj)))            viscAh_ZLthD(i,j)=(leithD2fac*grdDiv)*L3
513            grdDiv=max(grdDiv,            viscA4_ZLthD(i,j)=(leithD4fac*grdDiv)*L5
      &     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)))  
   
           viscAh_ZLth(i,j)=(viscC2leith*grdVrt  
      &                     +(viscC2leithD*grdDiv))*L3  
           viscA4_ZLth(i,j)=0.125 _d 0*(viscC4leith*grdVrt  
      &                     +(viscC4leithD*grdDiv))*L5  
           viscAh_ZLthD(i,j)=((viscC2leithD*grdDiv))*L3  
           viscA4_ZLthD(i,j)=0.125 _d 0*((viscC4leithD*grdDiv))*L5  
514           ELSE           ELSE
515            viscAh_ZLth(i,j)=0. _d 0            viscAh_ZLth(i,j)=0. _d 0
516            viscA4_ZLth(i,j)=0. _d 0            viscA4_ZLth(i,j)=0. _d 0
# Line 430  C but this approximation will work on cu Line 526  C but this approximation will work on cu
526            viscA4_ZSmg(i,j)=smag4fac*L2*viscAh_ZSmg(i,j)            viscA4_ZSmg(i,j)=smag4fac*L2*viscAh_ZSmg(i,j)
527            viscAh_ZSmg(i,j)=smag2fac*viscAh_ZSmg(i,j)            viscAh_ZSmg(i,j)=smag2fac*viscAh_ZSmg(i,j)
528           ENDIF           ENDIF
529    #endif /* ALLOW_AUTODIFF_TAMC */
530    
531  C  Harmonic on Zeta points  C  Harmonic on Zeta points
532           Alin=viscAhZ+viscAhGrid*L2rdt           Alin=viscAhZ+viscAhGrid*L2rdt
# Line 465  C  BiHarmonic on Zeta points Line 562  C  BiHarmonic on Zeta points
562         CALL DIAGNOSTICS_FILL(viscA4_D,'VISCA4D ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscA4_D,'VISCA4D ',k,1,2,bi,bj,myThid)
563         CALL DIAGNOSTICS_FILL(viscAh_Z,'VISCAHZ ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscAh_Z,'VISCAHZ ',k,1,2,bi,bj,myThid)
564         CALL DIAGNOSTICS_FILL(viscA4_Z,'VISCA4Z ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscA4_Z,'VISCA4Z ',k,1,2,bi,bj,myThid)
565    #ifdef ALLOW_NONHYDROSTATIC
566           CALL DIAGNOSTICS_FILL(viscAh_W,'VISCAHW ',k,1,2,bi,bj,myThid)
567           CALL DIAGNOSTICS_FILL(viscA4_W,'VISCA4W ',k,1,2,bi,bj,myThid)
568    #endif
569    
570         CALL DIAGNOSTICS_FILL(viscAh_DMax,'VAHDMAX ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscAh_DMax,'VAHDMAX ',k,1,2,bi,bj,myThid)
571         CALL DIAGNOSTICS_FILL(viscA4_DMax,'VA4DMAX ',k,1,2,bi,bj,myThid)         CALL DIAGNOSTICS_FILL(viscA4_DMax,'VA4DMAX ',k,1,2,bi,bj,myThid)
Legend:
Removed from v.1.13  
changed lines
  Added in v.1.23
  ViewVC Help
Powered by ViewVC 1.1.22