/[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.33 by heimbach, Fri Apr 18 18:11:43 2008 UTC
# 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 34  C     This allows the same value of the Line 34  C     This allows the same value of the
34  C     for roughly similar results with biharmonic and harmonic  C     for roughly similar results with biharmonic and harmonic
35  C  C
36  C     LIMITERS -- limit min and max values of viscosities  C     LIMITERS -- limit min and max values of viscosities
37  C     viscAhRemax is min value for grid point harmonic Reynolds num  C     viscAhReMax is min value for grid point harmonic Reynolds num
38  C      harmonic viscosity>sqrt(2*KE)*L/viscAhRemax  C      harmonic viscosity>sqrt(2*KE)*L/viscAhReMax
39  C  C
40  C     viscA4Remax is min value for grid point biharmonic Reynolds num  C     viscA4ReMax is min value for grid point biharmonic Reynolds num
41  C      biharmonic viscosity>sqrt(2*KE)*L**3/8/viscA4Remax  C      biharmonic viscosity>sqrt(2*KE)*L**3/8/viscA4ReMax
42  C  C
43  C     viscAhgridmax is CFL stability limiter for harmonic viscosity  C     viscAhgridmax is CFL stability limiter for harmonic viscosity
44  C      harmonic viscosity<0.25*viscAhgridmax*L**2/deltaT  C      harmonic viscosity<0.25*viscAhgridmax*L**2/deltaT
# Line 47  C     viscA4gridmax is CFL stability lim Line 47  C     viscA4gridmax is CFL stability lim
47  C      biharmonic viscosity<viscA4gridmax*L**4/32/deltaT (approx)  C      biharmonic viscosity<viscA4gridmax*L**4/32/deltaT (approx)
48  C  C
49  C     viscAhgridmin and viscA4gridmin are lower limits for viscosity:  C     viscAhgridmin and viscA4gridmin are lower limits for viscosity:
50  C      harmonic viscosity>0.25*viscAhgridmax*L**2/deltaT  C       harmonic viscosity>0.25*viscAhgridmin*L**2/deltaT
51  C      biharmonic viscosity>viscA4gridmax*L**4/32/deltaT (approx)  C       biharmonic viscosity>viscA4gridmin*L**4/32/deltaT (approx)
52    
53    
54  C  C
55  C     RECOMMENDED VALUES  C     RECOMMENDED VALUES
56  C     viscC2Leith=?  C     viscC2Leith=1-3
57  C     viscC2LeithD=?  C     viscC2LeithD=1-3
58  C     viscC4Leith=?  C     viscC4Leith=1-3
59  C     viscC4LeithD=?  C     viscC4LeithD=1.5-3
60  C     viscC2smag=2.2-4 (Griffies and Hallberg,2000)  C     viscC2smag=2.2-4 (Griffies and Hallberg,2000)
61  C               0.2-0.9 (Smagorinsky,1993)  C               0.2-0.9 (Smagorinsky,1993)
62  C     viscC4smag=2.2-4 (Griffies and Hallberg,2000)  C     viscC4smag=2.2-4 (Griffies and Hallberg,2000)
63  C     viscAhRemax>=1, (<2 suppresses a computational mode)  C     viscAhReMax>=1, (<2 suppresses a computational mode)
64  C     viscA4Remax>=1, (<2 suppresses a computational mode)  C     viscA4ReMax>=1, (<2 suppresses a computational mode)
65  C     viscAhgridmax=1  C     viscAhgridmax=1
66  C     viscA4gridmax=1  C     viscA4gridmax=1
67  C     viscAhgrid<1  C     viscAhgrid<1
# Line 72  C     == Global variables == Line 74  C     == Global variables ==
74  #include "GRID.h"  #include "GRID.h"
75  #include "EEPARAMS.h"  #include "EEPARAMS.h"
76  #include "PARAMS.h"  #include "PARAMS.h"
77    #ifdef ALLOW_NONHYDROSTATIC
78    #include "NH_VARS.h"
79    #endif
80    #ifdef ALLOW_AUTODIFF_TAMC
81    #include "tamc.h"
82    #include "tamc_keys.h"
83    #endif /* ALLOW_AUTODIFF_TAMC */
84    
85  C     == Routine arguments ==  C     == Routine arguments ==
86        INTEGER bi,bj,k        INTEGER bi,bj,k
# Line 90  C     == Routine arguments == Line 99  C     == Routine arguments ==
99    
100  C     == Local variables ==  C     == Local variables ==
101        INTEGER I,J        INTEGER I,J
102    #ifdef ALLOW_NONHYDROSTATIC
103          INTEGER kp1
104    #endif
105          INTEGER lockey
106        _RL smag2fac, smag4fac        _RL smag2fac, smag4fac
107          _RL leith2fac, leith4fac
108          _RL leithD2fac, leithD4fac
109        _RL viscAhRe_max, viscA4Re_max        _RL viscAhRe_max, viscA4Re_max
110        _RL Alin,grdVrt,grdDiv, keZpt        _RL Alin,grdVrt,grdDiv, keZpt
111        _RL recip_dt,L2,L3,L4,L5,L2rdt,L4rdt        _RL recip_dt,L2,L3,L4,L5,L2rdt,L4rdt
112        _RL Uscl,U4scl        _RL Uscl,U4scl
113          _RL divDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
114          _RL divDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
115          _RL vrtDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
116          _RL vrtDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
117        _RL viscAh_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscAh_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
118        _RL viscAh_DMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscAh_DMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
119        _RL viscA4_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscA4_ZMax(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 115  C     == Local variables == Line 134  C     == Local variables ==
134        _RL viscAh_DSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscAh_DSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
135        _RL viscA4_ZSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscA4_ZSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
136        _RL viscA4_DSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL viscA4_DSmg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
137        LOGICAL calcLeith,calcSmag        LOGICAL calcLeith, calcSmag
138    
139    #ifdef ALLOW_AUTODIFF_TAMC
140              act1 = bi - myBxLo(myThid)
141              max1 = myBxHi(myThid) - myBxLo(myThid) + 1
142              act2 = bj - myByLo(myThid)
143              max2 = myByHi(myThid) - myByLo(myThid) + 1
144              act3 = myThid - 1
145              max3 = nTx*nTy
146              act4 = ikey_dynamics - 1
147              ikey = (act1 + 1) + act2*max1
148         &                      + act3*max1*max2
149         &                      + act4*max1*max2*max3
150              lockey = (ikey-1)*Nr + k
151    #endif /* ALLOW_AUTODIFF_TAMC */
152    
153    C--   Set flags which are used in this S/R and elsewhere :
154        useVariableViscosity=        useVariableViscosity=
155       &      (viscAhGrid.NE.0.)       &      (viscAhGrid.NE.0.)
156       &  .OR.(viscA4Grid.NE.0.)       &  .OR.(viscA4Grid.NE.0.)
# Line 136  C     == Local variables == Line 170  C     == Local variables ==
170       &  .OR.(viscC2leithD.NE.0.)       &  .OR.(viscC2leithD.NE.0.)
171       &  .OR.(viscC2smag.NE.0.)       &  .OR.(viscC2smag.NE.0.)
172    
       IF ((harmonic).and.(viscAhremax.ne.0.)) THEN  
         viscAhre_max=sqrt(2. _d 0)/viscAhRemax  
       ELSE  
         viscAhre_max=0. _d 0  
       ENDIF  
   
173        biharmonic=        biharmonic=
174       &      (viscA4.NE.0.)       &      (viscA4.NE.0.)
175       &  .OR.(viscA4D.NE.0.)       &  .OR.(viscA4D.NE.0.)
# Line 151  C     == Local variables == Line 179  C     == Local variables ==
179       &  .OR.(viscC4leithD.NE.0.)       &  .OR.(viscC4leithD.NE.0.)
180       &  .OR.(viscC4smag.NE.0.)       &  .OR.(viscC4smag.NE.0.)
181    
182        IF ((biharmonic).and.(viscA4remax.ne.0.)) THEN        IF (useVariableViscosity) THEN
183          viscA4re_max=0.125 _d 0*sqrt(2. _d 0)/viscA4Remax  C---- variable viscosity :
       ELSE  
         viscA4re_max=0. _d 0  
       ENDIF  
184    
185        calcleith=         IF ((harmonic).AND.(viscAhReMax.NE.0.)) THEN
186            viscAhRe_max=SQRT(2. _d 0)/viscAhReMax
187           ELSE
188            viscAhRe_max=0. _d 0
189           ENDIF
190    
191           IF ((biharmonic).AND.(viscA4ReMax.NE.0.)) THEN
192            viscA4Re_max=0.125 _d 0*SQRT(2. _d 0)/viscA4ReMax
193           ELSE
194            viscA4Re_max=0. _d 0
195           ENDIF
196    
197           calcLeith=
198       &      (viscC2leith.NE.0.)       &      (viscC2leith.NE.0.)
199       &  .OR.(viscC2leithD.NE.0.)       &  .OR.(viscC2leithD.NE.0.)
200       &  .OR.(viscC4leith.NE.0.)       &  .OR.(viscC4leith.NE.0.)
201       &  .OR.(viscC4leithD.NE.0.)       &  .OR.(viscC4leithD.NE.0.)
202    
203        calcsmag=         calcSmag=
204       &      (viscC2smag.NE.0.)       &      (viscC2smag.NE.0.)
205       &  .OR.(viscC4smag.NE.0.)       &  .OR.(viscC4smag.NE.0.)
206    
207        IF (deltaTmom.NE.0.) THEN         IF (deltaTmom.NE.0.) THEN
208         recip_dt=1. _d 0/deltaTmom          recip_dt=1. _d 0/deltaTmom
209        ELSE         ELSE
210         recip_dt=0. _d 0          recip_dt=0. _d 0
211        ENDIF         ENDIF
212    
213        IF (calcsmag) THEN         IF (calcSmag) THEN
214          smag2fac=(viscC2smag/pi)**2          smag2fac=(viscC2smag/pi)**2
215          smag4fac=0.125 _d 0*(viscC4smag/pi)**2          smag4fac=0.125 _d 0*(viscC4smag/pi)**2
216        ELSE         ELSE
217          smag2fac=0. _d 0          smag2fac=0. _d 0
218          smag4fac=0. _d 0          smag4fac=0. _d 0
219        ENDIF         ENDIF
220    
221           IF (calcLeith) THEN
222            IF (useFullLeith) THEN
223             leith2fac =(viscC2leith /pi)**6
224             leithD2fac=(viscC2leithD/pi)**6
225             leith4fac =0.015625 _d 0*(viscC4leith /pi)**6
226             leithD4fac=0.015625 _d 0*(viscC4leithD/pi)**6
227            ELSE
228             leith2fac =(viscC2leith /pi)**3
229             leithD2fac=(viscC2leithD/pi)**3
230             leith4fac =0.125 _d 0*(viscC4leith /pi)**3
231             leithD4fac=0.125 _d 0*(viscC4leithD/pi)**3
232            ENDIF
233           ELSE
234            leith2fac=0. _d 0
235            leith4fac=0. _d 0
236            leithD2fac=0. _d 0
237            leithD4fac=0. _d 0
238           ENDIF
239    
240    #ifdef ALLOW_AUTODIFF_TAMC
241    cphtest       IF ( calcLeith .OR. calcSmag ) THEN
242    cphtest        STOP 'calcLeith or calcSmag not implemented for ADJOINT'
243    cphtest       ENDIF
244    #endif
245           DO j=1-Oly,sNy+Oly
246            DO i=1-Olx,sNx+Olx
247              viscAh_D(i,j)=viscAhD
248              viscAh_Z(i,j)=viscAhZ
249              viscA4_D(i,j)=viscA4D
250              viscA4_Z(i,j)=viscA4Z
251    c
252              visca4_zsmg(i,j) = 0. _d 0
253              viscah_zsmg(i,j) = 0. _d 0
254    c
255              viscAh_Dlth(i,j) = 0. _d 0
256              viscA4_Dlth(i,j) = 0. _d 0
257              viscAh_DlthD(i,j)= 0. _d 0
258              viscA4_DlthD(i,j)= 0. _d 0
259    c
260              viscAh_DSmg(i,j) = 0. _d 0
261              viscA4_DSmg(i,j) = 0. _d 0
262    c
263              viscAh_ZLth(i,j) = 0. _d 0
264              viscA4_ZLth(i,j) = 0. _d 0
265              viscAh_ZLthD(i,j)= 0. _d 0
266              viscA4_ZLthD(i,j)= 0. _d 0
267            ENDDO
268           ENDDO
269    
270    C-     Initialise to zero gradient of vorticity & divergence:
271           DO j=1-Oly,sNy+Oly
272            DO i=1-Olx,sNx+Olx
273              divDx(i,j) = 0.
274              divDy(i,j) = 0.
275              vrtDx(i,j) = 0.
276              vrtDy(i,j) = 0.
277            ENDDO
278           ENDDO
279    
280           IF (calcLeith) THEN
281    C      horizontal gradient of horizontal divergence:
282    
283    C-       gradient in x direction:
284    cph-exch2#ifndef ALLOW_AUTODIFF_TAMC
285             IF (useCubedSphereExchange) THEN
286    C        to compute d/dx(hDiv), fill corners with appropriate values:
287               CALL FILL_CS_CORNER_TR_RL( .TRUE., .FALSE.,
288         &                                hDiv, bi,bj, myThid )
289             ENDIF
290    cph-exch2#endif
291             DO j=2-Oly,sNy+Oly-1
292              DO i=2-Olx,sNx+Olx-1
293                divDx(i,j) = (hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj)
294              ENDDO
295             ENDDO
296    
297    C-       gradient in y direction:
298    cph-exch2#ifndef ALLOW_AUTODIFF_TAMC
299             IF (useCubedSphereExchange) THEN
300    C        to compute d/dy(hDiv), fill corners with appropriate values:
301               CALL FILL_CS_CORNER_TR_RL(.FALSE., .FALSE.,
302         &                                hDiv, bi,bj, myThid )
303             ENDIF
304    cph-exch2#endif
305             DO j=2-Oly,sNy+Oly-1
306              DO i=2-Olx,sNx+Olx-1
307                divDy(i,j) = (hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj)
308              ENDDO
309             ENDDO
310    
311    C      horizontal gradient of vertical vorticity:
312    C-       gradient in x direction:
313             DO j=2-Oly,sNy+Oly
314              DO i=2-Olx,sNx+Olx-1
315                vrtDx(i,j) = (vort3(i+1,j)-vort3(i,j))
316         &                  *recip_DXG(i,j,bi,bj)
317         &                  *maskS(i,j,k,bi,bj)
318              ENDDO
319             ENDDO
320    C-       gradient in y direction:
321             DO j=2-Oly,sNy+Oly-1
322              DO i=2-Olx,sNx+Olx
323                vrtDy(i,j) = (vort3(i,j+1)-vort3(i,j))
324         &                  *recip_DYG(i,j,bi,bj)
325         &                  *maskW(i,j,k,bi,bj)
326              ENDDO
327             ENDDO
328    
329           ENDIF
330    
331    #ifdef    ALLOW_AUTODIFF_TAMC
332    cphCADJ STORE viscA4_ZSmg(:,:)
333    cphCADJ &     = comlev1_bibj_k , key=lockey, byte=isbyte
334    cphCADJ STORE viscAh_ZSmg(:,:)
335    cphCADJ &     = comlev1_bibj_k , key=lockey, byte=isbyte
336    #endif    /* ALLOW_AUTODIFF_TAMC */
337    
 C     - Viscosity  
       IF (useVariableViscosity) THEN  
338         DO j=2-Oly,sNy+Oly-1         DO j=2-Oly,sNy+Oly-1
339          DO i=2-Olx,sNx+Olx-1          DO i=2-Olx,sNx+Olx-1
340  CCCCCCCCCCCCCCC Divergence Point CalculationsCCCCCCCCCCCCCCCCCCCC  CCCCCCCCCCCCCCC Divergence Point CalculationsCCCCCCCCCCCCCCCCCCCC
341    
342  C These are (powers of) length scales  C These are (powers of) length scales
343           IF (useAreaViscLength) THEN           IF (useAreaViscLength) THEN
344            L2=rA(i,j,bi,bj)            L2=rA(i,j,bi,bj)
345              L4rdt=0.03125 _d 0*recip_dt*L2**2
346           ELSE           ELSE
347            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))
348              L4rdt=recip_dt/( 6. _d 0*(recip_DXF(I,J,bi,bj)**4
349         &                             +recip_DYF(I,J,bi,bj)**4)
350         &                   +8. _d 0*((recip_DXF(I,J,bi,bj)
351         &                             *recip_DYF(I,J,bi,bj))**2) )
352           ENDIF           ENDIF
353           L3=(L2**1.5)           L3=(L2**1.5)
354           L4=(L2**2)           L4=(L2**2)
355           L5=(L2**2.5)           L5=(L2*L3)
356    
357           L2rdt=0.25 _d 0*recip_dt*L2           L2rdt=0.25 _d 0*recip_dt*L2
358    
          IF (useAreaViscLength) THEN  
           L4rdt=0.125 _d 0*recip_dt*rA(i,j,bi,bj)**2  
          ELSE  
           L4rdt=recip_dt/( 6. _d 0*(recip_DXF(I,J,bi,bj)**4  
      &                            +recip_DYF(I,J,bi,bj)**4)  
      &                   +8. _d 0*((recip_DXF(I,J,bi,bj)  
      &                             *recip_DYF(I,J,bi,bj))**2) )  
          ENDIF  
   
