/[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.19 by baylor, Mon Oct 10 19:49:48 2005 UTC revision 1.37 by heimbach, Mon Apr 6 23:47:06 2009 UTC
# 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=1-3  C     viscC2Leith=1-3
57  C     viscC2LeithD=1-3  C     viscC2LeithD=1-3
58  C     viscC4Leith=1-3  C     viscC4Leith=1-3
59  C     viscC4LeithD=1.5-3  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_1, lockey_2
106        _RL smag2fac, smag4fac        _RL smag2fac, smag4fac
107        _RL leith2fac, leith4fac        _RL leith2fac, leith4fac
108        _RL leithD2fac, leithD4fac        _RL leithD2fac, leithD4fac
# Line 99  C     == Local variables == Line 112  C     == Local variables ==
112        _RL Uscl,U4scl        _RL Uscl,U4scl
113        _RL divDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL divDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
114        _RL divDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _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 119  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_1 = (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 140  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 155  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 :
184        ELSE  
185          viscA4re_max=0. _d 0         IF ((harmonic).AND.(viscAhReMax.NE.0.)) THEN
186        ENDIF          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=         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         IF (calcLeith) THEN
222          IF (useFullLeith) THEN          IF (useFullLeith) THEN
223           leith2fac =(viscC2leith /pi)**6           leith2fac =(viscC2leith /pi)**6
224           leithD2fac=(viscC2leithD/pi)**6           leithD2fac=(viscC2leithD/pi)**6
# Line 197  C     == Local variables == Line 230  C     == Local variables ==
230           leith4fac =0.125 _d 0*(viscC4leith /pi)**3           leith4fac =0.125 _d 0*(viscC4leith /pi)**3
231           leithD4fac=0.125 _d 0*(viscC4leithD/pi)**3           leithD4fac=0.125 _d 0*(viscC4leithD/pi)**3
232          ENDIF          ENDIF
233        ELSE         ELSE
234          leith2fac=0. _d 0          leith2fac=0. _d 0
235          leith4fac=0. _d 0          leith4fac=0. _d 0
236          leithD2fac=0. _d 0          leithD2fac=0. _d 0
237          leithD4fac=0. _d 0          leithD4fac=0. _d 0
238        ENDIF         ENDIF
239    
240  C     - Viscosity  #ifdef ALLOW_AUTODIFF_TAMC
241        IF (useVariableViscosity) THEN  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      horizontal gradient of horizontal divergence:  C-     Initialise to zero gradient of vorticity & divergence:
271         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
272          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
273            divDx(i,j) = 0.            divDx(i,j) = 0.
274            divDy(i,j) = 0.            divDy(i,j) = 0.
275              vrtDx(i,j) = 0.
276              vrtDy(i,j) = 0.
277          ENDDO          ENDDO
278         ENDDO         ENDDO
279         IF (calcleith) THEN  
280           IF (calcLeith) THEN
281    C      horizontal gradient of horizontal divergence:
282    
283  C-       gradient in x direction:  C-       gradient in x direction:
284  #ifndef ALLOW_AUTODIFF_TAMC  cph-exch2#ifndef ALLOW_AUTODIFF_TAMC
285           IF (useCubedSphereExchange) THEN           IF (useCubedSphereExchange) THEN
286  C        to compute d/dx(hDiv), fill corners with appropriate values:  C        to compute d/dx(hDiv), fill corners with appropriate values:
287             CALL FILL_CS_CORNER_TR_RL( .TRUE., hDiv, bi,bj, myThid )             CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
288         &                                hDiv, bi,bj, myThid )
289           ENDIF           ENDIF
290  #endif  cph-exch2#endif
291           DO j=2-Oly,sNy+Oly-1           DO j=2-Oly,sNy+Oly-1
292            DO i=2-Olx,sNx+Olx-1            DO i=2-Olx,sNx+Olx-1
293              divDx(i,j) = (hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj)              divDx(i,j) = (hDiv(i,j)-hDiv(i-1,j))*recip_DXC(i,j,bi,bj)
# Line 229  C        to compute d/dx(hDiv), fill cor Line 295  C        to compute d/dx(hDiv), fill cor
295           ENDDO           ENDDO
296    
297  C-       gradient in y direction:  C-       gradient in y direction:
298  #ifndef ALLOW_AUTODIFF_TAMC  cph-exch2#ifndef ALLOW_AUTODIFF_TAMC
299           IF (useCubedSphereExchange) THEN           IF (useCubedSphereExchange) THEN
300  C        to compute d/dy(hDiv), fill corners with appropriate values:  C        to compute d/dy(hDiv), fill corners with appropriate values:
301             CALL FILL_CS_CORNER_TR_RL(.FALSE., hDiv, bi,bj, myThid )             CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
302         &                                hDiv, bi,bj, myThid )
303           ENDIF           ENDIF
304  #endif  cph-exch2#endif
305           DO j=2-Oly,sNy+Oly-1           DO j=2-Oly,sNy+Oly-1
306            DO i=2-Olx,sNx+Olx-1            DO i=2-Olx,sNx+Olx-1
307              divDy(i,j) = (hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj)              divDy(i,j) = (hDiv(i,j)-hDiv(i,j-1))*recip_DYC(i,j,bi,bj)
308            ENDDO            ENDDO
309           ENDDO           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         ENDIF
330    
331         DO j=2-Oly,sNy+Oly-1         DO j=2-Oly,sNy+Oly-1
332          DO i=2-Olx,sNx+Olx-1          DO i=2-Olx,sNx+Olx-1
333  CCCCCCCCCCCCCCC Divergence Point CalculationsCCCCCCCCCCCCCCCCCCCC  CCCCCCCCCCCCCCC Divergence Point CalculationsCCCCCCCCCCCCCCCCCCCC
334    
335  C These are (powers of) length scales  #ifdef    ALLOW_AUTODIFF_TAMC
336    # ifndef AUTODIFF_DISABLE_LEITH
337               lockey_2 = i+olx + (sNx+2*olx)*(j+oly-1)
338         &                      + (sNx+2*olx)*(sNy+2*oly)*(lockey_1-1)
339    CADJ STORE viscA4_ZSmg(i,j)
340    CADJ &     = comlev1_mom_ijk_loop , key=lockey_2, byte=isbyte
341    CADJ STORE viscAh_ZSmg(i,j)
342    CADJ &     = comlev1_mom_ijk_loop , key=lockey_2, byte=isbyte
343    # endif
344    #endif    /* ALLOW_AUTODIFF_TAMC */
345    
346    C These are (powers of) length scales
347           IF (useAreaViscLength) THEN           IF (useAreaViscLength) THEN
348            L2=rA(i,j,bi,bj)            L2=rA(i,j,bi,bj)
349              L4rdt=0.03125 _d 0*recip_dt*L2**2
350           ELSE           ELSE
351            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))
352              L4rdt=recip_dt/( 6. _d 0*(recip_DXF(I,J,bi,bj)**4
353         &                             +recip_DYF(I,J,bi,bj)**4)
354         &                   +8. _d 0*((recip_DXF(I,J,bi,bj)
355         &                             *recip_DYF(I,J,bi,bj))**2) )
356           ENDIF           ENDIF
357           L3=(L2**1.5)           L3=(L2**1.5)
358           L4=(L2**2)           L4=(L2**2)
359           L5=(L2**2.5)           L5=(L2*L3)
360    
361           L2rdt=0.25 _d 0*recip_dt*L2           L2rdt=0.25 _d 0*recip_dt*L2
362    
          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  
   
