/[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.10 by jmc, Thu Sep 22 00:21:23 2005 UTC revision 1.20 by jmc, Wed Oct 12 20:24:22 2005 UTC
# Line 8  C $Name$ Line 8  C $Name$
8       I        bi,bj,k,       I        bi,bj,k,
9       O        viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,       O        viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,
10       O        harmonic,biharmonic,useVariableViscosity,       O        harmonic,biharmonic,useVariableViscosity,
11       I        hDiv,vort3,tension,strain,KE,hfacZ,       I        hDiv,vort3,tension,strain,KE,hFacZ,
12       I        myThid)       I        myThid)
13    
14        IMPLICIT NONE        IMPLICIT NONE
# Line 17  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 51  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 91  C     == Routine arguments == Line 91  C     == Routine arguments ==
91  C     == Local variables ==  C     == Local variables ==
92        INTEGER I,J        INTEGER I,J
93        _RL smag2fac, smag4fac        _RL smag2fac, smag4fac
94          _RL leith2fac, leith4fac
95          _RL leithD2fac, leithD4fac
96        _RL viscAhRe_max, viscA4Re_max        _RL viscAhRe_max, viscA4Re_max
97        _RL Alin,Alinmin,grdVrt,grdDiv        _RL Alin,grdVrt,grdDiv, keZpt
98        _RL recip_dt,L2,L3,L4,L5,L2rdt,L4rdt        _RL recip_dt,L2,L3,L4,L5,L2rdt,L4rdt
99        _RL Uscl,U4scl        _RL Uscl,U4scl
100          _RL divDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
101          _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 181  C     == Local variables == Line 187  C     == Local variables ==
187          smag4fac=0. _d 0          smag4fac=0. _d 0
188        ENDIF        ENDIF
189    
190          IF (calcleith) THEN
191            IF (useFullLeith) THEN
192             leith2fac =(viscC2leith /pi)**6
193             leithD2fac=(viscC2leithD/pi)**6
194             leith4fac =0.015625 _d 0*(viscC4leith /pi)**6
195             leithD4fac=0.015625 _d 0*(viscC4leithD/pi)**6
196            ELSE
197             leith2fac =(viscC2leith /pi)**3
198             leithD2fac=(viscC2leithD/pi)**3
199             leith4fac =0.125 _d 0*(viscC4leith /pi)**3
200             leithD4fac=0.125 _d 0*(viscC4leithD/pi)**3
201            ENDIF
202          ELSE
203            leith2fac=0. _d 0
204            leith4fac=0. _d 0
205            leithD2fac=0. _d 0
206            leithD4fac=0. _d 0
207          ENDIF
208    
209  C     - Viscosity  C     - Viscosity
210        IF (useVariableViscosity) THEN        IF (useVariableViscosity) THEN
211    
212    C-     Initialise to zero gradient of vorticity & divergence:
213           DO j=1-Oly,sNy+Oly
214            DO i=1-Olx,sNx+Olx
215              divDx(i,j) = 0.
216              divDy(i,j) = 0.
217              vrtDx(i,j) = 0.
218              vrtDy(i,j) = 0.
219            ENDDO
220           ENDDO
221    
222           IF (calcleith) THEN
223    C      horizontal gradient of horizontal divergence:
224    
225    C-       gradient in x direction:
226    #ifndef ALLOW_AUTODIFF_TAMC
227             IF (useCubedSphereExchange) THEN
228    C        to compute d/dx(hDiv), fill corners with appropriate values:
229               CALL FILL_CS_CORNER_TR_RL( .TRUE., hDiv, bi,bj, myThid )
230             ENDIF
231    #endif
232             DO j=2-Oly,sNy+Oly-1
233              DO i=2-Olx,sNx+Olx-1
234                divDx(i,j) = (hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj)
235              ENDDO
236             ENDDO
237    
238    C-       gradient in y direction:
239    #ifndef ALLOW_AUTODIFF_TAMC
240             IF (useCubedSphereExchange) THEN
241    C        to compute d/dy(hDiv), fill corners with appropriate values:
242               CALL FILL_CS_CORNER_TR_RL(.FALSE., hDiv, bi,bj, myThid )
243             ENDIF
244    #endif
245             DO j=2-Oly,sNy+Oly-1
246              DO i=2-Olx,sNx+Olx-1
247                divDy(i,j) = (hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj)
248              ENDDO
249             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
270    
271         DO j=2-Oly,sNy+Oly-1         DO j=2-Oly,sNy+Oly-1
272          DO i=2-Olx,sNx+Olx-1          DO i=2-Olx,sNx+Olx-1
273  CCCCCCCCCCCCCCC Divergence Point CalculationsCCCCCCCCCCCCCCCCCCCC  CCCCCCCCCCCCCCC Divergence Point CalculationsCCCCCCCCCCCCCCCCCCCC
274    
275  C These are (powers of) length scales  C These are (powers of) length scales
276           L2=2. _d 0/((recip_DXF(I,J,bi,bj)**2+recip_DYF(I,J,bi,bj)**2))           IF (useAreaViscLength) THEN
277              L2=rA(i,j,bi,bj)
278             ELSE
279              L2=2. _d 0/((recip_DXF(I,J,bi,bj)**2+recip_DYF(I,J,bi,bj)**2))
280             ENDIF
281           L3=(L2**1.5)           L3=(L2**1.5)
282           L4=(L2**2)           L4=(L2**2)
283           L5=(L2**2.5)           L5=(L2**2.5)
284    
285           L2rdt=0.25 _d 0*recip_dt*L2           L2rdt=0.25 _d 0*recip_dt*L2
286    
287           L4rdt=recip_dt/( 6. _d 0*(recip_DXF(I,J,bi,bj)**4           IF (useAreaViscLength) THEN
288              L4rdt=0.125 _d 0*recip_dt*rA(i,j,bi,bj)**2
289             ELSE
290              L4rdt=recip_dt/( 6. _d 0*(recip_DXF(I,J,bi,bj)**4
291       &                            +recip_DYF(I,J,bi,bj)**4)       &                            +recip_DYF(I,J,bi,bj)**4)
292       &                   +8. _d 0*((recip_DXF(I,J,bi,bj)       &                   +8. _d 0*((recip_DXF(I,J,bi,bj)
293       &                             *recip_DYF(I,J,bi,bj))**2) )       &                             *recip_DYF(I,J,bi,bj))**2) )
294             ENDIF
295    
296  C Velocity Reynolds Scale  C Velocity Reynolds Scale
297           Uscl=sqrt(KE(i,j)*L2)*viscAhRe_max           IF ( viscAhRe_max.GT.0. .AND. KE(i,j).GT.0. ) THEN
298           U4scl=sqrt(KE(i,j))*L3*viscA4Re_max             Uscl=sqrt(KE(i,j)*L2)*viscAhRe_max
299             ELSE
300               Uscl=0.
301             ENDIF
302             IF ( viscA4Re_max.GT.0. .AND. KE(i,j).GT.0. ) THEN
303               U4scl=sqrt(KE(i,j))*L3*viscA4Re_max
304             ELSE
305               U4scl=0.
306             ENDIF
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.
317            grdDiv=0.25 _d 0*(            grdDiv=0.25 _d 0*( (divDx(i+1,j)*divDx(i+1,j)
318       &     ((hDiv(i+1,j)-hDiv(i,j))*recip_DXC(i+1,j,bi,bj))**2       &                        + divDx(i,j)*divDx(i,j) )
319       &     +((hDiv(i,j+1)-hDiv(i,j))*recip_DYC(i,j+1,bi,bj))**2       &                     + (divDy(i,j+1)*divDy(i,j+1)
320       &     +((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)  
321    
322            viscAh_DLth(i,j)=            viscAh_DLth(i,j)=
323       &     sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3       &     sqrt(leith2fac*grdVrt+leithD2fac*grdDiv)*L3
324            viscA4_DLth(i,j)=0.125 _d 0*            viscA4_DLth(i,j)=
325       &     sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5       &     sqrt(leith4fac*grdVrt+leithD4fac*grdDiv)*L5
326            viscAh_DLthd(i,j)=            viscAh_DLthd(i,j)=
327       &     sqrt(viscC2leithD**2*grdDiv)*L3       &     sqrt(leithD2fac*grdDiv)*L3
328            viscA4_DLthd(i,j)=0.125 _d 0*            viscA4_DLthd(i,j)=
329       &     sqrt(viscC4leithD**2*grdDiv)*L5       &     sqrt(leithD4fac*grdDiv)*L5
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)))  
   
           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)))  