359  C Velocity Reynolds Scale  C Velocity Reynolds Scale
360           IF ( viscAhRe_max.GT.0. .AND. KE(i,j).GT.0. ) THEN           IF ( viscAhRe_max.GT.0. .AND. KE(i,j).GT.0. ) THEN
361             Uscl=sqrt(KE(i,j)*L2)*viscAhRe_max             Uscl=SQRT(KE(i,j)*L2)*viscAhRe_max
362           ELSE           ELSE
363             Uscl=0.             Uscl=0.
364           ENDIF           ENDIF
365           IF ( viscA4Re_max.GT.0. .AND. KE(i,j).GT.0. ) THEN           IF ( viscA4Re_max.GT.0. .AND. KE(i,j).GT.0. ) THEN
366             U4scl=sqrt(KE(i,j))*L3*viscA4Re_max             U4scl=SQRT(KE(i,j))*L3*viscA4Re_max
367           ELSE           ELSE
368             U4scl=0.             U4scl=0.
369           ENDIF           ENDIF
370    
371           IF (useFullLeith.and.calcleith) THEN  cph-leith#ifndef ALLOW_AUTODIFF_TAMC
372    #ifndef AUTODIFF_DISABLE_LEITH
373             IF (useFullLeith.AND.calcLeith) THEN
374  C This is the vector magnitude of the vorticity gradient squared  C This is the vector magnitude of the vorticity gradient squared
375            grdVrt=0.25 _d 0*(            grdVrt=0.25 _d 0*( (vrtDx(i,j+1)*vrtDx(i,j+1)
376       &     ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2       &                        + vrtDx(i,j)*vrtDx(i,j) )
377       &     +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2       &                     + (vrtDy(i+1,j)*vrtDy(i+1,j)
378       &     +((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)  
379    
380  C This is the vector magnitude of grad (div.v) squared  C This is the vector magnitude of grad (div.v) squared
381  C Using it in Leith serves to damp instabilities in w.  C Using it in Leith serves to damp instabilities in w.
382            grdDiv=0.25 _d 0*(            grdDiv=0.25 _d 0*( (divDx(i+1,j)*divDx(i+1,j)
383       &     ((hDiv(i+1,j)-hDiv(i,j))*recip_DXC(i+1,j,bi,bj))**2       &                        + divDx(i,j)*divDx(i,j) )
384       &     +((hDiv(i,j+1)-hDiv(i,j))*recip_DYC(i,j+1,bi,bj))**2       &                     + (divDy(i,j+1)*divDy(i,j+1)
385       &     +((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)  
386    
387            viscAh_DLth(i,j)=            viscAh_DLth(i,j)=
388       &     sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3       &     SQRT(leith2fac*grdVrt+leithD2fac*grdDiv)*L3
389            viscA4_DLth(i,j)=0.125 _d 0*            viscA4_DLth(i,j)=
390       &     sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5       &     SQRT(leith4fac*grdVrt+leithD4fac*grdDiv)*L5
391            viscAh_DLthd(i,j)=            viscAh_DLthd(i,j)=
392       &     sqrt(viscC2leithD**2*grdDiv)*L3       &     SQRT(leithD2fac*grdDiv)*L3
393            viscA4_DLthd(i,j)=0.125 _d 0*            viscA4_DLthd(i,j)=
394       &     sqrt(viscC4leithD**2*grdDiv)*L5       &     SQRT(leithD4fac*grdDiv)*L5
395           ELSEIF (calcleith) THEN           ELSEIF (calcLeith) THEN
396  C but this approximation will work on cube  C but this approximation will work on cube
397  c (and differs by as much as 4X)  c (and differs by as much as 4X)
398            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)) )
399            grdVrt=max(grdVrt,            grdVrt=MAX( grdVrt, ABS(vrtDy(i+1,j)) )
400       &     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)))  
401    
402  c This approximation is good to the same order as above...  c This approximation is good to the same order as above...
403            viscAh_Dlth(i,j)=            grdDiv=MAX( ABS(divDx(i+1,j)), ABS(divDx(i,j)) )
404       &      (viscC2leith*grdVrt+(viscC2leithD*grdDiv))*L3            grdDiv=MAX( grdDiv, ABS(divDy(i,j+1)) )
405            viscA4_Dlth(i,j)=0.125 _d 0*            grdDiv=MAX( grdDiv, ABS(divDy(i,j))   )
406       &      (viscC4leith*grdVrt+(viscC4leithD*grdDiv))*L5  
407            viscAh_DlthD(i,j)=            viscAh_Dlth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3
408       &      ((viscC2leithD*grdDiv))*L3            viscA4_Dlth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5
409            viscA4_DlthD(i,j)=0.125 _d 0*            viscAh_DlthD(i,j)=((leithD2fac*grdDiv))*L3
410       &      ((viscC4leithD*grdDiv))*L5            viscA4_DlthD(i,j)=((leithD4fac*grdDiv))*L5
411           ELSE           ELSE
412            viscAh_Dlth(i,j)=0. _d 0            viscAh_Dlth(i,j)=0. _d 0
413            viscA4_Dlth(i,j)=0. _d 0            viscA4_Dlth(i,j)=0. _d 0
# Line 281  c This approximation is good to the same Line 415  c This approximation is good to the same
415            viscA4_DlthD(i,j)=0. _d 0            viscA4_DlthD(i,j)=0. _d 0
416           ENDIF           ENDIF
417    
418           IF (calcsmag) THEN           IF (calcSmag) THEN
419            viscAh_DSmg(i,j)=L2            viscAh_DSmg(i,j)=L2
420       &       *sqrt(tension(i,j)**2       &       *SQRT(tension(i,j)**2
421       &       +0.25 _d 0*(strain(i+1, j )**2+strain( i ,j+1)**2       &       +0.25 _d 0*(strain(i+1, j )**2+strain( i ,j+1)**2
422       &                  +strain(i  , j )**2+strain(i+1,j+1)**2))       &                  +strain(i  , j )**2+strain(i+1,j+1)**2))
423            viscA4_DSmg(i,j)=smag4fac*L2*viscAh_DSmg(i,j)            viscA4_DSmg(i,j)=smag4fac*L2*viscAh_DSmg(i,j)
# Line 292  c This approximation is good to the same Line 426  c This approximation is good to the same
426            viscAh_DSmg(i,j)=0. _d 0            viscAh_DSmg(i,j)=0. _d 0
427            viscA4_DSmg(i,j)=0. _d 0            viscA4_DSmg(i,j)=0. _d 0
428           ENDIF           ENDIF
429    #endif /* AUTODIFF_DISABLE_LEITH */
430    
431  C  Harmonic on Div.u points  C  Harmonic on Div.u points
432           Alin=viscAhD+viscAhGrid*L2rdt           Alin=viscAhD+viscAhGrid*L2rdt
433       &          +viscAh_DLth(i,j)+viscAh_DSmg(i,j)       &          +viscAh_DLth(i,j)+viscAh_DSmg(i,j)
434           viscAh_DMin(i,j)=max(viscAhGridMin*L2rdt,Uscl)           viscAh_DMin(i,j)=MAX(viscAhGridMin*L2rdt,Uscl)
435           viscAh_D(i,j)=max(viscAh_DMin(i,j),Alin)           viscAh_D(i,j)=MAX(viscAh_DMin(i,j),Alin)
436           viscAh_DMax(i,j)=min(viscAhGridMax*L2rdt,viscAhMax)           viscAh_DMax(i,j)=MIN(viscAhGridMax*L2rdt,viscAhMax)
437           viscAh_D(i,j)=min(viscAh_DMax(i,j),viscAh_D(i,j))           viscAh_D(i,j)=MIN(viscAh_DMax(i,j),viscAh_D(i,j))
438    
439  C  BiHarmonic on Div.u points  C  BiHarmonic on Div.u points
440           Alin=viscA4D+viscA4Grid*L4rdt           Alin=viscA4D+viscA4Grid*L4rdt
441       &          +viscA4_DLth(i,j)+viscA4_DSmg(i,j)       &          +viscA4_DLth(i,j)+viscA4_DSmg(i,j)
442           viscA4_DMin(i,j)=max(viscA4GridMin*L4rdt,U4scl)           viscA4_DMin(i,j)=MAX(viscA4GridMin*L4rdt,U4scl)
443           viscA4_D(i,j)=max(viscA4_DMin(i,j),Alin)           viscA4_D(i,j)=MAX(viscA4_DMin(i,j),Alin)
444           viscA4_DMax(i,j)=min(viscA4GridMax*L4rdt,viscA4Max)           viscA4_DMax(i,j)=MIN(viscA4GridMax*L4rdt,viscA4Max)
445           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))
446    
447    #ifdef ALLOW_NONHYDROSTATIC
448    C--   Pass Viscosities to calc_gw, if constant, not necessary
449    
450             kp1  = MIN(k+1,Nr)
451    
452             IF ( k.EQ.1 ) THEN
453    C     Prepare for next level (next call)
454              viscAh_W(i,j,kp1,bi,bj)=0.5*viscAh_D(i,j)
455              viscA4_W(i,j,kp1,bi,bj)=0.5*viscA4_D(i,j)
456    
457    C     These values dont get used
458              viscAh_W(i,j,k,bi,bj)=viscAh_D(i,j)
459              viscA4_W(i,j,k,bi,bj)=viscA4_D(i,j)
460    
461             ELSEIF ( k.EQ.Nr ) THEN
462              viscAh_W(i,j,k,bi,bj)=viscAh_W(i,j,k,bi,bj)+0.5*viscAh_D(i,j)
463              viscA4_W(i,j,k,bi,bj)=viscA4_W(i,j,k,bi,bj)+0.5*viscA4_D(i,j)
464    
465             ELSE
466    C     Prepare for next level (next call)
467              viscAh_W(i,j,kp1,bi,bj)=0.5*viscAh_D(i,j)
468              viscA4_W(i,j,kp1,bi,bj)=0.5*viscA4_D(i,j)
469    
470    C     Note that previous call of this function has already added half.
471              viscAh_W(i,j,k,bi,bj)=viscAh_W(i,j,k,bi,bj)+0.5*viscAh_D(i,j)
472              viscA4_W(i,j,k,bi,bj)=viscA4_W(i,j,k,bi,bj)+0.5*viscA4_D(i,j)
473    
474             ENDIF
475    #endif /* ALLOW_NONHYDROSTATIC */
476    
477  CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC  CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC
478  C These are (powers of) length scales  C These are (powers of) length scales
479           IF (useAreaViscLength) THEN           IF (useAreaViscLength) THEN
480            L2=rAz(i,j,bi,bj)            L2=rAz(i,j,bi,bj)
481              L4rdt=0.125 _d 0*recip_dt*rAz(i,j,bi,bj)**2
482           ELSE           ELSE
483            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))
484              L4rdt=recip_dt/
485         &     ( 6. _d 0*(recip_DXV(I,J,bi,bj)**4+recip_DYU(I,J,bi,bj)**4)
486         &      +8. _d 0*((recip_DXV(I,J,bi,bj)*recip_DYU(I,J,bi,bj))**2))
487           ENDIF           ENDIF
488    
489           L3=(L2**1.5)           L3=(L2**1.5)
490           L4=(L2**2)           L4=(L2**2)
491           L5=(L2**2.5)           L5=(L2*L3)
492    
493           L2rdt=0.25 _d 0*recip_dt*L2           L2rdt=0.25 _d 0*recip_dt*L2
          IF (useAreaViscLength) THEN  
           L4rdt=0.125 _d 0*recip_dt*rAz(i,j,bi,bj)**2  
          ELSE  
           L4rdt=recip_dt/  
      &     ( 6. _d 0*(recip_DXV(I,J,bi,bj)**4+recip_DYU(I,J,bi,bj)**4)  
      &      +8. _d 0*((recip_DXV(I,J,bi,bj)*recip_DYU(I,J,bi,bj))**2))  
          ENDIF  
