/[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.11 by baylor, Thu Sep 22 14:00:59 2005 UTC revision 1.21 by heimbach, Thu Nov 24 00:06:37 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    #ifdef ALLOW_AUTODIFF_TAMC
210           DO j=1-Oly,sNy+Oly
211            DO i=1-Olx,sNx+Olx
212              viscAh_D(i,j)=viscAhD
213              viscAh_Z(i,j)=viscAhZ
214              viscA4_D(i,j)=viscA4D
215              viscA4_Z(i,j)=viscA4Z
216    c
217              visca4_zsmg(i,j) = 0. _d 0
218              viscah_zsmg(i,j) = 0. _d 0
219    c
220              viscAh_Dlth(i,j) = 0. _d 0
221              viscA4_Dlth(i,j) = 0. _d 0
222              viscAh_DlthD(i,j)= 0. _d 0
223              viscA4_DlthD(i,j)= 0. _d 0
224    c
225              viscAh_DSmg(i,j) = 0. _d 0
226              viscA4_DSmg(i,j) = 0. _d 0
227    c
228              viscAh_ZLth(i,j) = 0. _d 0
229              viscA4_ZLth(i,j) = 0. _d 0
230              viscAh_ZLthD(i,j)= 0. _d 0
231              viscA4_ZLthD(i,j)= 0. _d 0
232            ENDDO
233           ENDDO
234    #endif
235    
236    
237    
238  C     - Viscosity  C     - Viscosity
239        IF (useVariableViscosity) THEN        IF (useVariableViscosity) THEN
240    cph(
241    #ifndef ALLOW_AUTODIFF_TAMC
242    cph)
243    
244    C-     Initialise to zero gradient of vorticity & divergence:
245           DO j=1-Oly,sNy+Oly
246            DO i=1-Olx,sNx+Olx
247              divDx(i,j) = 0.
248              divDy(i,j) = 0.
249              vrtDx(i,j) = 0.
250              vrtDy(i,j) = 0.
251            ENDDO
252           ENDDO
253    
254           IF (calcleith) THEN
255    C      horizontal gradient of horizontal divergence:
256    
257    C-       gradient in x direction:
258    #ifndef ALLOW_AUTODIFF_TAMC
259             IF (useCubedSphereExchange) THEN
260    C        to compute d/dx(hDiv), fill corners with appropriate values:
261               CALL FILL_CS_CORNER_TR_RL( .TRUE., hDiv, bi,bj, myThid )
262             ENDIF
263    #endif
264             DO j=2-Oly,sNy+Oly-1
265              DO i=2-Olx,sNx+Olx-1
266                divDx(i,j) = (hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj)
267              ENDDO
268             ENDDO
269    
270    C-       gradient in y direction:
271    #ifndef ALLOW_AUTODIFF_TAMC
272             IF (useCubedSphereExchange) THEN
273    C        to compute d/dy(hDiv), fill corners with appropriate values:
274               CALL FILL_CS_CORNER_TR_RL(.FALSE., hDiv, bi,bj, myThid )
275             ENDIF
276    #endif
277             DO j=2-Oly,sNy+Oly-1
278              DO i=2-Olx,sNx+Olx-1
279                divDy(i,j) = (hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj)
280              ENDDO
281             ENDDO
282    
283    C      horizontal gradient of vertical vorticity:
284    C-       gradient in x direction:
285             DO j=2-Oly,sNy+Oly
286              DO i=2-Olx,sNx+Olx-1
287                vrtDx(i,j) = (vort3(i+1,j)-vort3(i,j))
288         &                  *recip_DXG(i,j,bi,bj)
289         &                  *maskS(i,j,k,bi,bj)
290              ENDDO
291             ENDDO
292    C-       gradient in y direction:
293             DO j=2-Oly,sNy+Oly-1
294              DO i=2-Olx,sNx+Olx
295                vrtDy(i,j) = (vort3(i,j+1)-vort3(i,j))
296         &                  *recip_DYG(i,j,bi,bj)
297         &                  *maskW(i,j,k,bi,bj)
298              ENDDO
299             ENDDO
300    
301           ENDIF
302    
303         DO j=2-Oly,sNy+Oly-1         DO j=2-Oly,sNy+Oly-1
304          DO i=2-Olx,sNx+Olx-1          DO i=2-Olx,sNx+Olx-1
305  CCCCCCCCCCCCCCC Divergence Point CalculationsCCCCCCCCCCCCCCCCCCCC  CCCCCCCCCCCCCCC Divergence Point CalculationsCCCCCCCCCCCCCCCCCCCC
306    
307  C These are (powers of) length scales  C These are (powers of) length scales
308           IF (useAreaViscLength) THEN           IF (useAreaViscLength) THEN
309            L2=Ra(i,j,bi,bj)            L2=rA(i,j,bi,bj)
310           ELSE           ELSE
311            L2=2. _d 0/((recip_DXF(I,J,bi,bj)**2+recip_DYF(I,J,bi,bj)**2))            L2=2. _d 0/((recip_DXF(I,J,bi,bj)**2+recip_DYF(I,J,bi,bj)**2))
312           ENDIF           ENDIF
# Line 200  C These are (powers of) length scales Line 317  C These are (powers of) length scales
317           L2rdt=0.25 _d 0*recip_dt*L2           L2rdt=0.25 _d 0*recip_dt*L2
318    
319           IF (useAreaViscLength) THEN           IF (useAreaViscLength) THEN
320            L4rdt=0.125 _d 0*recip_dt*Ra(i,j,bi,bj)**2            L4rdt=0.125 _d 0*recip_dt*rA(i,j,bi,bj)**2
321           ELSE           ELSE
322            L4rdt=recip_dt/( 6. _d 0*(recip_DXF(I,J,bi,bj)**4            L4rdt=recip_dt/( 6. _d 0*(recip_DXF(I,J,bi,bj)**4
323       &                            +recip_DYF(I,J,bi,bj)**4)       &                            +recip_DYF(I,J,bi,bj)**4)
# Line 209  C These are (powers of) length scales Line 326  C These are (powers of) length scales
326           ENDIF           ENDIF
327    
328  C Velocity Reynolds Scale  C Velocity Reynolds Scale
329           Uscl=sqrt(KE(i,j)*L2)*viscAhRe_max           IF ( viscAhRe_max.GT.0. .AND. KE(i,j).GT.0. ) THEN
330           U4scl=sqrt(KE(i,j))*L3*viscA4Re_max             Uscl=sqrt(KE(i,j)*L2)*viscAhRe_max
331             ELSE
332               Uscl=0.
333             ENDIF
334             IF ( viscA4Re_max.GT.0. .AND. KE(i,j).GT.0. ) THEN
335               U4scl=sqrt(KE(i,j))*L3*viscA4Re_max
336             ELSE
337               U4scl=0.
338             ENDIF
339    
340           IF (useFullLeith.and.calcleith) THEN           IF (useFullLeith.and.calcleith) THEN
341  C This is the vector magnitude of the vorticity gradient squared  C This is the vector magnitude of the vorticity gradient squared
342            grdVrt=0.25 _d 0*(            grdVrt=0.25 _d 0*( (vrtDx(i,j+1)*vrtDx(i,j+1)
343       &     ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2       &                        + vrtDx(i,j)*vrtDx(i,j) )
344       &     +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2       &                     + (vrtDy(i+1,j)*vrtDy(i+1,j)
345       &     +((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)  
346    
347  C This is the vector magnitude of grad (div.v) squared  C This is the vector magnitude of grad (div.v) squared
348  C Using it in Leith serves to damp instabilities in w.  C Using it in Leith serves to damp instabilities in w.
349            grdDiv=0.25 _d 0*(            grdDiv=0.25 _d 0*( (divDx(i+1,j)*divDx(i+1,j)
350       &     ((hDiv(i+1,j)-hDiv(i,j))*recip_DXC(i+1,j,bi,bj))**2       &                        + divDx(i,j)*divDx(i,j) )
351       &     +((hDiv(i,j+1)-hDiv(i,j))*recip_DYC(i,j+1,bi,bj))**2       &                     + (divDy(i,j+1)*divDy(i,j+1)
352       &     +((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)  
353    
354            viscAh_DLth(i,j)=            viscAh_DLth(i,j)=
355       &     sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3       &     sqrt(leith2fac*grdVrt+leithD2fac*grdDiv)*L3
356            viscA4_DLth(i,j)=0.125 _d 0*            viscA4_DLth(i,j)=
357       &     sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5       &     sqrt(leith4fac*grdVrt+leithD4fac*grdDiv)*L5
358            viscAh_DLthd(i,j)=            viscAh_DLthd(i,j)=
359       &     sqrt(viscC2leithD**2*grdDiv)*L3       &     sqrt(leithD2fac*grdDiv)*L3
360            viscA4_DLthd(i,j)=0.125 _d 0*            viscA4_DLthd(i,j)=
361       &     sqrt(viscC4leithD**2*grdDiv)*L5       &     sqrt(leithD4fac*grdDiv)*L5
362           ELSEIF (calcleith) THEN           ELSEIF (calcleith) THEN
363  C but this approximation will work on cube  C but this approximation will work on cube
364  c (and differs by as much as 4X)  c (and differs by as much as 4X)
365            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)) )
366            grdVrt=max(grdVrt,            grdVrt=max( grdVrt, abs(vrtDy(i+1,j)) )
367       &     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)))  
368    
369  c This approximation is good to the same order as above...  c This approximation is good to the same order as above...
370            viscAh_Dlth(i,j)=            grdDiv=max( abs(divDx(i+1,j)), abs(divDx(i,j)) )
371       &      (viscC2leith*grdVrt+(viscC2leithD*grdDiv))*L3            grdDiv=max( grdDiv, abs(divDy(i,j+1)) )
372            viscA4_Dlth(i,j)=0.125 _d 0*            grdDiv=max( grdDiv, abs(divDy(i,j))   )
373       &      (viscC4leith*grdVrt+(viscC4leithD*grdDiv))*L5  
374            viscAh_DlthD(i,j)=            viscAh_Dlth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3
375       &      ((viscC2leithD*grdDiv))*L3            viscA4_Dlth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5
376            viscA4_DlthD(i,j)=0.125 _d 0*            viscAh_DlthD(i,j)=((leithD2fac*grdDiv))*L3
377       &      ((viscC4leithD*grdDiv))*L5            viscA4_DlthD(i,j)=((leithD4fac*grdDiv))*L5
378           ELSE           ELSE
379            viscAh_Dlth(i,j)=0. _d 0            viscAh_Dlth(i,j)=0. _d 0
380            viscA4_Dlth(i,j)=0. _d 0            viscA4_Dlth(i,j)=0. _d 0
# Line 273  c This approximation is good to the same Line 382  c This approximation is good to the same
382            viscA4_DlthD(i,j)=0. _d 0            viscA4_DlthD(i,j)=0. _d 0
383           ENDIF           ENDIF
384    
385    
386           IF (calcsmag) THEN           IF (calcsmag) THEN
387            viscAh_DSmg(i,j)=L2            viscAh_DSmg(i,j)=L2
388       &       *sqrt(tension(i,j)**2       &       *sqrt(tension(i,j)**2
# Line 304  C  BiHarmonic on Div.u points Line 414  C  BiHarmonic on Div.u points
414  CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC  CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC
415  C These are (powers of) length scales  C These are (powers of) length scales
416           IF (useAreaViscLength) THEN           IF (useAreaViscLength) THEN
417            L2=RaZ(i,j,bi,bj)            L2=rAz(i,j,bi,bj)
418           ELSE           ELSE
419            L2=2._d 0/((recip_DXV(I,J,bi,bj)**2+recip_DYU(I,J,bi,bj)**2))            L2=2. _d 0/((recip_DXV(I,J,bi,bj)**2+recip_DYU(I,J,bi,bj)**2))
420           ENDIF           ENDIF
421    
422           L3=(L2**1.5)           L3=(L2**1.5)
# Line 315  C These are (powers of) length scales Line 425  C These are (powers of) length scales
425    
426           L2rdt=0.25 _d 0*recip_dt*L2           L2rdt=0.25 _d 0*recip_dt*L2
427           IF (useAreaViscLength) THEN           IF (useAreaViscLength) THEN
428            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
429           ELSE           ELSE
430            L4rdt=recip_dt/            L4rdt=recip_dt/
431       &     ( 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)
432       &      +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))
433           ENDIF           ENDIF
434    
435  C Velocity Reynolds Scale  C Velocity Reynolds Scale (Pb here at CS-grid corners !)
436           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
437       &         *L2)*viscAhRe_max             keZpt=0.25 _d 0*( (KE(i,j)+KE(i-1,j-1))
438           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)) )
439       &         *L3*viscA4Re_max             IF ( keZpt.GT.0. ) THEN
440                 Uscl = sqrt(keZpt*L2)*viscAhRe_max
441                 U4scl= sqrt(keZpt)*L3*viscA4Re_max
442               ELSE
443                 Uscl =0.
444                 U4scl=0.
445               ENDIF
446             ELSE
447               Uscl =0.
448               U4scl=0.
449             ENDIF
450    
451  C This is the vector magnitude of the vorticity gradient squared  C This is the vector magnitude of the vorticity gradient squared
452           IF (useFullLeith.and.calcleith) THEN           IF (useFullLeith.and.calcleith) THEN
453            grdVrt=0.25 _d 0*(            grdVrt=0.25 _d 0*( (vrtDx(i-1,j)*vrtDx(i-1,j)
454       &     ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2       &                        + vrtDx(i,j)*vrtDx(i,j) )
455       &     +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2       &                     + (vrtDy(i,j-1)*vrtDy(i,j-1)
456       &     +((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)  
457    
458  C This is the vector magnitude of grad(div.v) squared  C This is the vector magnitude of grad(div.v) squared
459            grdDiv=0.25 _d 0*(            grdDiv=0.25 _d 0*( (divDx(i,j-1)*divDx(i,j-1)
460       &     ((hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj))**2       &                        + divDx(i,j)*divDx(i,j) )
461       &     +((hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj))**2       &                     + (divDy(i-1,j)*divDy(i-1,j)
462       &     +((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)  
463    
464            viscAh_ZLth(i,j)=            viscAh_ZLth(i,j)=
465       &     sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3       &     sqrt(leith2fac*grdVrt+leithD2fac*grdDiv)*L3
466            viscA4_ZLth(i,j)=0.125 _d 0*            viscA4_ZLth(i,j)=
467       &     sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5       &     sqrt(leith4fac*grdVrt+leithD4fac*grdDiv)*L5
468            viscAh_ZLthD(i,j)=            viscAh_ZLthD(i,j)=
469       &     sqrt(viscC2leithD**2*grdDiv)*L3       &     sqrt(leithD2fac*grdDiv)*L3
470            viscA4_ZLthD(i,j)=0.125 _d 0*            viscA4_ZLthD(i,j)=
471       &     sqrt(viscC4leithD**2*grdDiv)*L5       &     sqrt(leithD4fac*grdDiv)*L5
472    
473           ELSEIF (calcleith) THEN           ELSEIF (calcleith) THEN
474  C but this approximation will work on cube (and differs by 4X)  C but this approximation will work on cube (and differs by 4X)
475            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)) )
476            grdVrt=max(grdVrt,            grdVrt=max( grdVrt, abs(vrtDy(i,j-1)) )
477       &     abs((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj)))            grdVrt=max( grdVrt, abs(vrtDy(i,j))   )
478            grdVrt=max(grdVrt,  
479       &     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)) )
480            grdVrt=max(grdVrt,            grdDiv=max( grdDiv, abs(divDy(i,j))   )
481       &     abs((vort3(i,j-1)-vort3(i,j))*recip_DYG(i,j-1,bi,bj)))            grdDiv=max( grdDiv, abs(divDy(i-1,j)) )
482    
483            grdDiv=abs((hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj))            viscAh_ZLth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3
484            grdDiv=max(grdDiv,            viscA4_ZLth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5
485       &     abs((hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj)))            viscAh_ZLthD(i,j)=(leithD2fac*grdDiv)*L3
486            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  
487           ELSE           ELSE
488            viscAh_ZLth(i,j)=0. _d 0            viscAh_ZLth(i,j)=0. _d 0
489            viscA4_ZLth(i,j)=0. _d 0            viscA4_ZLth(i,j)=0. _d 0
# Line 409  C  BiHarmonic on Zeta points Line 517  C  BiHarmonic on Zeta points
517           viscA4_Z(i,j)=min(viscA4_ZMax(i,j),viscA4_Z(i,j))           viscA4_Z(i,j)=min(viscA4_ZMax(i,j),viscA4_Z(i,j))
518          ENDDO          ENDDO
519         ENDDO         ENDDO
520    cph(
521    #else
522           STOP 'useVariableViscosity not implemented for ADJOINT'
523    #endif /* ndef ALLOW_AUTODIFF_TAMC */
524    cph)
525        ELSE        ELSE
526         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
527          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
Legend:
Removed from v.1.11  
changed lines
  Added in v.1.21
  ViewVC Help
Powered by ViewVC 1.1.22