363  C Velocity Reynolds Scale  C Velocity Reynolds Scale
364           IF ( viscAhRe_max.GT.0. .AND. KE(i,j).GT.0. ) THEN           IF ( viscAhRe_max.GT.0. .AND. KE(i,j).GT.0. ) THEN
365             Uscl=sqrt(KE(i,j)*L2)*viscAhRe_max             Uscl=SQRT(KE(i,j)*L2)*viscAhRe_max
366           ELSE           ELSE
367             Uscl=0.             Uscl=0.
368           ENDIF           ENDIF
369           IF ( viscA4Re_max.GT.0. .AND. KE(i,j).GT.0. ) THEN           IF ( viscA4Re_max.GT.0. .AND. KE(i,j).GT.0. ) THEN
370             U4scl=sqrt(KE(i,j))*L3*viscA4Re_max             U4scl=SQRT(KE(i,j))*L3*viscA4Re_max
371           ELSE           ELSE
372             U4scl=0.             U4scl=0.
373           ENDIF           ENDIF
374    
375           IF (useFullLeith.and.calcleith) THEN  cph-leith#ifndef ALLOW_AUTODIFF_TAMC
376    #ifndef AUTODIFF_DISABLE_LEITH
377             IF (useFullLeith.AND.calcLeith) THEN
378  C This is the vector magnitude of the vorticity gradient squared  C This is the vector magnitude of the vorticity gradient squared
379            grdVrt=0.25 _d 0*(            grdVrt=0.25 _d 0*( (vrtDx(i,j+1)*vrtDx(i,j+1)
380       &     ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2       &                        + vrtDx(i,j)*vrtDx(i,j) )
381       &     +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2       &                     + (vrtDy(i+1,j)*vrtDy(i+1,j)
382       &     +((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)  
383    
384  C This is the vector magnitude of grad (div.v) squared  C This is the vector magnitude of grad (div.v) squared
385  C Using it in Leith serves to damp instabilities in w.  C Using it in Leith serves to damp instabilities in w.
# Line 297  C Using it in Leith serves to damp insta Line 389  C Using it in Leith serves to damp insta
389       &                        + divDy(i,j)*divDy(i,j) )  )       &                        + divDy(i,j)*divDy(i,j) )  )
390    
391            viscAh_DLth(i,j)=            viscAh_DLth(i,j)=
392       &     sqrt(leith2fac*grdVrt+leithD2fac*grdDiv)*L3       &     SQRT(leith2fac*grdVrt+leithD2fac*grdDiv)*L3
393            viscA4_DLth(i,j)=            viscA4_DLth(i,j)=
394       &     sqrt(leith4fac*grdVrt+leithD4fac*grdDiv)*L5       &     SQRT(leith4fac*grdVrt+leithD4fac*grdDiv)*L5
395            viscAh_DLthd(i,j)=            viscAh_DLthd(i,j)=
396       &     sqrt(leithD2fac*grdDiv)*L3       &     SQRT(leithD2fac*grdDiv)*L3
397            viscA4_DLthd(i,j)=            viscA4_DLthd(i,j)=
398       &     sqrt(leithD4fac*grdDiv)*L5       &     SQRT(leithD4fac*grdDiv)*L5
399           ELSEIF (calcleith) THEN           ELSEIF (calcLeith) THEN
400  C but this approximation will work on cube  C but this approximation will work on cube
401  c (and differs by as much as 4X)  c (and differs by as much as 4X)
402            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)) )
403            grdVrt=max(grdVrt,            grdVrt=MAX( grdVrt, ABS(vrtDy(i+1,j)) )
404       &     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=max( abs(divDx(i+1,j)), abs(divDx(i,j)) )  
           grdDiv=max( grdDiv, abs(divDy(i,j+1)) )  
           grdDiv=max( grdDiv, abs(divDy(i,j))   )  
405    
406  c This approximation is good to the same order as above...  c This approximation is good to the same order as above...
407              grdDiv=MAX( ABS(divDx(i+1,j)), ABS(divDx(i,j)) )
408              grdDiv=MAX( grdDiv, ABS(divDy(i,j+1)) )
409              grdDiv=MAX( grdDiv, ABS(divDy(i,j))   )
410    
411            viscAh_Dlth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3            viscAh_Dlth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3
412            viscA4_Dlth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5            viscA4_Dlth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5
413            viscAh_DlthD(i,j)=((leithD2fac*grdDiv))*L3            viscAh_DlthD(i,j)=((leithD2fac*grdDiv))*L3
# Line 331  c This approximation is good to the same Line 419  c This approximation is good to the same
419            viscA4_DlthD(i,j)=0. _d 0            viscA4_DlthD(i,j)=0. _d 0
420           ENDIF           ENDIF
421    
422           IF (calcsmag) THEN           IF (calcSmag) THEN
423            viscAh_DSmg(i,j)=L2            viscAh_DSmg(i,j)=L2
424       &       *sqrt(tension(i,j)**2       &       *SQRT(tension(i,j)**2
425       &       +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
426       &                  +strain(i  , j )**2+strain(i+1,j+1)**2))       &                  +strain(i  , j )**2+strain(i+1,j+1)**2))
427            viscA4_DSmg(i,j)=smag4fac*L2*viscAh_DSmg(i,j)            viscA4_DSmg(i,j)=smag4fac*L2*viscAh_DSmg(i,j)
# Line 342  c This approximation is good to the same Line 430  c This approximation is good to the same
430            viscAh_DSmg(i,j)=0. _d 0            viscAh_DSmg(i,j)=0. _d 0
431            viscA4_DSmg(i,j)=0. _d 0            viscA4_DSmg(i,j)=0. _d 0
432           ENDIF           ENDIF
433    #endif /* AUTODIFF_DISABLE_LEITH */
434    
435  C  Harmonic on Div.u points  C  Harmonic on Div.u points
436           Alin=viscAhD+viscAhGrid*L2rdt           Alin=viscAhD+viscAhGrid*L2rdt
437       &          +viscAh_DLth(i,j)+viscAh_DSmg(i,j)       &          +viscAh_DLth(i,j)+viscAh_DSmg(i,j)
438           viscAh_DMin(i,j)=max(viscAhGridMin*L2rdt,Uscl)           viscAh_DMin(i,j)=MAX(viscAhGridMin*L2rdt,Uscl)
439           viscAh_D(i,j)=max(viscAh_DMin(i,j),Alin)           viscAh_D(i,j)=MAX(viscAh_DMin(i,j),Alin)
440           viscAh_DMax(i,j)=min(viscAhGridMax*L2rdt,viscAhMax)           viscAh_DMax(i,j)=MIN(viscAhGridMax*L2rdt,viscAhMax)
441           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))
442    
443  C  BiHarmonic on Div.u points  C  BiHarmonic on Div.u points
444           Alin=viscA4D+viscA4Grid*L4rdt           Alin=viscA4D+viscA4Grid*L4rdt
445       &          +viscA4_DLth(i,j)+viscA4_DSmg(i,j)       &          +viscA4_DLth(i,j)+viscA4_DSmg(i,j)
446           viscA4_DMin(i,j)=max(viscA4GridMin*L4rdt,U4scl)           viscA4_DMin(i,j)=MAX(viscA4GridMin*L4rdt,U4scl)
447           viscA4_D(i,j)=max(viscA4_DMin(i,j),Alin)           viscA4_D(i,j)=MAX(viscA4_DMin(i,j),Alin)
448           viscA4_DMax(i,j)=min(viscA4GridMax*L4rdt,viscA4Max)           viscA4_DMax(i,j)=MIN(viscA4GridMax*L4rdt,viscA4Max)
449           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))
450    
451    #ifdef ALLOW_NONHYDROSTATIC
452    C--   Pass Viscosities to calc_gw, if constant, not necessary
453    
454             kp1  = MIN(k+1,Nr)
455    
456             IF ( k.EQ.1 ) THEN
457    C     Prepare for next level (next call)
458              viscAh_W(i,j,kp1,bi,bj)=0.5*viscAh_D(i,j)
459              viscA4_W(i,j,kp1,bi,bj)=0.5*viscA4_D(i,j)
460    
461    C     These values dont get used
462              viscAh_W(i,j,k,bi,bj)=viscAh_D(i,j)
463              viscA4_W(i,j,k,bi,bj)=viscA4_D(i,j)
464    
465             ELSEIF ( k.EQ.Nr ) THEN
466              viscAh_W(i,j,k,bi,bj)=viscAh_W(i,j,k,bi,bj)+0.5*viscAh_D(i,j)
467              viscA4_W(i,j,k,bi,bj)=viscA4_W(i,j,k,bi,bj)+0.5*viscA4_D(i,j)
468    
469             ELSE
470    C     Prepare for next level (next call)
471              viscAh_W(i,j,kp1,bi,bj)=0.5*viscAh_D(i,j)
472              viscA4_W(i,j,kp1,bi,bj)=0.5*viscA4_D(i,j)
473    
474    C     Note that previous call of this function has already added half.
475              viscAh_W(i,j,k,bi,bj)=viscAh_W(i,j,k,bi,bj)+0.5*viscAh_D(i,j)
476              viscA4_W(i,j,k,bi,bj)=viscA4_W(i,j,k,bi,bj)+0.5*viscA4_D(i,j)
477    
478             ENDIF
479    #endif /* ALLOW_NONHYDROSTATIC */
480    
481  CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC  CCCCCCCCCCCCC Vorticity Point CalculationsCCCCCCCCCCCCCCCCCC
482  C These are (powers of) length scales  C These are (powers of) length scales
483           IF (useAreaViscLength) THEN           IF (useAreaViscLength) THEN
484            L2=rAz(i,j,bi,bj)            L2=rAz(i,j,bi,bj)
485              L4rdt=0.125 _d 0*recip_dt*rAz(i,j,bi,bj)**2
486           ELSE           ELSE
487            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))
488              L4rdt=recip_dt/
489         &     ( 6. _d 0*(recip_DXV(I,J,bi,bj)**4+recip_DYU(I,J,bi,bj)**4)
490         &      +8. _d 0*((recip_DXV(I,J,bi,bj)*recip_DYU(I,J,bi,bj))**2))
491           ENDIF           ENDIF
492    
493           L3=(L2**1.5)           L3=(L2**1.5)
494           L4=(L2**2)           L4=(L2**2)
495           L5=(L2**2.5)           L5=(L2*L3)
496    
497           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  
498    
499  C Velocity Reynolds Scale (Pb here at CS-grid corners !)  C Velocity Reynolds Scale (Pb here at CS-grid corners !)
500           IF ( viscAhRe_max.GT.0. .OR. viscA4Re_max.GT.0. ) THEN           IF ( viscAhRe_max.GT.0. .OR. viscA4Re_max.GT.0. ) THEN
501             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))
502       &                      +(KE(i-1,j)+KE(i,j-1)) )       &                      +(KE(i-1,j)+KE(i,j-1)) )
503             IF ( keZpt.GT.0. ) THEN             IF ( keZpt.GT.0. ) THEN
504               Uscl = sqrt(keZpt*L2)*viscAhRe_max               Uscl = SQRT(keZpt*L2)*viscAhRe_max
505               U4scl= sqrt(keZpt)*L3*viscA4Re_max               U4scl= SQRT(keZpt)*L3*viscA4Re_max
506             ELSE             ELSE
507               Uscl =0.               Uscl =0.
508               U4scl=0.               U4scl=0.
# Line 396  C Velocity Reynolds Scale (Pb here at CS Line 512  C Velocity Reynolds Scale (Pb here at CS
512             U4scl=0.             U4scl=0.
513           ENDIF           ENDIF
514    
515    #ifndef AUTODIFF_DISABLE_LEITH
516  C This is the vector magnitude of the vorticity gradient squared  C This is the vector magnitude of the vorticity gradient squared
517           IF (useFullLeith.and.calcleith) THEN           IF (useFullLeith.AND.calcLeith) THEN
518            grdVrt=0.25 _d 0*(            grdVrt=0.25 _d 0*( (vrtDx(i-1,j)*vrtDx(i-1,j)
519       &     ((vort3(i+1,j)-vort3(i,j))*recip_DXG(i,j,bi,bj))**2       &                        + vrtDx(i,j)*vrtDx(i,j) )
520       &     +((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj))**2       &                     + (vrtDy(i,j-1)*vrtDy(i,j-1)
521       &     +((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)  
522    
523  C This is the vector magnitude of grad(div.v) squared  C This is the vector magnitude of grad(div.v) squared
524            grdDiv=0.25 _d 0*( (divDx(i,j-1)*divDx(i,j-1)            grdDiv=0.25 _d 0*( (divDx(i,j-1)*divDx(i,j-1)
# Line 411  C This is the vector magnitude of grad(d Line 527  C This is the vector magnitude of grad(d
527       &                        + divDy(i,j)*divDy(i,j) )  )       &                        + divDy(i,j)*divDy(i,j) )  )
528    
529            viscAh_ZLth(i,j)=            viscAh_ZLth(i,j)=
530       &     sqrt(leith2fac*grdVrt+leithD2fac*grdDiv)*L3       &     SQRT(leith2fac*grdVrt+leithD2fac*grdDiv)*L3
531            viscA4_ZLth(i,j)=            viscA4_ZLth(i,j)=
532       &     sqrt(leith4fac*grdVrt+leithD4fac*grdDiv)*L5       &     SQRT(leith4fac*grdVrt+leithD4fac*grdDiv)*L5
533            viscAh_ZLthD(i,j)=            viscAh_ZLthD(i,j)=
534       &     sqrt(leithD2fac*grdDiv)*L3       &     SQRT(leithD2fac*grdDiv)*L3
535            viscA4_ZLthD(i,j)=            viscA4_ZLthD(i,j)=
536       &     sqrt(leithD4fac*grdDiv)*L5       &     SQRT(leithD4fac*grdDiv)*L5
537    
538           ELSEIF (calcleith) THEN           ELSEIF (calcLeith) THEN
539  C but this approximation will work on cube (and differs by 4X)  C but this approximation will work on cube (and differs by 4X)
540            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)) )
541            grdVrt=max(grdVrt,            grdVrt=MAX( grdVrt, ABS(vrtDy(i,j-1)) )
542       &     abs((vort3(i,j+1)-vort3(i,j))*recip_DYG(i,j,bi,bj)))            grdVrt=MAX( grdVrt, ABS(vrtDy(i,j))   )
543            grdVrt=max(grdVrt,  
544       &     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)) )
545            grdVrt=max(grdVrt,            grdDiv=MAX( grdDiv, ABS(divDy(i,j))   )
546       &     abs((vort3(i,j-1)-vort3(i,j))*recip_DYG(i,j-1,bi,bj)))            grdDiv=MAX( grdDiv, ABS(divDy(i-1,j)) )
   
           grdDiv=max( abs(divDx(i,j)), abs(divDx(i,j-1)) )  
           grdDiv=max( grdDiv, abs(divDy(i,j))   )  
           grdDiv=max( grdDiv, abs(divDy(i-1,j)) )  