336    
337  c This approximation is good to the same order as above...  c This approximation is good to the same order as above...
338            viscAh_Dlth(i,j)=            grdDiv=max( abs(divDx(i+1,j)), abs(divDx(i,j)) )
339       &      (viscC2leith*grdVrt+(viscC2leithD*grdDiv))*L3            grdDiv=max( grdDiv, abs(divDy(i,j+1)) )
340            viscA4_Dlth(i,j)=0.125 _d 0*            grdDiv=max( grdDiv, abs(divDy(i,j))   )
341       &      (viscC4leith*grdVrt+(viscC4leithD*grdDiv))*L5  
342            viscAh_DlthD(i,j)=            viscAh_Dlth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3
343       &      ((viscC2leithD*grdDiv))*L3            viscA4_Dlth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5
344            viscA4_DlthD(i,j)=0.125 _d 0*            viscAh_DlthD(i,j)=((leithD2fac*grdDiv))*L3
345       &      ((viscC4leithD*grdDiv))*L5            viscA4_DlthD(i,j)=((leithD4fac*grdDiv))*L5
346           ELSE           ELSE
347            viscAh_Dlth(i,j)=0. _d 0            viscAh_Dlth(i,j)=0. _d 0
348            viscA4_Dlth(i,j)=0. _d 0            viscA4_Dlth(i,j)=0. _d 0
# Line 295  C  BiHarmonic on Div.u points Line 380  C  BiHarmonic on Div.u points
380    
381  CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC  CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC
382  C These are (powers of) length scales  C These are (powers of) length scales
383           L2=2. _d 0/((recip_DXV(I,J,bi,bj)**2+recip_DYU(I,J,bi,bj)**2))           IF (useAreaViscLength) THEN
384              L2=rAz(i,j,bi,bj)
385             ELSE
386              L2=2. _d 0/((recip_DXV(I,J,bi,bj)**2+recip_DYU(I,J,bi,bj)**2))
387             ENDIF
388    
389           L3=(L2**1.5)           L3=(L2**1.5)
390           L4=(L2**2)           L4=(L2**2)
391           L5=(L2**2.5)           L5=(L2**2.5)
392    
393           L2rdt=0.25 _d 0*recip_dt*L2           L2rdt=0.25 _d 0*recip_dt*L2
394           L4rdt=recip_dt/           IF (useAreaViscLength) THEN
395       &     ( 6. _d 0*(recip_DXF(I,J,bi,bj)**4+recip_DYF(I,J,bi,bj)**4)            L4rdt=0.125 _d 0*recip_dt*rAz(i,j,bi,bj)**2
396       &      +8. _d 0*((recip_DXF(I,J,bi,bj)*recip_DYF(I,J,bi,bj))**2))           ELSE
397              L4rdt=recip_dt/
398         &     ( 6. _d 0*(recip_DXV(I,J,bi,bj)**4+recip_DYU(I,J,bi,bj)**4)
399         &      +8. _d 0*((recip_DXV(I,J,bi,bj)*recip_DYU(I,J,bi,bj))**2))
400             ENDIF
401    
402  C Velocity Reynolds Scale  C Velocity Reynolds Scale (Pb here at CS-grid corners !)
403           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
404       &         *L2)*viscAhRe_max             keZpt=0.25 _d 0*( (KE(i,j)+KE(i-1,j-1))
405           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)) )
406       &         *L3*viscA4Re_max             IF ( keZpt.GT.0. ) THEN
407                 Uscl = sqrt(keZpt*L2)*viscAhRe_max
408                 U4scl= sqrt(keZpt)*L3*viscA4Re_max
409               ELSE
410                 Uscl =0.
411                 U4scl=0.
412               ENDIF
413             ELSE
414               Uscl =0.
415               U4scl=0.
416             ENDIF
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*(            grdDiv=0.25 _d 0*( (divDx(i,j-1)*divDx(i,j-1)
427       &     ((hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj))**2       &                        + divDx(i,j)*divDx(i,j) )
428       &     +((hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj))**2       &                     + (divDy(i-1,j)*divDy(i-1,j)
429       &     +((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)  
430    
431            viscAh_ZLth(i,j)=            viscAh_ZLth(i,j)=
432       &     sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3       &     sqrt(leith2fac*grdVrt+leithD2fac*grdDiv)*L3
433            viscA4_ZLth(i,j)=0.125 _d 0*            viscA4_ZLth(i,j)=
434       &     sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5       &     sqrt(leith4fac*grdVrt+leithD4fac*grdDiv)*L5
435            viscAh_ZLthD(i,j)=            viscAh_ZLthD(i,j)=
436       &     sqrt(viscC2leithD**2*grdDiv)*L3       &     sqrt(leithD2fac*grdDiv)*L3
437            viscA4_ZLthD(i,j)=0.125 _d 0*            viscA4_ZLthD(i,j)=
438       &     sqrt(viscC4leithD**2*grdDiv)*L5       &     sqrt(leithD4fac*grdDiv)*L5
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))   )
445            grdVrt=max(grdVrt,  
446       &     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)) )
447            grdVrt=max(grdVrt,            grdDiv=max( grdDiv, abs(divDy(i,j))   )
448       &     abs((vort3(i,j-1)-vort3(i,j))*recip_DYG(i,j-1,bi,bj)))            grdDiv=max( grdDiv, abs(divDy(i-1,j)) )
449    
450            grdDiv=abs((hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj))            viscAh_ZLth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3
451            grdDiv=max(grdDiv,            viscA4_ZLth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5
452       &     abs((hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj)))            viscAh_ZLthD(i,j)=(leithD2fac*grdDiv)*L3
453            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  
454           ELSE           ELSE
455            viscAh_ZLth(i,j)=0. _d 0            viscAh_ZLth(i,j)=0. _d 0
456            viscA4_ZLth(i,j)=0. _d 0            viscA4_ZLth(i,j)=0. _d 0
Legend:
Removed from v.1.10  
changed lines
  Added in v.1.20
  ViewVC Help
Powered by ViewVC 1.1.22