494    
495  C Velocity Reynolds Scale (Pb here at CS-grid corners !)  C Velocity Reynolds Scale (Pb here at CS-grid corners !)
496           IF ( viscAhRe_max.GT.0. .OR. viscA4Re_max.GT.0. ) THEN           IF ( viscAhRe_max.GT.0. .OR. viscA4Re_max.GT.0. ) THEN
497             keZpt=0.25 _d 0*( (KE(i,j)+KE(i-1,j-1))             keZpt=0.25 _d 0*( (KE(i,j)+KE(i-1,j-1))
498       &                      +(KE(i-1,j)+KE(i,j-1)) )       &                      +(KE(i-1,j)+KE(i,j-1)) )
499             IF ( keZpt.GT.0. ) THEN             IF ( keZpt.GT.0. ) THEN
500               Uscl = sqrt(keZpt*L2)*viscAhRe_max               Uscl = SQRT(keZpt*L2)*viscAhRe_max
501               U4scl= sqrt(keZpt)*L3*viscA4Re_max               U4scl= SQRT(keZpt)*L3*viscA4Re_max
502             ELSE             ELSE
503               Uscl =0.               Uscl =0.
504               U4scl=0.               U4scl=0.
# Line 346  C Velocity Reynolds Scale (Pb here at CS Line 508  C Velocity Reynolds Scale (Pb here at CS
508             U4scl=0.             U4scl=0.
509           ENDIF           ENDIF
510    
511    #ifndef AUTODIFF_DISABLE_LEITH
512  C This is the vector magnitude of the vorticity gradient squared  C This is the vector magnitude of the vorticity gradient squared
513           IF (useFullLeith.and.calcleith) THEN           IF (useFullLeith.AND.calcLeith) THEN
514            grdVrt=0.25 _d 0*(            grdVrt=0.25 _d 0*( (vrtDx(i-1,j)*vrtDx(i-1,j)
515       &     ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2       &                        + vrtDx(i,j)*vrtDx(i,j) )
516       &     +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2       &                     + (vrtDy(i,j-1)*vrtDy(i,j-1)
517       &     +((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)  
518    
519  C This is the vector magnitude of grad(div.v) squared  C This is the vector magnitude of grad(div.v) squared
520            grdDiv=0.25 _d 0*(            grdDiv=0.25 _d 0*( (divDx(i,j-1)*divDx(i,j-1)
521       &     ((hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj))**2       &                        + divDx(i,j)*divDx(i,j) )
522       &     +((hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj))**2       &                     + (divDy(i-1,j)*divDy(i-1,j)
523       &     +((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)  
524    
525            viscAh_ZLth(i,j)=            viscAh_ZLth(i,j)=
526       &     sqrt(viscC2leith**2*grdVrt+viscC2leithD**2*grdDiv)*L3       &     SQRT(leith2fac*grdVrt+leithD2fac*grdDiv)*L3
527            viscA4_ZLth(i,j)=0.125 _d 0*            viscA4_ZLth(i,j)=
528       &     sqrt(viscC4leith**2*grdVrt+viscC4leithD**2*grdDiv)*L5       &     SQRT(leith4fac*grdVrt+leithD4fac*grdDiv)*L5
529            viscAh_ZLthD(i,j)=            viscAh_ZLthD(i,j)=
530       &     sqrt(viscC2leithD**2*grdDiv)*L3       &     SQRT(leithD2fac*grdDiv)*L3
531            viscA4_ZLthD(i,j)=0.125 _d 0*            viscA4_ZLthD(i,j)=
532       &     sqrt(viscC4leithD**2*grdDiv)*L5       &     SQRT(leithD4fac*grdDiv)*L5
533    
534           ELSEIF (calcleith) THEN           ELSEIF (calcLeith) THEN
535  C but this approximation will work on cube (and differs by 4X)  C but this approximation will work on cube (and differs by 4X)
536            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)) )
537            grdVrt=max(grdVrt,            grdVrt=MAX( grdVrt, ABS(vrtDy(i,j-1)) )
538       &     abs((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj)))            grdVrt=MAX( grdVrt, ABS(vrtDy(i,j))   )
539            grdVrt=max(grdVrt,  
540       &     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)) )
541            grdVrt=max(grdVrt,            grdDiv=MAX( grdDiv, ABS(divDy(i,j))   )
542       &     abs((vort3(i,j-1)-vort3(i,j))*recip_DYG(i,j-1,bi,bj)))            grdDiv=MAX( grdDiv, ABS(divDy(i-1,j)) )
543    
544            grdDiv=abs((hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj))            viscAh_ZLth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3
545            grdDiv=max(grdDiv,            viscA4_ZLth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5
546       &     abs((hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj)))            viscAh_ZLthD(i,j)=(leithD2fac*grdDiv)*L3
547            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  
548           ELSE           ELSE
549            viscAh_ZLth(i,j)=0. _d 0            viscAh_ZLth(i,j)=0. _d 0
550            viscA4_ZLth(i,j)=0. _d 0            viscA4_ZLth(i,j)=0. _d 0
# Line 401  C but this approximation will work on cu Line 552  C but this approximation will work on cu
552            viscA4_ZLthD(i,j)=0. _d 0            viscA4_ZLthD(i,j)=0. _d 0
553           ENDIF           ENDIF
554    
555           IF (calcsmag) THEN           IF (calcSmag) THEN
556            viscAh_ZSmg(i,j)=L2            viscAh_ZSmg(i,j)=L2
557       &      *sqrt(strain(i,j)**2       &      *SQRT(strain(i,j)**2
558       &        +0.25 _d 0*(tension( i , j )**2+tension( i ,j-1)**2       &        +0.25 _d 0*(tension( i , j )**2+tension( i ,j-1)**2
559       &                   +tension(i-1, j )**2+tension(i-1,j-1)**2))       &                   +tension(i-1, j )**2+tension(i-1,j-1)**2))
560            viscA4_ZSmg(i,j)=smag4fac*L2*viscAh_ZSmg(i,j)            viscA4_ZSmg(i,j)=smag4fac*L2*viscAh_ZSmg(i,j)
561            viscAh_ZSmg(i,j)=smag2fac*viscAh_ZSmg(i,j)            viscAh_ZSmg(i,j)=smag2fac*viscAh_ZSmg(i,j)
562           ENDIF           ENDIF
563    #endif /* AUTODIFF_DISABLE_LEITH */
564    
565  C  Harmonic on Zeta points  C  Harmonic on Zeta points
566           Alin=viscAhZ+viscAhGrid*L2rdt           Alin=viscAhZ+viscAhGrid*L2rdt
567       &           +viscAh_ZLth(i,j)+viscAh_ZSmg(i,j)       &           +viscAh_ZLth(i,j)+viscAh_ZSmg(i,j)
568           viscAh_ZMin(i,j)=max(viscAhGridMin*L2rdt,Uscl)           viscAh_ZMin(i,j)=MAX(viscAhGridMin*L2rdt,Uscl)
569           viscAh_Z(i,j)=max(viscAh_ZMin(i,j),Alin)           viscAh_Z(i,j)=MAX(viscAh_ZMin(i,j),Alin)
570           viscAh_ZMax(i,j)=min(viscAhGridMax*L2rdt,viscAhMax)           viscAh_ZMax(i,j)=MIN(viscAhGridMax*L2rdt,viscAhMax)
571           viscAh_Z(i,j)=min(viscAh_ZMax(i,j),viscAh_Z(i,j))           viscAh_Z(i,j)=MIN(viscAh_ZMax(i,j),viscAh_Z(i,j))
572    
573  C  BiHarmonic on Zeta points  C  BiHarmonic on Zeta points
574           Alin=viscA4Z+viscA4Grid*L4rdt           Alin=viscA4Z+viscA4Grid*L4rdt
575       &           +viscA4_ZLth(i,j)+viscA4_ZSmg(i,j)       &           +viscA4_ZLth(i,j)+viscA4_ZSmg(i,j)
576           viscA4_ZMin(i,j)=max(viscA4GridMin*L4rdt,U4scl)           viscA4_ZMin(i,j)=MAX(viscA4GridMin*L4rdt,U4scl)
577           viscA4_Z(i,j)=max(viscA4_ZMin(i,j),Alin)           viscA4_Z(i,j)=MAX(viscA4_ZMin(i,j),Alin)
578           viscA4_ZMax(i,j)=min(viscA4GridMax*L4rdt,viscA4Max)           viscA4_ZMax(i,j)=MIN(viscA4GridMax*L4rdt,viscA4Max)
579           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))
580          ENDDO          ENDDO
581         ENDDO         ENDDO
582    
583        ELSE        ELSE
584    C---- use constant viscosity (useVariableViscosity=F):
585    
586         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
587          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
588           viscAh_D(i,j)=viscAhD           viscAh_D(i,j)=viscAhD
# Line 436  C  BiHarmonic on Zeta points Line 591  C  BiHarmonic on Zeta points
591           viscA4_Z(i,j)=viscA4Z           viscA4_Z(i,j)=viscA4Z
592          ENDDO          ENDDO
593         ENDDO         ENDDO
594    
595    C---- variable/constant viscosity : end if/else block
596        ENDIF        ENDIF
597    
598  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
Legend:
Removed from v.1.15  
changed lines
  Added in v.1.33
  ViewVC Help
Powered by ViewVC 1.1.22