547    
548            viscAh_ZLth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3            viscAh_ZLth(i,j)=(leith2fac*grdVrt+(leithD2fac*grdDiv))*L3
549            viscA4_ZLth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5            viscA4_ZLth(i,j)=(leith4fac*grdVrt+(leithD4fac*grdDiv))*L5
# Line 444  C but this approximation will work on cu Line 556  C but this approximation will work on cu
556            viscA4_ZLthD(i,j)=0. _d 0            viscA4_ZLthD(i,j)=0. _d 0
557           ENDIF           ENDIF
558    
559           IF (calcsmag) THEN           IF (calcSmag) THEN
560            viscAh_ZSmg(i,j)=L2            viscAh_ZSmg(i,j)=L2
561       &      *sqrt(strain(i,j)**2       &      *SQRT(strain(i,j)**2
562       &        +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
563       &                   +tension(i-1, j )**2+tension(i-1,j-1)**2))       &                   +tension(i-1, j )**2+tension(i-1,j-1)**2))
564            viscA4_ZSmg(i,j)=smag4fac*L2*viscAh_ZSmg(i,j)            viscA4_ZSmg(i,j)=smag4fac*L2*viscAh_ZSmg(i,j)
565            viscAh_ZSmg(i,j)=smag2fac*viscAh_ZSmg(i,j)            viscAh_ZSmg(i,j)=smag2fac*viscAh_ZSmg(i,j)
566           ENDIF           ENDIF
567    #endif /* AUTODIFF_DISABLE_LEITH */
568    
569  C  Harmonic on Zeta points  C  Harmonic on Zeta points
570           Alin=viscAhZ+viscAhGrid*L2rdt           Alin=viscAhZ+viscAhGrid*L2rdt
571       &           +viscAh_ZLth(i,j)+viscAh_ZSmg(i,j)       &           +viscAh_ZLth(i,j)+viscAh_ZSmg(i,j)
572           viscAh_ZMin(i,j)=max(viscAhGridMin*L2rdt,Uscl)           viscAh_ZMin(i,j)=MAX(viscAhGridMin*L2rdt,Uscl)
573           viscAh_Z(i,j)=max(viscAh_ZMin(i,j),Alin)           viscAh_Z(i,j)=MAX(viscAh_ZMin(i,j),Alin)
574           viscAh_ZMax(i,j)=min(viscAhGridMax*L2rdt,viscAhMax)           viscAh_ZMax(i,j)=MIN(viscAhGridMax*L2rdt,viscAhMax)
575           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))
576    
577  C  BiHarmonic on Zeta points  C  BiHarmonic on Zeta points
578           Alin=viscA4Z+viscA4Grid*L4rdt           Alin=viscA4Z+viscA4Grid*L4rdt
579       &           +viscA4_ZLth(i,j)+viscA4_ZSmg(i,j)       &           +viscA4_ZLth(i,j)+viscA4_ZSmg(i,j)
580           viscA4_ZMin(i,j)=max(viscA4GridMin*L4rdt,U4scl)           viscA4_ZMin(i,j)=MAX(viscA4GridMin*L4rdt,U4scl)
581           viscA4_Z(i,j)=max(viscA4_ZMin(i,j),Alin)           viscA4_Z(i,j)=MAX(viscA4_ZMin(i,j),Alin)
582           viscA4_ZMax(i,j)=min(viscA4GridMax*L4rdt,viscA4Max)           viscA4_ZMax(i,j)=MIN(viscA4GridMax*L4rdt,viscA4Max)
583           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))
584          ENDDO          ENDDO
585         ENDDO         ENDDO
586    
587        ELSE        ELSE
588    C---- use constant viscosity (useVariableViscosity=F):
589    
590         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
591          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
592           viscAh_D(i,j)=viscAhD           viscAh_D(i,j)=viscAhD
# Line 479  C  BiHarmonic on Zeta points Line 595  C  BiHarmonic on Zeta points
595           viscA4_Z(i,j)=viscA4Z           viscA4_Z(i,j)=viscA4Z
596          ENDDO          ENDDO
597         ENDDO         ENDDO
598    
599    C---- variable/constant viscosity : end if/else block
600        ENDIF        ENDIF
601    
602  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
Legend:
Removed from v.1.19  
changed lines
  Added in v.1.37
  ViewVC Help
Powered by ViewVC 1.1.22